Graph-level prior samplers and a progress bar for the correction-table build

DESCRIPTION

@@ -28,14 +28,15 @@ Copyright: Includes datasets 'ADHD' and 'Boredom', which are licensed under CC-B
 License: GPL (>= 2)
 URL: https://Bayesian-Graphical-Modelling-Lab.github.io/bgms/, https://github.com/Bayesian-Graphical-Modelling-Lab/bgms
 BugReports: https://github.com/Bayesian-Graphical-Modelling-Lab/bgms/issues
-Imports: 
-    Rcpp (>= 1.0.7), 
-    RcppParallel, 
-    Rdpack, 
+Imports:
+    Rcpp (>= 1.0.7),
+    RcppParallel,
+    Rdpack,
     S7,
-    methods, 
+    methods,
     lifecycle,
-    stats
+    stats,
+    utils
 RdMacros: Rdpack
 LinkingTo: 
     Rcpp, 

NAMESPACE

@@ -73,6 +73,8 @@ export(gamma_prior)
 export(mrfSampler)
 export(normal_prior)
 export(sample_ggm_prior)
+export(sample_graph_prior)
+export(sample_sbm_prior)
 export(sbm_prior)
 export(simulate_mrf)
 export(summarize_zratio_diagnostics)

NEWS.md

@@ -12,6 +12,9 @@
 
 ## New features
 
+* `sample_graph_prior()`: draws edge-inclusion indicators, together with any edge-prior hyperparameters, from the graph level of the spike-and-slab prior. Under the hierarchical specification the draw is ancestral and exact; under the joint specification it runs the zero-data prior chain, so the draws carry the per-graph normalizer tilt. Hyperparameters can be fixed instead of sampled (`theta` for Bernoulli and Beta-Bernoulli priors, `allocations` plus `block_probs` for the Stochastic-Block prior).
+* `sample_sbm_prior()`: ancestral draws of block allocations and pair-inclusion probabilities from the MFM-SBM edge-prior hyperprior.
+* The edge-prior correction table now announces itself only when it actually builds (a cache hit is silent) and states that the build is one-time and cached. In interactive sessions it draws a progress bar for both serial and parallel builds (a parallel build advances the bar as each batch of forked chains completes). The bar follows `display_progress` (the same control as the sampler's bar), not the advisory `bgms.verbose` flag.
 * Gaussian graphical models (GGM): `bgm(x, variable_type = "continuous")` fits a GGM with Bayesian edge selection. Sampling uses NUTS on a free-element Cholesky (theta-space) parameterization of the precision matrix, which keeps the precision matrix positive-definite by construction; adaptive-metropolis is also available.
 * GGM Gibbs sampler: `update_method = "gibbs"` fits a GGM with a conjugate row-by-row update of the precision matrix, with or without edge selection. It needs no step-size or proposal tuning and supports a Normal or Cauchy prior on the edges. Continuous data only.
 * Gibbs warmup staging: with `update_method = "gibbs"` and edge selection, the first 15% of the warmup runs the full model so the precision matrix settles, and edge selection is active for the remaining 85%. Previously edge selection only started at the first retained iteration, so the graph's equilibration happened inside the retained samples. Both windows scale with the warmup budget; the warmup default is unchanged.

R/bgm.R

@@ -286,7 +286,9 @@
 #' @param chains Integer. Number of parallel chains to run. Default: \code{4}.
 #'
 #' @param cores Integer. Number of CPU cores for parallel execution.
-#'   Default: \code{parallel::detectCores()}.
+#'   Sampling uses \code{min(cores, chains)}; some computations outside of
+#'   sampling (such as building the edge-selection prior correction table)
+#'   use all \code{cores}. Default: \code{parallel::detectCores()}.
 #'
 #' @param seed Optional integer. Random seed for reproducibility. Must be a
 #'   single non-negative integer.

R/correction_tables.R

@@ -153,6 +153,92 @@ normalize_builder_cores = function(cores) {
 }
 
 
+# ------------------------------------------------------------------
+# new_correction_progress (internal)
+# ------------------------------------------------------------------
+# In-place progress bar matching the bgms MCMC bar (see the theme and
+# bar-width logic in src/utils/progress_manager.cpp): tortoise-shell
+# brackets, a heavy horizontal rule filled in blue and empty in gray,
+# a sub-cell partial glyph, and a "prefix: <bar> cur/tot (xx.x%)"
+# layout redrawn with a carriage return. Unicode + ANSI colour when the
+# session is UTF-8, an ASCII "[=== ]" fallback otherwise. Returns
+# update(current) and close(); the caller decides whether to draw it.
+# ------------------------------------------------------------------
+new_correction_progress = function(total, prefix = "Correction table") {
+  unicode = isTRUE(l10n_info()[["UTF-8"]])
+  is_rstudio = Sys.getenv("RSTUDIO") == "1"
+
+  # Bar width, mirroring progress_manager.cpp so the bar does not wrap.
+  console_width = if(is_rstudio) {
+    max(0L, as.integer(getOption("width", 80L))) + 3L
+  } else {
+    80L
+  }
+  line_width = if(is_rstudio) {
+    max(10L, min(console_width - 25L, 70L))
+  } else {
+    70L
+  }
+  bar_width = if(line_width <= 5L) {
+    0L
+  } else if(line_width < 20L) {
+    line_width - 10L
+  } else if(line_width < 40L) {
+    line_width - 15L
+  } else if(line_width > 70L) {
+    40L
+  } else {
+    line_width - 30L
+  }
+  if(is_rstudio) {
+    bar_width = if(bar_width > 30L) bar_width - 20L else 10L
+  }
+
+  if(unicode) {
+    # Tortoise-shell brackets (U+2997/U+2998), heavy horizontal rule
+    # (U+2501) filled/empty, heavy sub-cell (U+257A); ANSI 38;5;73 blue
+    # and 37 gray, matching progress_manager.cpp.
+    lhs = "\u2997"
+    rhs = "\u2998"
+    filled = "\u001b[38;5;73m\u2501\u001b[39m"
+    partial_more = filled
+    partial_less = "\u001b[37m\u257a\u001b[39m"
+    empty = "\u001b[37m\u2501\u001b[39m"
+  } else {
+    lhs = "["
+    rhs = "]"
+    filled = "="
+    partial_more = " "
+    partial_less = " "
+    empty = " "
+  }
+
+  draw = function(current) {
+    frac = if(total > 0L) current / total else 1
+    exact = frac * bar_width
+    n_filled = min(as.integer(exact), bar_width)
+    bar = strrep(filled, n_filled)
+    if(n_filled < bar_width) {
+      part = exact - n_filled
+      if(part > 0) {
+        bar = paste0(bar, if(part > 0.5) partial_more else partial_less)
+        n_filled = n_filled + 1L
+      }
+    }
+    if(n_filled < bar_width) {
+      bar = paste0(bar, strrep(empty, bar_width - n_filled))
+    }
+    cat(sprintf(
+      "\r%s: %s%s%s %d/%d (%.1f%%)", prefix, lhs, bar, rhs,
+      current, total, 100 * frac
+    ))
+    utils::flush.console()
+  }
+
+  list(update = draw, close = function() cat("\n"))
+}
+
+
 # ------------------------------------------------------------------
 # sweep_prior_edge_density (internal)
 # ------------------------------------------------------------------
@@ -166,7 +252,8 @@ sweep_prior_edge_density = function(p, theta, delta,
                                     n_samples = 2000L, n_warmup = 500L,
                                     n_seeds = 3L,
                                     update_method = "gibbs",
-                                    cores = 1L, base_seed = 1L) {
+                                    cores = 1L, base_seed = 1L,
+                                    show_progress = FALSE) {
   cores = normalize_builder_cores(cores)
   num_pairs = p * (p - 1) / 2
   cells = expand.grid(theta = theta, seed = seq_len(n_seeds))
@@ -185,13 +272,42 @@ sweep_prior_edge_density = function(p, theta, delta,
     mean(draws$edge_indicators)
   }
 
+  # The bar follows display_progress (like the sampler's own bar), not the
+  # advisory bgms.verbose flag. It redraws in place with a carriage return, so
+  # it is drawn only in an interactive session; batch runs (scripts, R CMD
+  # check, the test suite) rely on the one-line build announcement in
+  # ggm_correction_table(). A parallel sweep forks the cells in batches of
+  # `cores` and advances the bar as each batch of chains completes.
+  n_cells = nrow(cells)
+  draw_bar = show_progress && interactive()
+  pb = NULL
+  if(draw_bar) {
+    pb = new_correction_progress(n_cells)
+    on.exit(pb$close(), add = TRUE)
+    pb$update(0L)
+  }
   edens_cells = if(cores > 1L) {
-    unlist(parallel::mclapply(
-      seq_len(nrow(cells)), one_cell,
-      mc.cores = cores, mc.preschedule = FALSE
-    ))
+    out = numeric(n_cells)
+    done = 0L
+    for(batch in split(seq_len(n_cells), ceiling(seq_len(n_cells) / cores))) {
+      out[batch] = as.numeric(unlist(parallel::mclapply(
+        batch, one_cell,
+        mc.cores = cores, mc.preschedule = FALSE
+      )))
+      done = done + length(batch)
+      if(draw_bar) {
+        pb$update(done)
+      }
+    }
+    out
   } else {
-    vapply(seq_len(nrow(cells)), one_cell, numeric(1))
+    vapply(seq_len(n_cells), function(k) {
+      v = one_cell(k)
+      if(draw_bar) {
+        pb$update(k)
+      }
+      v
+    }, numeric(1))
   }
   if(anyNA(edens_cells)) {
     stop("Correction-table sweep: a prior chain failed.")
@@ -242,7 +358,7 @@ build_ggm_correction_table = function(
   precision_scale_prior = gamma_prior(shape = 1, eta = 1),
   n_grid = 120L, n_samples = 2000L, n_warmup = 500L, n_seeds = 3L,
   update_method = c("gibbs", "adaptive-metropolis"),
-  cores = 1L, base_seed = 1L
+  cores = 1L, base_seed = 1L, show_progress = FALSE
 ) {
   update_method = match.arg(update_method)
   if(is.null(delta)) {
@@ -255,7 +371,8 @@ build_ggm_correction_table = function(
     interaction_prior = interaction_prior,
     precision_scale_prior = precision_scale_prior,
     n_samples = n_samples, n_warmup = n_warmup, n_seeds = n_seeds,
-    update_method = update_method, cores = cores, base_seed = base_seed
+    update_method = update_method, cores = cores, base_seed = base_seed,
+    show_progress = show_progress
   )
 
   table = correction_table_from_edens(
@@ -371,18 +488,20 @@ ggm_edge_prior_correction = function(prior, sampler, num_variables,
   interaction_prior = correction_interaction_prior(prior)
   precision_scale_prior = correction_scale_prior(prior)
 
-  if(isTRUE(sampler$verbose)) {
-    message(
-      "Edge-prior correction: building or loading the normalizing-constant ",
-      "table for this model (cached across fits)."
-    )
-  }
+  # The build announcement follows the advisory bgms.verbose flag; the progress
+  # bar follows display_progress (progress_type 0 is "none"), matching the
+  # sampler's own bar.
+  show_progress = is.null(sampler$progress_type) ||
+    !identical(as.integer(sampler$progress_type), 0L)
+
   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
+    cores = sampler$cores,
+    verbose = isTRUE(sampler$verbose),
+    show_progress = show_progress
   )
   if(identical(prior$edge_prior, "Stochastic-Block") &&
     is.null(table$fprime)) {
@@ -420,7 +539,8 @@ ggm_correction_table = function(
   precision_scale_prior = gamma_prior(shape = 1, eta = 1),
   n_grid = 120L, n_samples = 2000L, n_warmup = 500L, n_seeds = 3L,
   update_method = c("gibbs", "adaptive-metropolis"),
-  cores = 1L, base_seed = 1L, refresh = FALSE
+  cores = 1L, base_seed = 1L, refresh = FALSE, verbose = FALSE,
+  show_progress = FALSE
 ) {
   update_method = match.arg(update_method)
   if(is.null(delta)) {
@@ -451,13 +571,19 @@ ggm_correction_table = function(
     }
   }
 
+  if(verbose) {
+    message(
+      "Building the edge-selection prior correction table (one-time for ",
+      "this model size and prior; cached for later fits)."
+    )
+  }
   table = build_ggm_correction_table(
     p = p, delta = delta,
     interaction_prior = interaction_prior,
     precision_scale_prior = precision_scale_prior,
     n_grid = n_grid, n_samples = n_samples, n_warmup = n_warmup,
     n_seeds = n_seeds, update_method = update_method,
-    cores = cores, base_seed = base_seed
+    cores = cores, base_seed = base_seed, show_progress = show_progress
   )
 
   if(use_cache) {

R/sample_ggm_prior.R

@@ -354,7 +354,8 @@ sample_ggm_prior = function(
       p = p, delta = delta,
       interaction_prior = interaction_prior,
       precision_scale_prior = precision_scale_prior,
-      update_method = "gibbs"
+      update_method = "gibbs",
+      verbose = isTRUE(verbose)
     )
     correction = correction_list_from_table(table, ep$edge_prior)
   }