Skip to content

fix(zratio): truncation guards, ties-mean slope, cache and build wiring #207

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
MaartenMarsman merged 1 commit into from
Oct 21, 2026
Merged

fix(zratio): truncation guards, ties-mean slope, cache and build wiring #207

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
129 changes: 85 additions & 44 deletions R/correction_tables.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,7 @@
# so the per-pair curve f(theta) = log C(theta) / E is the cumulative
# trapezoid of the per-pair score, up to an additive constant that cancels in
# every ratio. The implied per-edge slope at local density d is
# f'(d) = log[(d / (1 - d)) * (1 - theta) / theta], and the per-edge read-off
# table is fed(theta, d) = log(1 - theta + theta * exp(f'(d))).
# f'(d) = log[(d / (1 - d)) * (1 - theta) / theta].
#
# The sweep runs bgms's own tilted prior sampler (sample_ggm_prior,
# spec = "joint"), so the tables are consistent with the deployed model by
Expand Down Expand Up @@ -69,18 +68,17 @@ ggm_correction_theta_grid = function(n_grid = 120L, lower = 0.002,
# increasing in theta), then clamp away from {0, 1} for the logs
# - per-pair score -> trapezoid -> f(theta) (per-pair log C)
# - implied per-edge slope fprime vs local density
# - fed(theta, d) read-off grid
#
# The integrated f/logC curve uses the full theta grid (the score has
# no boundary singularity in the estimate). The slope curve is
# fprime = logit(d) - logit(theta), so near the boundaries the logit
# amplifies Monte-Carlo noise in d and the resolution floor 0.5/E can
# exceed the true density; fprime and the fed table are therefore
# built only from grid points with theta <= fprime_theta_cap whose
# density is resolvable (inside (0.5/E, 1 - 0.5/E)), and extended as
# constants beyond the covered density range. With a single pair the
# resolvable window is empty, so the slope pieces are returned as
# NULL; the logC curve does not depend on them.
# exceed the true density; fprime is therefore built only from grid
# points with theta <= fprime_theta_cap whose density is resolvable
# (inside (0.5/E, 1 - 0.5/E)), and extended as constants beyond the
# covered density range. With a single pair the resolvable window is
# empty, so the slope pieces are returned as NULL; the logC curve does
# not depend on them.
#
# Returns the table as a list; logC is the whole-graph curve
# num_pairs * f used by the beta-binomial draw. fprime is tabulated
Expand All @@ -102,26 +100,23 @@ correction_table_from_edens = function(theta, edens_raw, num_pairs,
edens > 0.5 / num_pairs & edens < 1 - 0.5 / num_pairs
fprime_density = NULL
fprime = NULL
fed_theta = NULL
fed_density = NULL
fed = NULL
if(any(keep)) {
fprime_density = edens[keep]
fprime = log(
(fprime_density / (1 - fprime_density)) * (1 - theta[keep]) / theta[keep]
dens_keep = edens[keep]
fprime_keep = log(
(dens_keep / (1 - dens_keep)) * (1 - theta[keep]) / theta[keep]
)
stopifnot(all(is.finite(fprime)))

fed_theta = seq(0.005, 0.995, length.out = 120L)
fed_density = seq(min(fprime_density), max(fprime_density),
length.out = 60L
)
fprime_at = function(d) {
stats::approx(fprime_density, fprime, d, rule = 2, ties = mean)$y
}
fed = outer(fed_theta, fed_density, function(tt, dd) {
log(1 - tt + tt * exp(fprime_at(dd)))
})
stopifnot(all(is.finite(fprime_keep)))

# The isotonic repair pools a crash dip into a run of grid points with one
# fitted density while theta keeps rising, so fprime falls along the run.
# The C++ reader (SBMCorrection::fprime_at) brackets by binary search and
# takes the last point of a tied run. Aggregate to unique densities, mean
# slope per run, so the abscissae are strictly increasing. isoreg emits
# bit-identical values within a run, so exact-value grouping is safe.
grp = match(dens_keep, unique(dens_keep))
fprime_density = unique(dens_keep)
fprime = as.numeric(tapply(fprime_keep, grp, mean))
stopifnot(all(diff(fprime_density) > 0))
}

list(
Expand All @@ -134,9 +129,6 @@ correction_table_from_edens = function(theta, edens_raw, num_pairs,
fprime_density = fprime_density,
fprime = fprime,
fprime_theta_cap = fprime_theta_cap,
fed_theta = fed_theta,
fed_density = fed_density,
fed = fed,
num_pairs = num_pairs
)
}
Expand Down Expand Up @@ -447,6 +439,22 @@ correction_scale_prior = function(prior) {
}


# ------------------------------------------------------------------
# warn_sbm_slope_unresolvable (internal)
# ------------------------------------------------------------------
# The single wording for "this cell has no slope curve", raised by the
# pre-sweep check and by the post-build check.
# ------------------------------------------------------------------
warn_sbm_slope_unresolvable = function() {
warning(
"The Stochastic-Block updates are run without the ",
"normalizing-constant correction: the slope curve is not resolvable ",
"for this model cell (a single tilted pair).",
call. = FALSE
)
}


# ------------------------------------------------------------------
# ggm_edge_prior_correction (internal)
# ------------------------------------------------------------------
Expand Down Expand Up @@ -488,6 +496,15 @@ ggm_edge_prior_correction = function(prior, sampler, num_variables,
return(NULL)
}

# With two continuous variables the tilted subgraph has one pair, so the
# resolvable density window (0.5/E, 1 - 0.5/E) is empty and the slope curve
# is NULL for certain (correction_table_from_edens). Warn before the sweep;
# the post-build check below covers the general case.
if(identical(prior$edge_prior, "Stochastic-Block") && num_continuous < 3L) {
warn_sbm_slope_unresolvable()
return(NULL)
}

interaction_prior = correction_interaction_prior(prior)
precision_scale_prior = correction_scale_prior(prior)

Expand All @@ -497,23 +514,21 @@ ggm_edge_prior_correction = function(prior, sampler, num_variables,
show_progress = is.null(sampler$progress_type) ||
!identical(as.integer(sampler$progress_type), 0L)

# The table build runs before the chains launch, on an idle machine, so it
# sizes its own core count and announces itself from the verbose option
# (same policy as prior_pip_table()).
table = ggm_correction_table(
p = num_continuous, delta = prior$delta,
interaction_prior = interaction_prior,
precision_scale_prior = precision_scale_prior,
update_method = "gibbs",
cores = sampler$cores,
verbose = isTRUE(sampler$verbose),
cores = normalize_parallel_cores(parallel::detectCores()),
verbose = isTRUE(getOption("bgms.verbose", TRUE)),
show_progress = show_progress
)
if(identical(prior$edge_prior, "Stochastic-Block") &&
is.null(table$fprime)) {
warning(
"The Stochastic-Block updates are run without the ",
"normalizing-constant correction: the slope curve is not resolvable ",
"for this model cell (a single tilted pair).",
call. = FALSE
)
warn_sbm_slope_unresolvable()
return(NULL)
}
is_continuous = NULL
Expand Down Expand Up @@ -553,18 +568,26 @@ correction_cache_dir = function() {
# earlier version's table out of the shared cache directory. Files
# under an older key are simply never read again.
#
# base_seed seeds every chain in the sweep, so it is part of the key.
# Prefix v2: v1 tables have tied fprime abscissae and are not served to
# this reader.
#
# @param cell Cell identity from ggm_correction_cell().
# @param n_grid,n_samples,n_warmup,n_seeds,update_method Builder settings.
# @param base_seed Seed offset for the sweep chains.
#
# Returns: the file name, including the .rds extension.
# ------------------------------------------------------------------
ggm_correction_table_key = function(cell, n_grid, n_samples, n_warmup,
n_seeds, update_method) {
n_seeds, update_method, base_seed) {
sprintf(
"ggm_ctable_v1_%s_q%d_delta%.8g_eta%.8g_%s_shape%.8g_g%d_ns%d_nw%d_sd%d_%s.rds",
paste0(
"ggm_ctable_v2_%s_q%d_delta%.8g_eta%.8g_%s_shape%.8g",
"_g%d_ns%d_nw%d_sd%d_bs%d_%s.rds"
),
as.character(utils::packageVersion("bgms")),
cell$q, cell$delta, cell$eta, cell$slab_family, cell$scale_shape,
n_grid, n_samples, n_warmup, n_seeds, update_method
n_grid, n_samples, n_warmup, n_seeds, base_seed, update_method
)
}

Expand Down Expand Up @@ -598,7 +621,7 @@ ggm_correction_table = function(
p, delta, interaction_prior, precision_scale_prior
)
key = ggm_correction_table_key(
cell, n_grid, n_samples, n_warmup, n_seeds, update_method
cell, n_grid, n_samples, n_warmup, n_seeds, update_method, base_seed
)
cache_dir = correction_cache_dir()
cache_file = file.path(cache_dir, key)
Expand Down Expand Up @@ -626,8 +649,26 @@ ggm_correction_table = function(
)

if(use_cache) {
dir.create(dirname(cache_file), recursive = TRUE, showWarnings = FALSE)
saveRDS(table, cache_file)
# The read side tolerates an unusable cache; so does the write side. A
# failed write is one warning saying the table could not be cached.
# saveRDS warns ("cannot open compressed file") before it errors, so both
# conditions are caught.
report_failure = function(cond) {
warning(
"The correction table was built but could not be cached to ",
cache_file, " (", conditionMessage(cond),
"); later fits of this model cell will rebuild it.",
call. = FALSE
)
}
tryCatch(
{
dir.create(dirname(cache_file), recursive = TRUE, showWarnings = FALSE)
saveRDS(table, cache_file)
},
error = report_failure,
warning = report_failure
)
}
table
}
27 changes: 21 additions & 6 deletions R/zratio_surfaces.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -373,8 +373,14 @@ zratio_anchor_grids_empty = function(cap) {
# the cap cannot drift between them.
.zratio_surface_size_cap = 80L

# `verbose` governs the build announcement and the completion line. It
# defaults to the bgms.verbose option; zratio_attach_surface passes the fit's
# own flag.
zratio_build_surfaces = function(zc, max_size = .zratio_surface_size_cap,
cores = 1L, seed0 = 700000L) {
cores = 1L, seed0 = 700000L,
verbose = isTRUE(
getOption("bgms.verbose", TRUE)
)) {
if(zc$alpha < .zratio_surface_shape_lo ||
zc$alpha > .zratio_surface_shape_hi) {
return(NULL)
Expand Down Expand Up @@ -479,7 +485,7 @@ zratio_build_surfaces = function(zc, max_size = .zratio_surface_size_cap,
# the one-time cost so the pre-chain pause is not silent. n_workers is the
# parallelism actually used, which drops to 1 on the Windows serial tier.
n_workers = if(use_psock || .Platform$OS.type == "unix") cores else 1L
verbose = isTRUE(getOption("bgms.verbose", TRUE))
verbose = isTRUE(verbose)
if(verbose) {
message(
"Building normalizing-constant corrections for the hierarchical ",
Expand Down Expand Up @@ -511,8 +517,13 @@ zratio_build_surfaces = function(zc, max_size = .zratio_surface_size_cap,
mc.cores = mc, mc.preschedule = FALSE
)
}
cn_rows = do.call(rbind, res[jobs$fam == "cn"])
bip_rows = do.call(rbind, res[jobs$fam == "bip"])
# An anchor job returns a one-row data frame, or NULL when the oracle did
# not converge; rbind drops the NULLs. A worker that died comes back from
# mclapply as a try-error, which rbind would coerce into a character matrix,
# so keep only the data frames.
is_row = vapply(res, is.data.frame, logical(1))
cn_rows = do.call(rbind, res[is_row & jobs$fam == "cn"])
bip_rows = do.call(rbind, res[is_row & jobs$fam == "bip"])

if(is.null(cn_rows) || is.null(bip_rows)) {
return(NULL)
Expand Down Expand Up @@ -656,13 +667,17 @@ zratio_spec_list = function(zc, gauge_sweeps) {
# Build the Option-B surface for cell `zc`, sized on `size` variables, and, on a
# successful build, attach it to the `zratio` spec; otherwise keep whatever
# route the spec already carries -- isolated-edge past the validated shape
# range, additive elsewhere -- and message the reason when `verbose`.
# range, additive elsewhere -- and message the reason when `verbose`. The same
# flag governs the build's own announcement.
# Centralizes the size cap, cores policy, and fence message the GGM, mixed, and
# prior sampler paths share. Returns the (possibly surface-carrying) `zratio`
# list.
zratio_attach_surface = function(zratio, zc, size, cores, verbose = FALSE) {
max_size = min(size, .zratio_surface_size_cap)
surf = zratio_build_surfaces(zc, max_size = max_size, cores = cores)
surf = zratio_build_surfaces(
zc,
max_size = max_size, cores = cores, verbose = verbose
)
if(!is.null(surf)) {
zratio$surface = surf
} else if(isTRUE(verbose)) {
Expand Down
76 changes: 51 additions & 25 deletions R/zratio_tables.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -454,6 +454,23 @@ zratio_eta = function(pairwise_scale, scale_rate, scale_eta = NA_real_) {
.zratio_constants_shape_lo = 0.5
.zratio_constants_shape_hi = 20

# Retained-edge-mass tolerances for the truncation guard below, one per
# channel. Each is a ratio of a pair integral at the end of the fixed c-grid to
# its value at the start, so it measures how much mass the grid fails to cover.
#
# I_spike carries no slab dependence -- the spike channel never sees the slab
# -- and its decay tracks the grid width directly: 3.0e-01 at eta = 0.05,
# 2.4e-02 at eta = 0.1, 1.2e-06 at eta = 0.25, below 1e-14 past eta = 0.5, the
# same figures under either family. One tolerance serves it.
#
# G depends on the slab. Under the normal slab it decays as sharply (1e-34 at
# eta = 1). Under the Cauchy slab the slab's polynomial tail leaves a floor of
# 5.7e-04 to 3.1e-03 across the certified band (eta >= 0.25, delta 0.1 to 3)
# that no grid width removes, so its tolerance sits above that floor; width
# is still detected there through I_spike.
.zratio_truncation_tol_spike = 1e-6
.zratio_truncation_tol_g = c(normal = 1e-6, cauchy = 1e-2)

# Session cache for zratio_constants: the constant set is deterministic per
# (delta, eta, alpha, slab) cell, and one fit resolves the same cell more
# than once (sampler dispatch and diagnostics assembly).
Expand Down Expand Up @@ -495,31 +512,40 @@ zratio_constants = function(delta, eta, alpha = 1, slab = "normal") {
)
}
pair = zratio_pair_integrals(delta, sigma, beta, slab, alpha)
if(abs(alpha - 1) > 1e-12) {
if(alpha < .zratio_constants_shape_lo ||
alpha > .zratio_constants_shape_hi) {
warning(
"Z-ratio constants at diagonal shape ", format(alpha),
" lie outside the certified range [",
format(.zratio_constants_shape_lo), ", ",
format(.zratio_constants_shape_hi),
"]; the fixed quadrature grids are not scored there and may lose ",
"accuracy.",
call. = FALSE
)
}
tail_g = pair$gv[length(pair$gv)] / pair$gv[1]
tail_i = pair$ispike(max(pair$cg)) / pair$ispike(0)
if(!is.finite(tail_g) || !is.finite(tail_i) ||
tail_g > 1e-6 || tail_i > 1e-6) {
warning(
"Z-ratio pair integrals retain visible mass at the grid edge ",
"(delta = ", format(delta), ", eta = ", format(eta),
", shape = ", format(alpha), "); ",
"the saddle tables may be truncated.",
call. = FALSE
)
}
# The certified-range warning is shape-specific: the grids were scored at
# alpha = 1 by construction, so only a non-exponential diagonal can leave the
# certified band.
if(abs(alpha - 1) > 1e-12 &&
(alpha < .zratio_constants_shape_lo ||
alpha > .zratio_constants_shape_hi)) {
warning(
"Z-ratio constants at diagonal shape ", format(alpha),
" lie outside the certified range [",
format(.zratio_constants_shape_lo), ", ",
format(.zratio_constants_shape_hi),
"]; the fixed quadrature grids are not scored there and may lose ",
"accuracy.",
call. = FALSE
)
}
# The truncation check runs for every cell. It is the only detector that the
# fixed c-grid (cmax in zratio_pair_integrals) is wide enough for this
# (delta, eta, alpha): the pair tables clamp to their edge value when read
# past it, so retained edge mass means the saddle tables integrate a plateau
# instead of a decaying tail. The two channels have their own tolerances;
# see .zratio_truncation_tol_g above.
tail_g = pair$gv[length(pair$gv)] / pair$gv[1]
tail_i = pair$ispike(max(pair$cg)) / pair$ispike(0)
tol_g = .zratio_truncation_tol_g[[slab]]
if(!is.finite(tail_g) || !is.finite(tail_i) ||
tail_g > tol_g || tail_i > .zratio_truncation_tol_spike) {
warning(
"Z-ratio pair integrals retain visible mass at the grid edge ",
"(delta = ", format(delta), ", eta = ", format(eta),
", shape = ", format(alpha), "); ",
"the saddle tables may be truncated.",
call. = FALSE
)
}
grid = zratio_saddle_grid(pair)
w12 = zratio_node_channel(delta, sigma, beta, slab, alpha)
Expand Down
Loading
Loading