From e80a2cf20b0bb17558073cfd416276674bee1a23 Mon Sep 17 00:00:00 2001
From: Maarten Marsman <m.marsman@uva.nl>
Date: Mon, 6 Jul 2026 10:02:22 +0200
Subject: [PATCH 1/2] Add graph-level prior samplers and quiet the
 correction-table build

sample_graph_prior() draws edge indicators and edge-prior
hyperparameters from the graph level of the spike-and-slab prior. The
hierarchical specification samples ancestrally (hyperparameters from
their prior, then independent pair flips), which is exact. The joint
specification runs the zero-data prior chain and discards the precision
draws, so the returned graphs carry the per-graph normalizer tilt and
the hyperparameter updates apply the normalizing-constant correction.
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() draws block allocations and pair-inclusion
probabilities ancestrally from the MFM-SBM hyperprior, wrapping the
internal helpers that until now only initialized the joint chain.

The correction-table build message printed before the cache check, so
every fit read 'building or loading' even when the table loaded
instantly. The message now prints only on an actual build, states that
the build is one-time and cached, and the sweep reports progress: a
progress bar on one core, a cell-and-core count when parallel.

Stage 1 of the prior-sampler redesign
(dev/plans/active/2026-07-06_prior-sampler-redesign_DESIGN.md).
---
 DESCRIPTION                              |  13 +-
 NAMESPACE                                |   2 +
 NEWS.md                                  |   3 +
 R/correction_tables.R                    | 131 ++++++-
 R/sample_ggm_prior.R                     |   3 +-
 R/sample_graph_prior.R                   | 427 +++++++++++++++++++++++
 man/sample_graph_prior.Rd                | 148 ++++++++
 man/sample_sbm_prior.Rd                  |  52 +++
 tests/testthat/test-correction-tables.R  |  32 ++
 tests/testthat/test-mixed-correction.R   |  27 +-
 tests/testthat/test-sample-graph-prior.R | 152 ++++++++
 tests/testthat/test-sample-sbm-prior.R   |  59 ++++
 12 files changed, 1021 insertions(+), 28 deletions(-)
 create mode 100644 R/sample_graph_prior.R
 create mode 100644 man/sample_graph_prior.Rd
 create mode 100644 man/sample_sbm_prior.Rd
 create mode 100644 tests/testthat/test-sample-graph-prior.R
 create mode 100644 tests/testthat/test-sample-sbm-prior.R

diff --git a/DESCRIPTION b/DESCRIPTION
index a5c0ddd0..701d9a2f 100644
--- a/DESCRIPTION
+++ b/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, 
diff --git a/NAMESPACE b/NAMESPACE
index bfb05b15..c7e82c2f 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/NEWS.md b/NEWS.md
index 0f668e1c..467de928 100644
--- a/NEWS.md
+++ b/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), states that the build is one-time and cached, and reports progress (a progress bar on one core, a cell-and-core count otherwise).
 * 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.
diff --git a/R/correction_tables.R b/R/correction_tables.R
index 555d9326..3eb8fd55 100644
--- a/R/correction_tables.R
+++ b/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,
+                                    verbose = FALSE) {
   cores = normalize_builder_cores(cores)
   num_pairs = p * (p - 1) / 2
   cells = expand.grid(theta = theta, seed = seq_len(n_seeds))
@@ -185,11 +272,29 @@ sweep_prior_edge_density = function(p, theta, delta,
     mean(draws$edge_indicators)
   }
 
+  # A progress bar is drawn only in an interactive session; in batch runs
+  # (scripts, R CMD check, the test suite) the one-line build announcement
+  # in ggm_correction_table() is the only output.
+  show_bar = verbose && interactive() && cores == 1L
   edens_cells = if(cores > 1L) {
+    if(verbose) {
+      message(sprintf(
+        "Sweeping %d prior cells on %d cores.", nrow(cells), cores
+      ))
+    }
     unlist(parallel::mclapply(
       seq_len(nrow(cells)), one_cell,
       mc.cores = cores, mc.preschedule = FALSE
     ))
+  } else if(show_bar) {
+    pb = new_correction_progress(nrow(cells))
+    on.exit(pb$close(), add = TRUE)
+    pb$update(0L)
+    vapply(seq_len(nrow(cells)), function(k) {
+      v = one_cell(k)
+      pb$update(k)
+      v
+    }, numeric(1))
   } else {
     vapply(seq_len(nrow(cells)), one_cell, numeric(1))
   }
@@ -242,7 +347,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, verbose = FALSE
 ) {
   update_method = match.arg(update_method)
   if(is.null(delta)) {
@@ -255,7 +360,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,
+    verbose = verbose
   )
 
   table = correction_table_from_edens(
@@ -371,18 +477,13 @@ 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)."
-    )
-  }
   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)
   )
   if(identical(prior$edge_prior, "Stochastic-Block") &&
     is.null(table$fprime)) {
@@ -420,7 +521,7 @@ 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
 ) {
   update_method = match.arg(update_method)
   if(is.null(delta)) {
@@ -451,13 +552,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, verbose = verbose
   )
 
   if(use_cache) {
diff --git a/R/sample_ggm_prior.R b/R/sample_ggm_prior.R
index 18cffb4e..7311abd6 100644
--- a/R/sample_ggm_prior.R
+++ b/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)
   }
diff --git a/R/sample_graph_prior.R b/R/sample_graph_prior.R
new file mode 100644
index 00000000..b15dae19
--- /dev/null
+++ b/R/sample_graph_prior.R
@@ -0,0 +1,427 @@
+# ==============================================================================
+# Graph-level prior samplers
+# ==============================================================================
+#
+# sample_graph_prior(): draws (hyperparameters, edge indicators) from the
+# graph level of the spike-and-slab prior, under either specification of
+# how the precision prior composes with the graph:
+#
+#   hierarchical: p(hyper) . pi(Gamma | hyper) -- ancestral and exact
+#   joint:        p(hyper) . q(Gamma | hyper), with
+#                 q(Gamma | hyper) proportional to Z(Gamma) pi(Gamma | hyper)
+#                 -- the determinant-tilted law, sampled by the zero-data
+#                 (K, Gamma) chain with K discarded from the output
+#
+# sample_sbm_prior(): ancestral draws of (allocations, pair probabilities)
+# from the MFM-SBM edge-prior hyperprior.
+#
+# Edge indicators are returned in row-major upper-triangle order (i < j),
+# matching the K_offdiag column order of sample_ggm_prior().
+# ==============================================================================
+
+
+# Run fn() inside an RNG scope keyed to `seed`, restoring the caller's RNG
+# state on exit (same idiom as ggm_prior_ancestral_indicators).
+with_graph_prior_seed = function(seed, fn) {
+  has_seed = exists(".Random.seed", envir = globalenv(), inherits = FALSE)
+  if(has_seed) {
+    old_seed = get(".Random.seed", envir = globalenv(), inherits = FALSE)
+    on.exit(assign(".Random.seed", old_seed, envir = globalenv()), add = TRUE)
+  } else {
+    on.exit(
+      if(exists(".Random.seed", envir = globalenv(), inherits = FALSE)) {
+        rm(list = ".Random.seed", envir = globalenv())
+      },
+      add = TRUE
+    )
+  }
+  set.seed(seed)
+  fn()
+}
+
+
+# Row-major upper-triangle pair indices (i < j): a (E x 2) matrix with
+# E = p(p-1)/2 rows, ordered (1,2), (1,3), ..., (1,p), (2,3), ...
+graph_pair_indices = function(p) {
+  ii = rep(seq_len(p - 1L), times = (p - 1L):1L)
+  jj = unlist(lapply(seq_len(p - 1L), function(i) (i + 1L):p))
+  cbind(ii, jj)
+}
+
+
+# Resolve the optional conditioning arguments to a fixed p x p pair
+# probability matrix, or NULL when no conditioning was requested.
+graph_prior_conditioning = function(ep, p, theta, allocations, block_probs) {
+  if(!is.null(theta) && (!is.null(allocations) || !is.null(block_probs))) {
+    stop(
+      "Supply either 'theta' or ('allocations', 'block_probs') to condition ",
+      "the graph prior, not both."
+    )
+  }
+  if(!is.null(theta)) {
+    if(identical(ep$edge_prior, "Stochastic-Block")) {
+      stop(
+        "Condition the Stochastic-Block prior with 'allocations' and ",
+        "'block_probs'; 'theta' conditions the Bernoulli and Beta-Bernoulli ",
+        "priors."
+      )
+    }
+    if(!is.numeric(theta) || length(theta) != 1L || is.na(theta) ||
+      theta <= 0 || theta >= 1) {
+      stop("'theta' must be a single numeric in (0, 1).")
+    }
+    return(matrix(theta, p, p))
+  }
+  if(!is.null(allocations) || !is.null(block_probs)) {
+    if(is.null(allocations) || is.null(block_probs)) {
+      stop(
+        "Conditioning the Stochastic-Block prior needs both 'allocations' ",
+        "and 'block_probs'."
+      )
+    }
+    if(!identical(ep$edge_prior, "Stochastic-Block")) {
+      stop(
+        "'allocations'/'block_probs' condition the Stochastic-Block prior; ",
+        "use edge_prior = sbm_prior()."
+      )
+    }
+    allocations = as.integer(allocations)
+    if(length(allocations) != p || anyNA(allocations) ||
+      any(allocations < 1L)) {
+      stop("'allocations' must be a length-p vector of 1-based block labels.")
+    }
+    if(!is.matrix(block_probs) || nrow(block_probs) != ncol(block_probs) ||
+      !isSymmetric(unname(block_probs)) ||
+      any(block_probs <= 0) || any(block_probs >= 1)) {
+      stop(
+        "'block_probs' must be a symmetric matrix with entries in (0, 1)."
+      )
+    }
+    if(max(allocations) > nrow(block_probs)) {
+      stop("'allocations' labels exceed the size of 'block_probs'.")
+    }
+    prob = matrix(0.5, p, p)
+    for(i in seq_len(p - 1L)) {
+      for(j in (i + 1L):p) {
+        prob[i, j] = block_probs[allocations[i], allocations[j]]
+        prob[j, i] = prob[i, j]
+      }
+    }
+    return(prob)
+  }
+  NULL
+}
+
+
+#' @title Sample from the Graph Prior
+#'
+#' @description
+#' Draws edge-inclusion indicators, together with any edge-prior
+#' hyperparameters, from the graph level of the spike-and-slab prior used by
+#' \code{\link{bgm}} for models with continuous variables. The \code{spec}
+#' argument selects how the precision prior composes with the graph:
+#' \itemize{
+#'   \item \code{"hierarchical"} (default): the graph marginal is exactly the
+#'     edge prior, \eqn{p(\mathrm{hyper}) \, \pi(\Gamma \mid \mathrm{hyper})}.
+#'     Sampling is ancestral and exact: hyperparameters from their prior,
+#'     then independent pair flips.
+#'   \item \code{"joint"}: the graph marginal is reweighted by the per-graph
+#'     normalizer of the determinant-tilted precision prior,
+#'     \eqn{q(\Gamma \mid \mathrm{hyper}) \propto Z(\Gamma) \,
+#'     \pi(\Gamma \mid \mathrm{hyper})}. Sampling runs the zero-data
+#'     \eqn{(K, \Gamma)} chain of \code{\link{sample_ggm_prior}} and discards
+#'     \eqn{K}; with a Beta-Bernoulli or Stochastic-Block prior the
+#'     hyperparameter updates apply the normalizing-constant correction, so
+#'     the first call for a model cell may build the correction table
+#'     (cached across fits).
+#' }
+#'
+#' @details
+#' The optional conditioning arguments fix the edge-prior hyperparameters
+#' instead of sampling them: \code{theta} fixes the inclusion probability of
+#' a Bernoulli or Beta-Bernoulli prior, and \code{allocations} plus
+#' \code{block_probs} fix the block structure of a Stochastic-Block prior
+#' (reducing it to independent pair flips at the given block probabilities).
+#'
+#' Under \code{spec = "joint"} the tilted graph law depends on the precision
+#' prior through its normalizer, so \code{interaction_prior},
+#' \code{precision_scale_prior}, and \code{delta} are part of the graph law;
+#' they are ignored under \code{spec = "hierarchical"}.
+#'
+#' @param p Integer. Number of nodes (\eqn{p \ge 2}).
+#' @param n_samples Integer. Number of prior draws.
+#' @param edge_prior An edge prior specification object:
+#'   \code{\link{bernoulli_prior}()}, \code{\link{beta_bernoulli_prior}()},
+#'   or \code{\link{sbm_prior}()}. Default \code{bernoulli_prior(0.5)}.
+#' @param spec One of \code{"hierarchical"} (default) or \code{"joint"}.
+#' @param interaction_prior A \code{\link{cauchy_prior}()} or
+#'   \code{\link{normal_prior}()} for the pairwise (slab) part of the
+#'   precision prior. Used only when \code{spec = "joint"}.
+#' @param precision_scale_prior A \code{\link{gamma_prior}()} or
+#'   \code{\link{exponential_prior}()} for the precision diagonal. Used only
+#'   when \code{spec = "joint"}.
+#' @param delta Non-negative numeric or \code{NULL} (default): determinant
+#'   tilt exponent; \code{NULL} resolves to \eqn{0.5 \log(p)}. Used only
+#'   when \code{spec = "joint"}.
+#' @param theta Optional numeric in (0, 1): fix the inclusion probability of
+#'   a Bernoulli or Beta-Bernoulli edge prior instead of sampling it.
+#' @param allocations Optional integer vector of length \code{p} with 1-based
+#'   block labels: fix the Stochastic-Block allocation instead of sampling it.
+#'   Requires \code{block_probs}.
+#' @param block_probs Optional symmetric matrix with entries in (0, 1): the
+#'   block-pair inclusion probabilities that go with \code{allocations}.
+#' @param n_warmup Integer. Warmup iterations of the zero-data chain. Used
+#'   only when \code{spec = "joint"}. Default \code{2e3}.
+#' @param seed Integer. Seed for the draw; the caller's RNG state is
+#'   restored on exit.
+#' @param verbose Logical. Print progress of the zero-data chain and of a
+#'   correction-table build. Default \code{TRUE}.
+#'
+#' @return A list with elements:
+#'   \describe{
+#'     \item{\code{edge_indicators}}{Integer matrix
+#'       (\code{n_samples x p(p-1)/2}) of edge-inclusion indicators, columns
+#'       in row-major upper-triangle order (matching
+#'       \code{sample_ggm_prior()}'s \code{K_offdiag}).}
+#'     \item{\code{pair_names}}{Character vector labeling the columns as
+#'       \code{"i-j"}.}
+#'     \item{\code{theta}}{Only with an unconditioned
+#'       \code{beta_bernoulli_prior()}: numeric vector of sampled inclusion
+#'       probabilities.}
+#'     \item{\code{allocations}}{Only with an unconditioned
+#'       \code{sbm_prior()}: integer matrix (\code{n_samples x p}) of sampled
+#'       block allocations.}
+#'     \item{\code{spec}, \code{edge_prior}, \code{p}}{The specification,
+#'       edge-prior family, and node count of the draw.}
+#'   }
+#'
+#' @examples
+#' # Hierarchical spec: the graph marginal is exactly the edge prior.
+#' g = sample_graph_prior(
+#'   p = 6, n_samples = 200,
+#'   edge_prior = bernoulli_prior(0.3), seed = 11
+#' )
+#' mean(g$edge_indicators) # about 0.3
+#'
+#' # Beta-Bernoulli: inclusion probabilities are sampled alongside.
+#' g = sample_graph_prior(
+#'   p = 6, n_samples = 200,
+#'   edge_prior = beta_bernoulli_prior(2, 4), seed = 11
+#' )
+#' mean(g$theta) # about 1/3
+#'
+#' \donttest{
+#' # Joint spec: the graph law carries the per-graph normalizer Z(Gamma).
+#' g = sample_graph_prior(
+#'   p = 6, n_samples = 500,
+#'   edge_prior = bernoulli_prior(0.3), spec = "joint",
+#'   interaction_prior = normal_prior(scale = 0.5),
+#'   precision_scale_prior = gamma_prior(shape = 1, rate = 2),
+#'   seed = 11, verbose = FALSE
+#' )
+#' mean(g$edge_indicators) # shifted away from 0.3 by the Z(Gamma) tilt
+#' }
+#' @seealso \code{\link{sample_ggm_prior}}, \code{\link{sample_sbm_prior}},
+#'   \code{\link{bernoulli_prior}}, \code{\link{beta_bernoulli_prior}},
+#'   \code{\link{sbm_prior}}, \code{\link{bgm}}
+#'
+#' @export
+sample_graph_prior = function(
+  p,
+  n_samples,
+  edge_prior = bernoulli_prior(0.5),
+  spec = c("hierarchical", "joint"),
+  interaction_prior = cauchy_prior(scale = 2.5),
+  precision_scale_prior = gamma_prior(shape = 1, eta = 1),
+  delta = NULL,
+  theta = NULL,
+  allocations = NULL,
+  block_probs = NULL,
+  n_warmup = 2e3,
+  seed = 1L,
+  verbose = TRUE
+) {
+  spec = match.arg(spec)
+  if(!is.numeric(p) || length(p) != 1L || is.na(p) || p < 2 || p != round(p)) {
+    stop("'p' must be a single integer >= 2.")
+  }
+  p = as.integer(p)
+  if(!is.numeric(n_samples) || length(n_samples) != 1L || is.na(n_samples) ||
+    n_samples < 1 || n_samples != round(n_samples)) {
+    stop("'n_samples' must be a single positive integer.")
+  }
+  n_samples = as.integer(n_samples)
+
+  ep = unpack_indicator_prior(edge_prior, num_variables = p)
+  cond_prob = graph_prior_conditioning(ep, p, theta, allocations, block_probs)
+
+  pairs = graph_pair_indices(p)
+  pair_names = paste0(pairs[, 1L], "-", pairs[, 2L])
+
+  if(spec == "joint") {
+    eff_edge_prior = if(is.null(cond_prob)) {
+      edge_prior
+    } else {
+      bernoulli_prior(inclusion_probability = cond_prob)
+    }
+    draws = sample_ggm_prior(
+      p = p, n_samples = n_samples, n_warmup = as.integer(n_warmup),
+      interaction_prior = interaction_prior,
+      precision_scale_prior = precision_scale_prior,
+      delta = delta,
+      spec = "joint",
+      edge_prior = eff_edge_prior,
+      update_method = "gibbs",
+      seed = as.integer(seed),
+      verbose = verbose
+    )
+    out = list(
+      edge_indicators = draws$edge_indicators,
+      pair_names = pair_names,
+      theta = draws$theta,
+      allocations = draws$allocations,
+      spec = spec,
+      edge_prior = ep$edge_prior,
+      p = p
+    )
+    return(out)
+  }
+
+  # Hierarchical spec: ancestral. Hyperparameters from their prior, then
+  # independent pair flips given the implied pair probabilities.
+  res = with_graph_prior_seed(seed, function() {
+    n_edges = nrow(pairs)
+    theta_out = NULL
+    alloc_out = NULL
+
+    if(!is.null(cond_prob) || identical(ep$edge_prior, "Bernoulli")) {
+      prob = if(is.null(cond_prob)) ep$inclusion_probability else cond_prob
+      pv = prob[pairs]
+      ru = matrix(runif(n_samples * n_edges), n_samples, n_edges)
+      indicators = 1L * (ru < matrix(pv, n_samples, n_edges, byrow = TRUE))
+    } else if(identical(ep$edge_prior, "Beta-Bernoulli")) {
+      theta_out = rbeta(
+        n_samples, ep$beta_bernoulli_alpha, ep$beta_bernoulli_beta
+      )
+      ru = matrix(runif(n_samples * n_edges), n_samples, n_edges)
+      indicators = 1L * (ru < theta_out)
+    } else {
+      alloc_out = matrix(0L, n_samples, p)
+      indicators = matrix(0L, n_samples, n_edges)
+      for(s in seq_len(n_samples)) {
+        z = ancestral_mfm_sbm_partition(p, ep$lambda, ep$dirichlet_alpha)
+        prob = ancestral_sbm_pair_probabilities(
+          z,
+          ep$beta_bernoulli_alpha, ep$beta_bernoulli_beta,
+          ep$beta_bernoulli_alpha_between, ep$beta_bernoulli_beta_between
+        )
+        alloc_out[s, ] = as.integer(z)
+        indicators[s, ] = as.integer(runif(n_edges) < prob[pairs])
+      }
+    }
+    storage.mode(indicators) = "integer"
+    list(indicators = indicators, theta = theta_out, allocations = alloc_out)
+  })
+
+  list(
+    edge_indicators = res$indicators,
+    pair_names = pair_names,
+    theta = res$theta,
+    allocations = res$allocations,
+    spec = spec,
+    edge_prior = ep$edge_prior,
+    p = p
+  )
+}
+
+
+#' @title Sample from the Stochastic-Block Edge-Prior Hyperprior
+#'
+#' @description
+#' Ancestral draws from the MFM-SBM hyperprior of
+#' \code{\link{sbm_prior}()}: a partition of the nodes into blocks
+#' (shifted-Poisson number of components, Dirichlet-weighted allocation) and
+#' Beta-distributed within- and between-block edge-inclusion probabilities.
+#' These are the hyperparameters that \code{\link{sample_graph_prior}} and
+#' \code{\link{bgm}} integrate over when the edge prior is a Stochastic-Block
+#' prior.
+#'
+#' @param p Integer. Number of nodes (\eqn{p \ge 2}).
+#' @param n_samples Integer. Number of prior draws.
+#' @param edge_prior An \code{\link{sbm_prior}()} object. Default
+#'   \code{sbm_prior()}.
+#' @param seed Integer. Seed for the draw; the caller's RNG state is
+#'   restored on exit.
+#'
+#' @return A list with elements:
+#'   \describe{
+#'     \item{\code{allocations}}{Integer matrix (\code{n_samples x p}) of
+#'       block labels.}
+#'     \item{\code{pair_probability}}{Numeric matrix
+#'       (\code{n_samples x p(p-1)/2}) of implied pair-inclusion
+#'       probabilities, columns in row-major upper-triangle order.}
+#'     \item{\code{num_blocks}}{Integer vector: number of occupied blocks
+#'       per draw.}
+#'     \item{\code{pair_names}}{Character vector labeling the pair columns
+#'       as \code{"i-j"}.}
+#'     \item{\code{p}}{The node count.}
+#'   }
+#'
+#' @examples
+#' draws = sample_sbm_prior(p = 8, n_samples = 100, seed = 4)
+#' table(draws$num_blocks)
+#' range(draws$pair_probability)
+#' @seealso \code{\link{sbm_prior}}, \code{\link{sample_graph_prior}},
+#'   \code{\link{bgm}}
+#'
+#' @export
+sample_sbm_prior = function(p, n_samples, edge_prior = sbm_prior(),
+                            seed = 1L) {
+  if(!is.numeric(p) || length(p) != 1L || is.na(p) || p < 2 || p != round(p)) {
+    stop("'p' must be a single integer >= 2.")
+  }
+  p = as.integer(p)
+  if(!is.numeric(n_samples) || length(n_samples) != 1L || is.na(n_samples) ||
+    n_samples < 1 || n_samples != round(n_samples)) {
+    stop("'n_samples' must be a single positive integer.")
+  }
+  n_samples = as.integer(n_samples)
+
+  ep = unpack_indicator_prior(edge_prior, num_variables = p)
+  if(!identical(ep$edge_prior, "Stochastic-Block")) {
+    stop("'edge_prior' must be an sbm_prior() object.")
+  }
+
+  pairs = graph_pair_indices(p)
+  pair_names = paste0(pairs[, 1L], "-", pairs[, 2L])
+
+  res = with_graph_prior_seed(seed, function() {
+    allocations = matrix(0L, n_samples, p)
+    pair_probability = matrix(0, n_samples, nrow(pairs))
+    num_blocks = integer(n_samples)
+    for(s in seq_len(n_samples)) {
+      z = ancestral_mfm_sbm_partition(p, ep$lambda, ep$dirichlet_alpha)
+      prob = ancestral_sbm_pair_probabilities(
+        z,
+        ep$beta_bernoulli_alpha, ep$beta_bernoulli_beta,
+        ep$beta_bernoulli_alpha_between, ep$beta_bernoulli_beta_between
+      )
+      allocations[s, ] = as.integer(z)
+      pair_probability[s, ] = prob[pairs]
+      num_blocks[s] = length(unique(z))
+    }
+    list(
+      allocations = allocations,
+      pair_probability = pair_probability,
+      num_blocks = num_blocks
+    )
+  })
+
+  list(
+    allocations = res$allocations,
+    pair_probability = res$pair_probability,
+    num_blocks = res$num_blocks,
+    pair_names = pair_names,
+    p = p
+  )
+}
diff --git a/man/sample_graph_prior.Rd b/man/sample_graph_prior.Rd
new file mode 100644
index 00000000..fce4a00c
--- /dev/null
+++ b/man/sample_graph_prior.Rd
@@ -0,0 +1,148 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sample_graph_prior.R
+\name{sample_graph_prior}
+\alias{sample_graph_prior}
+\title{Sample from the Graph Prior}
+\usage{
+sample_graph_prior(
+  p,
+  n_samples,
+  edge_prior = bernoulli_prior(0.5),
+  spec = c("hierarchical", "joint"),
+  interaction_prior = cauchy_prior(scale = 2.5),
+  precision_scale_prior = gamma_prior(shape = 1, eta = 1),
+  delta = NULL,
+  theta = NULL,
+  allocations = NULL,
+  block_probs = NULL,
+  n_warmup = 2000,
+  seed = 1L,
+  verbose = TRUE
+)
+}
+\arguments{
+\item{p}{Integer. Number of nodes (\eqn{p \ge 2}).}
+
+\item{n_samples}{Integer. Number of prior draws.}
+
+\item{edge_prior}{An edge prior specification object:
+\code{\link{bernoulli_prior}()}, \code{\link{beta_bernoulli_prior}()},
+or \code{\link{sbm_prior}()}. Default \code{bernoulli_prior(0.5)}.}
+
+\item{spec}{One of \code{"hierarchical"} (default) or \code{"joint"}.}
+
+\item{interaction_prior}{A \code{\link{cauchy_prior}()} or
+\code{\link{normal_prior}()} for the pairwise (slab) part of the
+precision prior. Used only when \code{spec = "joint"}.}
+
+\item{precision_scale_prior}{A \code{\link{gamma_prior}()} or
+\code{\link{exponential_prior}()} for the precision diagonal. Used only
+when \code{spec = "joint"}.}
+
+\item{delta}{Non-negative numeric or \code{NULL} (default): determinant
+tilt exponent; \code{NULL} resolves to \eqn{0.5 \log(p)}. Used only
+when \code{spec = "joint"}.}
+
+\item{theta}{Optional numeric in (0, 1): fix the inclusion probability of
+a Bernoulli or Beta-Bernoulli edge prior instead of sampling it.}
+
+\item{allocations}{Optional integer vector of length \code{p} with 1-based
+block labels: fix the Stochastic-Block allocation instead of sampling it.
+Requires \code{block_probs}.}
+
+\item{block_probs}{Optional symmetric matrix with entries in (0, 1): the
+block-pair inclusion probabilities that go with \code{allocations}.}
+
+\item{n_warmup}{Integer. Warmup iterations of the zero-data chain. Used
+only when \code{spec = "joint"}. Default \code{2e3}.}
+
+\item{seed}{Integer. Seed for the draw; the caller's RNG state is
+restored on exit.}
+
+\item{verbose}{Logical. Print progress of the zero-data chain and of a
+correction-table build. Default \code{TRUE}.}
+}
+\value{
+A list with elements:
+\describe{
+\item{\code{edge_indicators}}{Integer matrix
+(\code{n_samples x p(p-1)/2}) of edge-inclusion indicators, columns
+in row-major upper-triangle order (matching
+\code{sample_ggm_prior()}'s \code{K_offdiag}).}
+\item{\code{pair_names}}{Character vector labeling the columns as
+\code{"i-j"}.}
+\item{\code{theta}}{Only with an unconditioned
+\code{beta_bernoulli_prior()}: numeric vector of sampled inclusion
+probabilities.}
+\item{\code{allocations}}{Only with an unconditioned
+\code{sbm_prior()}: integer matrix (\code{n_samples x p}) of sampled
+block allocations.}
+\item{\code{spec}, \code{edge_prior}, \code{p}}{The specification,
+edge-prior family, and node count of the draw.}
+}
+}
+\description{
+Draws edge-inclusion indicators, together with any edge-prior
+hyperparameters, from the graph level of the spike-and-slab prior used by
+\code{\link{bgm}} for models with continuous variables. The \code{spec}
+argument selects how the precision prior composes with the graph:
+\itemize{
+\item \code{"hierarchical"} (default): the graph marginal is exactly the
+edge prior, \eqn{p(\mathrm{hyper}) \, \pi(\Gamma \mid \mathrm{hyper})}.
+Sampling is ancestral and exact: hyperparameters from their prior,
+then independent pair flips.
+\item \code{"joint"}: the graph marginal is reweighted by the per-graph
+normalizer of the determinant-tilted precision prior,
+\eqn{q(\Gamma \mid \mathrm{hyper}) \propto Z(\Gamma) \,
+    \pi(\Gamma \mid \mathrm{hyper})}. Sampling runs the zero-data
+\eqn{(K, \Gamma)} chain of \code{\link{sample_ggm_prior}} and discards
+\eqn{K}; with a Beta-Bernoulli or Stochastic-Block prior the
+hyperparameter updates apply the normalizing-constant correction, so
+the first call for a model cell may build the correction table
+(cached across fits).
+}
+}
+\details{
+The optional conditioning arguments fix the edge-prior hyperparameters
+instead of sampling them: \code{theta} fixes the inclusion probability of
+a Bernoulli or Beta-Bernoulli prior, and \code{allocations} plus
+\code{block_probs} fix the block structure of a Stochastic-Block prior
+(reducing it to independent pair flips at the given block probabilities).
+
+Under \code{spec = "joint"} the tilted graph law depends on the precision
+prior through its normalizer, so \code{interaction_prior},
+\code{precision_scale_prior}, and \code{delta} are part of the graph law;
+they are ignored under \code{spec = "hierarchical"}.
+}
+\examples{
+# Hierarchical spec: the graph marginal is exactly the edge prior.
+g = sample_graph_prior(
+  p = 6, n_samples = 200,
+  edge_prior = bernoulli_prior(0.3), seed = 11
+)
+mean(g$edge_indicators) # about 0.3
+
+# Beta-Bernoulli: inclusion probabilities are sampled alongside.
+g = sample_graph_prior(
+  p = 6, n_samples = 200,
+  edge_prior = beta_bernoulli_prior(2, 4), seed = 11
+)
+mean(g$theta) # about 1/3
+
+\donttest{
+# Joint spec: the graph law carries the per-graph normalizer Z(Gamma).
+g = sample_graph_prior(
+  p = 6, n_samples = 500,
+  edge_prior = bernoulli_prior(0.3), spec = "joint",
+  interaction_prior = normal_prior(scale = 0.5),
+  precision_scale_prior = gamma_prior(shape = 1, rate = 2),
+  seed = 11, verbose = FALSE
+)
+mean(g$edge_indicators) # shifted away from 0.3 by the Z(Gamma) tilt
+}
+}
+\seealso{
+\code{\link{sample_ggm_prior}}, \code{\link{sample_sbm_prior}},
+\code{\link{bernoulli_prior}}, \code{\link{beta_bernoulli_prior}},
+\code{\link{sbm_prior}}, \code{\link{bgm}}
+}
diff --git a/man/sample_sbm_prior.Rd b/man/sample_sbm_prior.Rd
new file mode 100644
index 00000000..85a2cce8
--- /dev/null
+++ b/man/sample_sbm_prior.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sample_graph_prior.R
+\name{sample_sbm_prior}
+\alias{sample_sbm_prior}
+\title{Sample from the Stochastic-Block Edge-Prior Hyperprior}
+\usage{
+sample_sbm_prior(p, n_samples, edge_prior = sbm_prior(), seed = 1L)
+}
+\arguments{
+\item{p}{Integer. Number of nodes (\eqn{p \ge 2}).}
+
+\item{n_samples}{Integer. Number of prior draws.}
+
+\item{edge_prior}{An \code{\link{sbm_prior}()} object. Default
+\code{sbm_prior()}.}
+
+\item{seed}{Integer. Seed for the draw; the caller's RNG state is
+restored on exit.}
+}
+\value{
+A list with elements:
+\describe{
+\item{\code{allocations}}{Integer matrix (\code{n_samples x p}) of
+block labels.}
+\item{\code{pair_probability}}{Numeric matrix
+(\code{n_samples x p(p-1)/2}) of implied pair-inclusion
+probabilities, columns in row-major upper-triangle order.}
+\item{\code{num_blocks}}{Integer vector: number of occupied blocks
+per draw.}
+\item{\code{pair_names}}{Character vector labeling the pair columns
+as \code{"i-j"}.}
+\item{\code{p}}{The node count.}
+}
+}
+\description{
+Ancestral draws from the MFM-SBM hyperprior of
+\code{\link{sbm_prior}()}: a partition of the nodes into blocks
+(shifted-Poisson number of components, Dirichlet-weighted allocation) and
+Beta-distributed within- and between-block edge-inclusion probabilities.
+These are the hyperparameters that \code{\link{sample_graph_prior}} and
+\code{\link{bgm}} integrate over when the edge prior is a Stochastic-Block
+prior.
+}
+\examples{
+draws = sample_sbm_prior(p = 8, n_samples = 100, seed = 4)
+table(draws$num_blocks)
+range(draws$pair_probability)
+}
+\seealso{
+\code{\link{sbm_prior}}, \code{\link{sample_graph_prior}},
+\code{\link{bgm}}
+}
diff --git a/tests/testthat/test-correction-tables.R b/tests/testthat/test-correction-tables.R
index a9361da2..f6228e3b 100644
--- a/tests/testthat/test-correction-tables.R
+++ b/tests/testthat/test-correction-tables.R
@@ -108,6 +108,38 @@ test_that("table build with cache round-trips and reuses the file", {
   expect_equal(tab1$cell$delta, 0.5 * log(4))
 })
 
+test_that("the progress bar renders the label, counts, and percentage", {
+  pb = new_correction_progress(120L, prefix = "Correction table")
+  mid = paste(capture.output(pb$update(70L)), collapse = "")
+  expect_match(mid, "Correction table:", fixed = TRUE)
+  expect_match(mid, "70/120", fixed = TRUE)
+  expect_match(mid, "58.3%", fixed = TRUE)
+  full = paste(capture.output(pb$update(120L)), collapse = "")
+  expect_match(full, "120/120 (100.0%)", fixed = TRUE)
+})
+
+test_that("the build announces itself once; a cache hit is silent", {
+  cache_dir = file.path(tempdir(), "bgms-ctable-msg-test")
+  unlink(cache_dir, recursive = TRUE)
+  old = options(bgms.correction_cache_dir = cache_dir)
+  on.exit(options(old), add = TRUE)
+
+  build_args = list(
+    p = 4, n_grid = 12L, n_samples = 100L, n_warmup = 100L, n_seeds = 1L,
+    update_method = "gibbs", verbose = TRUE
+  )
+  # capture.output silences the progress bar (stdout); messages pass through.
+  capture.output(
+    expect_message(
+      do.call(ggm_correction_table, build_args),
+      "Building the edge-selection prior correction table"
+    )
+  )
+  capture.output(
+    expect_no_message(do.call(ggm_correction_table, build_args))
+  )
+})
+
 test_that("gibbs and adaptive-metropolis sweeps agree on edge density", {
   skip_on_cran()
 
diff --git a/tests/testthat/test-mixed-correction.R b/tests/testthat/test-mixed-correction.R
index 00d0248e..c426418c 100644
--- a/tests/testthat/test-mixed-correction.R
+++ b/tests/testthat/test-mixed-correction.R
@@ -242,19 +242,28 @@ test_that("bgm() applies the correction to a mixed beta-bernoulli fit", {
   skip_on_cran()
   old = mixed_correction_cache()
   on.exit(options(old), add = TRUE)
+  # A fresh cache forces an actual build, the only case that announces
+  # itself; a cache hit is silent.
+  unlink(getOption("bgms.correction_cache_dir"), recursive = TRUE)
 
   d = mixed_smoke_data()
-  msgs = capture.output(
-    fit <- bgm(d$x,
-      variable_type = d$variable_type,
-      edge_prior = beta_bernoulli_prior(),
-      iter = 300, warmup = 300, chains = 1, cores = 1,
-      display_progress = "none", verbose = TRUE
+  # Capture stdout too, so an interactive test run does not leak the
+  # build progress bar; the announcement is a message (stderr).
+  msgs = NULL
+  capture.output(
+    msgs <- capture.output(
+      fit <- bgm(d$x,
+        variable_type = d$variable_type,
+        edge_prior = beta_bernoulli_prior(),
+        iter = 300, warmup = 300, chains = 1, cores = 1,
+        display_progress = "none", verbose = TRUE
+      ),
+      type = "message"
     ),
-    type = "message"
+    type = "output"
   )
 
-  expect_true(any(grepl("Edge-prior correction", msgs)))
+  expect_true(any(grepl("correction table", msgs)))
   expect_s3_class(fit, "bgms")
   th = fit$inclusion_parameter_samples
   expect_length(th, 1L)
@@ -294,7 +303,7 @@ test_that("bgm() skips the correction with one continuous variable", {
     type = "message"
   )
 
-  expect_false(any(grepl("Edge-prior correction", msgs)))
+  expect_false(any(grepl("correction table", msgs)))
   expect_s3_class(fit, "bgms")
 })
 
diff --git a/tests/testthat/test-sample-graph-prior.R b/tests/testthat/test-sample-graph-prior.R
new file mode 100644
index 00000000..08d377bc
--- /dev/null
+++ b/tests/testthat/test-sample-graph-prior.R
@@ -0,0 +1,152 @@
+# Tests for sample_graph_prior(): the hierarchical spec's ancestral graph
+# law, hyperparameter conditioning, and the joint spec's Z(Gamma)-tilted
+# negative control (mirroring test-hier-zratio-identity).
+
+test_that("hierarchical Bernoulli graph law matches its prior", {
+  g = sample_graph_prior(
+    p = 6, n_samples = 4000,
+    edge_prior = bernoulli_prior(0.3), seed = 11
+  )
+  expect_identical(dim(g$edge_indicators), c(4000L, 15L))
+  expect_true(all(g$edge_indicators %in% 0:1))
+  expect_lt(abs(mean(g$edge_indicators) - 0.3), 0.02)
+  expect_null(g$theta)
+  expect_null(g$allocations)
+  expect_identical(g$spec, "hierarchical")
+  expect_identical(g$edge_prior, "Bernoulli")
+  expect_identical(g$pair_names[1:2], c("1-2", "1-3"))
+})
+
+test_that("hierarchical Beta-Bernoulli samples theta and coheres with it", {
+  a = 2
+  b = 4
+  g = sample_graph_prior(
+    p = 6, n_samples = 4000,
+    edge_prior = beta_bernoulli_prior(a, b), seed = 5
+  )
+  expect_length(g$theta, 4000L)
+  expect_true(all(g$theta > 0 & g$theta < 1))
+  expect_lt(abs(mean(g$theta) - a / (a + b)), 0.02)
+  expect_lt(abs(var(g$theta) - a * b / ((a + b)^2 * (a + b + 1))), 0.01)
+  # The graph marginal integrates theta out to a/(a + b), and per draw the
+  # edge mean tracks the sampled theta.
+  expect_lt(abs(mean(g$edge_indicators) - a / (a + b)), 0.02)
+  expect_gt(cor(rowMeans(g$edge_indicators), g$theta), 0.7)
+})
+
+test_that("hierarchical SBM separates within- from between-block edges", {
+  g = sample_graph_prior(
+    p = 8, n_samples = 1000,
+    edge_prior = sbm_prior(
+      alpha = 6, beta = 1, alpha_between = 1, beta_between = 6
+    ),
+    seed = 7
+  )
+  expect_identical(dim(g$allocations), c(1000L, 8L))
+  expect_true(all(g$allocations >= 1L))
+
+  pairs = which(upper.tri(matrix(0, 8, 8)), arr.ind = TRUE)
+  pairs = pairs[order(pairs[, 1L], pairs[, 2L]), , drop = FALSE]
+  within = matrix(
+    g$allocations[, pairs[, 1L]] == g$allocations[, pairs[, 2L]],
+    nrow = 1000L
+  )
+  ind = g$edge_indicators
+  expect_gt(mean(ind[within]), mean(ind[!within]) + 0.3)
+})
+
+test_that("conditioning on theta fixes the graph law", {
+  g = sample_graph_prior(
+    p = 6, n_samples = 2000,
+    edge_prior = beta_bernoulli_prior(1, 1), theta = 0.8, seed = 3
+  )
+  expect_lt(abs(mean(g$edge_indicators) - 0.8), 0.02)
+  expect_null(g$theta)
+})
+
+test_that("conditioning on (allocations, block_probs) fixes the SBM law", {
+  z = c(1L, 1L, 1L, 2L, 2L, 2L)
+  bp = matrix(c(0.9, 0.05, 0.05, 0.9), 2, 2)
+  g = sample_graph_prior(
+    p = 6, n_samples = 2000,
+    edge_prior = sbm_prior(), allocations = z, block_probs = bp, seed = 9
+  )
+  pairs = which(upper.tri(matrix(0, 6, 6)), arr.ind = TRUE)
+  pairs = pairs[order(pairs[, 1L], pairs[, 2L]), , drop = FALSE]
+  within = z[pairs[, 1L]] == z[pairs[, 2L]]
+  expect_lt(abs(mean(g$edge_indicators[, within]) - 0.9), 0.03)
+  expect_lt(abs(mean(g$edge_indicators[, !within]) - 0.05), 0.03)
+  expect_null(g$allocations)
+})
+
+test_that("conditioning arguments are validated", {
+  expect_error(
+    sample_graph_prior(
+      p = 5, n_samples = 10,
+      edge_prior = sbm_prior(), theta = 0.5
+    ),
+    "allocations"
+  )
+  expect_error(
+    sample_graph_prior(
+      p = 5, n_samples = 10,
+      edge_prior = sbm_prior(), allocations = rep(1L, 5)
+    ),
+    "block_probs"
+  )
+  expect_error(
+    sample_graph_prior(
+      p = 5, n_samples = 10,
+      edge_prior = bernoulli_prior(0.5),
+      allocations = rep(1L, 5),
+      block_probs = matrix(0.5, 1, 1)
+    ),
+    "sbm_prior"
+  )
+  expect_error(
+    sample_graph_prior(p = 5, n_samples = 10, theta = 1.5),
+    "theta"
+  )
+  expect_error(
+    sample_graph_prior(
+      p = 5, n_samples = 10,
+      edge_prior = sbm_prior(),
+      allocations = c(1L, 1L, 2L, 2L, 3L),
+      block_probs = matrix(0.5, 2, 2)
+    ),
+    "exceed"
+  )
+})
+
+test_that("the same seed reproduces the draw and restores the RNG state", {
+  set.seed(123)
+  before = runif(1)
+  set.seed(123)
+  g1 = sample_graph_prior(p = 6, n_samples = 50, seed = 42)
+  after = runif(1)
+  g2 = sample_graph_prior(p = 6, n_samples = 50, seed = 42)
+  expect_identical(g1$edge_indicators, g2$edge_indicators)
+  expect_identical(before, after)
+})
+
+test_that("the joint spec produces the Z(Gamma)-tilted law, not the prior", {
+  skip_on_cran()
+  # Mirrors the negative control in test-hier-zratio-identity: at
+  # tau^2 = 2 the joint marginal separates cleanly from Bernoulli(0.3)
+  # (about 0.22), while the hierarchical law sits at 0.3.
+  g = sample_graph_prior(
+    p = 6, n_samples = 4000,
+    edge_prior = bernoulli_prior(0.3), spec = "joint",
+    interaction_prior = normal_prior(scale = 0.5),
+    precision_scale_prior = gamma_prior(shape = 1, rate = 2),
+    n_warmup = 1500, seed = 7, verbose = FALSE
+  )
+  expect_identical(dim(g$edge_indicators), c(4000L, 15L))
+  expect_gt(abs(mean(g$edge_indicators) - 0.3), 0.05)
+
+  h = sample_graph_prior(
+    p = 6, n_samples = 4000,
+    edge_prior = bernoulli_prior(0.3), seed = 7
+  )
+  expect_lt(abs(mean(h$edge_indicators) - 0.3), 0.02)
+})
diff --git a/tests/testthat/test-sample-sbm-prior.R b/tests/testthat/test-sample-sbm-prior.R
new file mode 100644
index 00000000..0ac2124f
--- /dev/null
+++ b/tests/testthat/test-sample-sbm-prior.R
@@ -0,0 +1,59 @@
+# Tests for sample_sbm_prior(): ancestral draws from the MFM-SBM
+# edge-prior hyperprior.
+
+test_that("sample_sbm_prior returns coherent draws", {
+  d = sample_sbm_prior(p = 8, n_samples = 500, seed = 4)
+  expect_identical(dim(d$allocations), c(500L, 8L))
+  expect_identical(dim(d$pair_probability), c(500L, 28L))
+  expect_length(d$num_blocks, 500L)
+  expect_true(all(d$allocations >= 1L))
+  expect_true(all(d$pair_probability > 0 & d$pair_probability < 1))
+  expect_true(all(d$num_blocks >= 1L & d$num_blocks <= 8L))
+  # num_blocks counts the occupied blocks of each allocation.
+  expect_identical(
+    d$num_blocks,
+    apply(d$allocations, 1, function(z) length(unique(z)))
+  )
+})
+
+test_that("pair probabilities follow the block structure", {
+  # Dense within (Beta(8, 1)), sparse between (Beta(1, 8)).
+  d = sample_sbm_prior(
+    p = 8, n_samples = 500,
+    edge_prior = sbm_prior(
+      alpha = 8, beta = 1, alpha_between = 1, beta_between = 8
+    ),
+    seed = 2
+  )
+  pairs = which(upper.tri(matrix(0, 8, 8)), arr.ind = TRUE)
+  pairs = pairs[order(pairs[, 1L], pairs[, 2L]), , drop = FALSE]
+  within = matrix(
+    d$allocations[, pairs[, 1L]] == d$allocations[, pairs[, 2L]],
+    nrow = 500L
+  )
+  expect_gt(mean(d$pair_probability[within]), 0.8)
+  expect_lt(mean(d$pair_probability[!within]), 0.2)
+})
+
+test_that("lambda shifts the number of blocks", {
+  few = sample_sbm_prior(
+    p = 10, n_samples = 400,
+    edge_prior = sbm_prior(lambda = 0.1), seed = 6
+  )
+  many = sample_sbm_prior(
+    p = 10, n_samples = 400,
+    edge_prior = sbm_prior(lambda = 8), seed = 6
+  )
+  expect_lt(mean(few$num_blocks), mean(many$num_blocks))
+})
+
+test_that("input validation and reproducibility", {
+  expect_error(
+    sample_sbm_prior(p = 8, n_samples = 10, edge_prior = bernoulli_prior()),
+    "sbm_prior"
+  )
+  expect_error(sample_sbm_prior(p = 1, n_samples = 10), "'p'")
+  d1 = sample_sbm_prior(p = 6, n_samples = 20, seed = 9)
+  d2 = sample_sbm_prior(p = 6, n_samples = 20, seed = 9)
+  expect_identical(d1, d2)
+})

From 879169019a5f1d701a19312d0805c02c10af87ee Mon Sep 17 00:00:00 2001
From: Maarten Marsman <m.marsman@uva.nl>
Date: Mon, 6 Jul 2026 12:56:46 +0200
Subject: [PATCH 2/2] Drive the correction-table progress bar off
 display_progress

Gate the build progress bar on display_progress (the sampler's own
control) instead of bgms.verbose, and draw it for both serial and
parallel builds; a parallel build advances the bar as each batch of
forked chains completes. The build announcement stays on bgms.verbose.

Remove the top-level options(bgms.verbose = FALSE) from
helper-fixtures.R so devtools::load_all no longer silences fit-time
output in development sessions; setup.R still quiets the test run.

Document the cores contract in bgm(): sampling uses min(cores, chains),
while computations outside sampling (the correction-table build) use
all cores.
---
 NEWS.md                                 |  2 +-
 R/bgm.R                                 |  4 +-
 R/correction_tables.R                   | 71 ++++++++++++++++---------
 man/bgm.Rd                              |  4 +-
 tests/testthat/helper-fixtures.R        |  6 ++-
 tests/testthat/test-correction-tables.R | 16 ++++++
 6 files changed, 72 insertions(+), 31 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 467de928..3087e7db 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -14,7 +14,7 @@
 
 * `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), states that the build is one-time and cached, and reports progress (a progress bar on one core, a cell-and-core count otherwise).
+* 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.
diff --git a/R/bgm.R b/R/bgm.R
index 92f5d6de..e00f81e8 100644
--- a/R/bgm.R
+++ b/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.
diff --git a/R/correction_tables.R b/R/correction_tables.R
index 3eb8fd55..18d94800 100644
--- a/R/correction_tables.R
+++ b/R/correction_tables.R
@@ -253,7 +253,7 @@ sweep_prior_edge_density = function(p, theta, delta,
                                     n_seeds = 3L,
                                     update_method = "gibbs",
                                     cores = 1L, base_seed = 1L,
-                                    verbose = FALSE) {
+                                    show_progress = FALSE) {
   cores = normalize_builder_cores(cores)
   num_pairs = p * (p - 1) / 2
   cells = expand.grid(theta = theta, seed = seq_len(n_seeds))
@@ -272,31 +272,42 @@ sweep_prior_edge_density = function(p, theta, delta,
     mean(draws$edge_indicators)
   }
 
-  # A progress bar is drawn only in an interactive session; in batch runs
-  # (scripts, R CMD check, the test suite) the one-line build announcement
-  # in ggm_correction_table() is the only output.
-  show_bar = verbose && interactive() && cores == 1L
-  edens_cells = if(cores > 1L) {
-    if(verbose) {
-      message(sprintf(
-        "Sweeping %d prior cells on %d cores.", nrow(cells), cores
-      ))
-    }
-    unlist(parallel::mclapply(
-      seq_len(nrow(cells)), one_cell,
-      mc.cores = cores, mc.preschedule = FALSE
-    ))
-  } else if(show_bar) {
-    pb = new_correction_progress(nrow(cells))
+  # 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)
-    vapply(seq_len(nrow(cells)), function(k) {
+  }
+  edens_cells = if(cores > 1L) {
+    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(n_cells), function(k) {
       v = one_cell(k)
-      pb$update(k)
+      if(draw_bar) {
+        pb$update(k)
+      }
       v
     }, numeric(1))
-  } else {
-    vapply(seq_len(nrow(cells)), one_cell, numeric(1))
   }
   if(anyNA(edens_cells)) {
     stop("Correction-table sweep: a prior chain failed.")
@@ -347,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, verbose = FALSE
+  cores = 1L, base_seed = 1L, show_progress = FALSE
 ) {
   update_method = match.arg(update_method)
   if(is.null(delta)) {
@@ -361,7 +372,7 @@ build_ggm_correction_table = function(
     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,
-    verbose = verbose
+    show_progress = show_progress
   )
 
   table = correction_table_from_edens(
@@ -477,13 +488,20 @@ ggm_edge_prior_correction = function(prior, sampler, num_variables,
   interaction_prior = correction_interaction_prior(prior)
   precision_scale_prior = correction_scale_prior(prior)
 
+  # 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,
-    verbose = isTRUE(sampler$verbose)
+    verbose = isTRUE(sampler$verbose),
+    show_progress = show_progress
   )
   if(identical(prior$edge_prior, "Stochastic-Block") &&
     is.null(table$fprime)) {
@@ -521,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, verbose = FALSE
+  cores = 1L, base_seed = 1L, refresh = FALSE, verbose = FALSE,
+  show_progress = FALSE
 ) {
   update_method = match.arg(update_method)
   if(is.null(delta)) {
@@ -564,7 +583,7 @@ ggm_correction_table = function(
     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, verbose = verbose
+    cores = cores, base_seed = base_seed, show_progress = show_progress
   )
 
   if(use_cache) {
diff --git a/man/bgm.Rd b/man/bgm.Rd
index 9b9ea71f..c2d70311 100644
--- a/man/bgm.Rd
+++ b/man/bgm.Rd
@@ -259,7 +259,9 @@ matrix. Default: \code{TRUE}.}
 \item{chains}{Integer. Number of parallel chains to run. Default: \code{4}.}
 
 \item{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()}.}
 
 \item{display_progress}{Character. Controls progress reporting during
 sampling. Options: \code{"per-chain"} (separate bar per chain),
diff --git a/tests/testthat/helper-fixtures.R b/tests/testthat/helper-fixtures.R
index 82f9ba73..f35dff31 100644
--- a/tests/testthat/helper-fixtures.R
+++ b/tests/testthat/helper-fixtures.R
@@ -49,8 +49,10 @@
 # Ensure bgms package is loaded
 library(bgms)
 
-# Suppress informational messages during tests
-options(bgms.verbose = FALSE)
+# Advisory output is quieted for the test run in setup.R, not here. Helper files
+# are sourced by devtools::load_all() (setup files are not), so setting the
+# option here would leak bgms.verbose = FALSE into interactive development
+# sessions and silence fit-time messages and progress bars.
 
 # ------------------------------------------------------------------------------
 # 1. Session-Cached Model Fixtures
diff --git a/tests/testthat/test-correction-tables.R b/tests/testthat/test-correction-tables.R
index f6228e3b..8813b46a 100644
--- a/tests/testthat/test-correction-tables.R
+++ b/tests/testthat/test-correction-tables.R
@@ -118,6 +118,22 @@ test_that("the progress bar renders the label, counts, and percentage", {
   expect_match(full, "120/120 (100.0%)", fixed = TRUE)
 })
 
+test_that("the parallel sweep matches the serial sweep cell for cell", {
+  skip_on_cran()
+  skip_on_os("windows")
+
+  args = list(
+    p = 4, theta = c(0.2, 0.5, 0.8), delta = 0.5 * log(4),
+    interaction_prior = cauchy_prior(scale = 2.5),
+    precision_scale_prior = gamma_prior(shape = 1, eta = 1),
+    n_samples = 200L, n_warmup = 100L, n_seeds = 2L, update_method = "gibbs"
+  )
+  serial = do.call(sweep_prior_edge_density, c(args, cores = 1L))
+  parallel = do.call(sweep_prior_edge_density, c(args, cores = 2L))
+
+  expect_equal(parallel$edens_raw, serial$edens_raw)
+})
+
 test_that("the build announces itself once; a cache hit is silent", {
   cache_dir = file.path(tempdir(), "bgms-ctable-msg-test")
   unlink(cache_dir, recursive = TRUE)