# Trust-gauge reference route (block_reference_logR): the block-local exact
# Monte-Carlo reference log R_e that the in-chain gauge compares against the
# deployed log J. These drive the engine reference in isolation.
#
# The p = 16 biased-fit detector blocks run in the weekly certification tier
# (T2, BGMS_RUN_CERTIFICATION) BY CONTRACT: maintainer ruling 2026-08-03
# (F-103 (iii)) — their statistic is a per-fit Monte-Carlo quantity with real
# seed-to-seed dispersion, so they are distributional assertions and belong
# with the Monte-Carlo machinery whatever their runtime (helper-tiers.R).
#
# harm_threshold history (F-103): first parked 2026-08-01 when the Linux
# nightly caught harm_pred 0.01071 against the then-threshold 0.01 (run
# 30715557912, report 12). Report 13 measured the healthy negative control
# across 12 seeds: 0.0004-0.0125, threshold inside the spread. Report 24 ran
# the sweeps program: the spread is the fit's, not the audit's (between-seed
# sd flat across an 8x audit increase, per-seed r = 0.94, the top seed
# CONVERGES to ~0.0108), so no sweep count separates it — the default sweep
# count is unchanged. Resolution 2026-08-03 (F-103 (i), maintainer):
# harm_threshold raised 0.01 -> 0.02 — above the healthy 12-seed max
# (0.01253 observed, ~0.0108 converged; 1.6-1.9x margin) and below the
# biased-kernel fixture (0.0431 at the pinned seed, seed-stable across
# audit sizes; 2.2x margin), near their geometric midpoint. The biased
# kernel's per-fit values 0.0012-0.0431 OVERLAP the healthy range at some
# seeds: a per-realization projection can be honestly small on a coarse
# kernel, and the overlap zone is exactly where an alarm could not be told
# from healthy dispersion. Whether the projection itself should be
# recalibrated is backlogged (F-103 (ii), plans/software/flagged-issues.md).
# Tables: reports/13, 24 + the F-103 row in dev/review-2026-08/FINDINGS.md.

test_that("the block reference is finite and lands near the deployed ratio", {
  skip_on_cran()
  p = 10
  delta = 0.5 * log(p)
  eta = 2
  zc = bgms:::zratio_constants(delta, eta = eta)
  set.seed(3)
  G = matrix(0L, p, p)
  ut = upper.tri(G)
  G[ut] = rbinom(sum(ut), 1, 0.6)
  G = G + t(G)
  diag(G) = 1L

  # First edge whose mediating block is non-trivial (m >= 2).
  edge = NULL
  for(i in seq_len(p - 1)) {
    for(j in (i + 1):p) {
      r = bgms:::zratio_test_reference(
        G, i, j, zc$addc, zc$tg, zc$ihat, zc$ghat, zc$wt, zc$psi0,
        delta, eta, 120L, 30L, 7L, FALSE
      )
      if(isTRUE(r$valid) && r$m >= 2) {
        edge = c(i, j)
        break
      }
    }
    if(!is.null(edge)) break
  }
  expect_false(is.null(edge))

  r = bgms:::zratio_test_reference(
    G, edge[1], edge[2], zc$addc, zc$tg, zc$ihat, zc$ghat, zc$wt, zc$psi0,
    delta, eta, 240L, 30L, 7L, FALSE
  )
  expect_true(r$ok)
  expect_true(is.finite(r$logR))
  expect_true(is.finite(r$mcse) && r$mcse > 0)
  # The deployed corrected ratio is a validated approximation of the exact
  # reference on an entangled edge: they agree to within a few MCSE.
  expect_lt(abs(r$logR - r$log_zratio), 0.05)
})

test_that("the reference MCSE shrinks with the draw count", {
  skip_on_cran()
  p = 10
  delta = 0.5 * log(p)
  eta = 2
  zc = bgms:::zratio_constants(delta, eta = eta)
  set.seed(3)
  G = matrix(0L, p, p)
  ut = upper.tri(G)
  G[ut] = rbinom(sum(ut), 1, 0.6)
  G = G + t(G)
  diag(G) = 1L

  # A high-degree endpoint pair has an entangled block.
  small = bgms:::zratio_test_reference(
    G, 9, 10, zc$addc, zc$tg, zc$ihat, zc$ghat, zc$wt, zc$psi0,
    delta, eta, 60L, 30L, 11L, FALSE
  )
  large = bgms:::zratio_test_reference(
    G, 9, 10, zc$addc, zc$tg, zc$ihat, zc$ghat, zc$wt, zc$psi0,
    delta, eta, 960L, 30L, 11L, FALSE
  )
  expect_true(small$valid && large$valid)
  # Two estimates from independent streams agree within combined MCSE.
  expect_lt(
    abs(small$logR - large$logR),
    5 * sqrt(small$mcse^2 + large$mcse^2)
  )
  expect_lt(large$mcse, small$mcse)
})

# Harm channel: standardized se (projected inclusion-probability distortion)
# on top of the flip-rate channel.

test_that("the harm channel computes the documented statistics", {
  gauge = list(
    flip_rate = 0.001, noise_floor = 1e-4, se_mean = 0.05, se_sd = 0.02,
    se_mcse = 0.002, n_ent = 100, n_ref = 25, n_capped = 0
  )
  chains = list(list(zratio = list(gauge = gauge)))

  # Beta-Bernoulli, exchangeable inclusion probabilities: the linearized
  # feedback gain is E * m / (theta (1 - theta) (a + b + E)).
  a = 9
  b = 1
  E = 190
  th = 0.9
  pip = rep(th, E)
  s = summarize_zratio_gauge(
    chains,
    verbose = FALSE,
    harm_inputs = list(pip = list(pip), a = a, b = b)
  )
  pc = s$per_chain
  expect_equal(pc$se_se, sqrt(0.002^2 + 0.02^2 / 25))
  m_bar = th * (1 - th)
  theta_hat = (a + E * th) / (a + b + E)
  gain = E * m_bar / (theta_hat * (1 - theta_hat) * (a + b + E))
  expect_equal(pc$amplification, 1 / (1 - min(gain, 0.98)))
  expect_equal(pc$harm_pred, abs(0.05) * m_bar * pc$amplification)
  # 0.05 nats, sensitivity 0.09, ~18x amplification: far over tolerance,
  # and resolved (|se_mean| = 0.05 > 2 * se_se).
  expect_true(pc$harm_flag)
  expect_true(s$flagged)

  # Fixed inclusion probability: no feedback, amplification 1.
  s_bern = summarize_zratio_gauge(
    chains,
    verbose = FALSE,
    harm_inputs = list(pip = list(pip), a = NULL, b = NULL)
  )
  expect_equal(s_bern$per_chain$amplification, 1)
  expect_equal(s_bern$per_chain$harm_pred, 0.05 * m_bar)

  # Unresolved se (large reference noise) must not flag regardless of size.
  gauge_noisy = gauge
  gauge_noisy$se_mcse = 0.2
  s_noisy = summarize_zratio_gauge(
    list(list(zratio = list(gauge = gauge_noisy))),
    verbose = FALSE,
    harm_inputs = list(pip = list(pip), a = a, b = b)
  )
  expect_false(s_noisy$per_chain$harm_flag)

  # No harm inputs: channel reports NA and never flags.
  s_off = summarize_zratio_gauge(chains, verbose = FALSE)
  expect_true(is.na(s_off$per_chain$harm_pred))
  expect_false(s_off$per_chain$harm_flag)

  # Unsupported edge priors disable the channel at the input builder.
  expect_null(zratio_harm_inputs(list(pip), "Stochastic-Block", a = 1, b = 1))
})

test_that("a flagged chain prints the remediation ladder in order", {
  withr::local_options(bgms.verbose = TRUE)
  gauge = list(
    flip_rate = 0.05, noise_floor = 1e-4, se_mean = 0.05, se_sd = 0.02,
    se_mcse = 0.002, n_ent = 9900, n_ref = 50, n_capped = 9850,
    pair_m = c(4L, 38L)
  )
  out = capture.output(
    summarize_zratio_gauge(list(list(zratio = list(gauge = gauge))))
  )
  txt = paste(out, collapse = " ")

  # Rung 1: what was measured, on which blocks, out of how many moves.
  expect_match(txt, "5.0% of edge-toggle decisions")
  expect_match(txt, "audited 50 of 9900 non-trivial edge moves")
  expect_match(txt, "mediating blocks 4-38 variables")
  # Rung 2: resolve the signal before changing the model.
  expect_match(txt, "Raise options\\(bgms.zratio_gauge_sweeps\\)")
  # Rung 3: the joint specification, named as a different inferential target.
  expect_match(txt, "targets a different model")
  expect_match(txt, "reweighted by the per-graph normalizer")
  # The ladder is ordered: measurement, then more sweeps, then the joint spec.
  expect_lt(
    regexpr("Raise options", txt, fixed = TRUE),
    regexpr("precision_graph_prior = \"joint\"", txt, fixed = TRUE)
  )
  expect_lt(
    regexpr("audited 50 of", txt, fixed = TRUE),
    regexpr("Raise options", txt, fixed = TRUE)
  )
})

test_that("a clean chain prints nothing", {
  withr::local_options(bgms.verbose = TRUE)
  gauge = list(
    flip_rate = 0.0001, noise_floor = 1e-4, se_mean = 0, se_sd = 0.02,
    se_mcse = 0.002, n_ent = 100, n_ref = 25, n_capped = 0, pair_m = c(3L, 5L)
  )
  s = summarize_zratio_gauge(list(list(zratio = list(gauge = gauge))))
  expect_false(s$flagged)
  expect_equal(
    capture.output(
      summarize_zratio_gauge(list(list(zratio = list(gauge = gauge))))
    ),
    character(0)
  )
  expect_equal(s$per_chain$block_lo, 3L)
  expect_equal(s$per_chain$block_hi, 5L)
})

# A 16-variable prior draw under a dense-leaning Beta-Bernoulli prior with no
# data: the feedback-amplified regime where the flip rate stays quiet but the
# projected distortion is first-order. The Gamma-diagonal shape selects which
# kernel the engine reaches, so both sides of the surface's deployment fence
# run the same fixture.
biased_evidence_free_fit = function(shape) {
  sample_ggm_prior(
    p = 16L, n_samples = 1200L, n_warmup = 500L,
    interaction_prior = normal_prior(scale = 0.5),
    precision_scale_prior = gamma_prior(shape = shape, rate = 6),
    spec = "hierarchical", edge_prior = beta_bernoulli_prior(9, 1),
    update_method = "gibbs",
    zratio_diagnostics = TRUE, seed = 7L, verbose = FALSE
  )
}

test_that("a known-biased evidence-free fit fires the harm channel", {
  skip_on_cran()
  skip_unless_certification()
  # Below the surface's deployment fence (.zratio_surface_shape_lo = 0.5) the
  # engine falls back to the additive-counts saddle, and that coarse kernel is
  # what this channel exists to police. The Gamma-shape constants are
  # themselves unscored there, which is what the cell warns about; the fixture
  # is chosen for the kernel it reaches, not as a certified cell.
  f = suppressWarnings(biased_evidence_free_fit(0.4))
  pc = f$zratio_diagnostics$per_chain
  expect_true(is.finite(pc$se_mcse) && pc$se_mcse > 0)
  expect_true(is.finite(pc$harm_pred))
  expect_gt(pc$amplification, 5)
  expect_gt(pc$harm_pred, f$zratio_diagnostics$harm_threshold)
  expect_true(pc$harm_flag)
  # Rung 1 stays quiet: the coarse kernel is invisible on the flip rate, which
  # is why the harm channel is a separate one.
  expect_lt(pc$flip_rate, 0.01)
})

test_that("the deployed surface holds that fixture under the harm threshold", {
  skip_on_cran()
  skip_unless_certification()
  # The same fixture at shape 2, which the surface's deployment range [0.5, 10]
  # covers. The amplification is a property of the prior and the fit, so it
  # stays first-order; the accurate kernel cuts the projected distortion about
  # five-fold and the flag correctly stays down. This is the harm channel's
  # negative control on a cell that is genuinely at risk.
  f = biased_evidence_free_fit(2)
  pc = f$zratio_diagnostics$per_chain
  expect_gt(pc$amplification, 5)
  expect_lt(pc$harm_pred, f$zratio_diagnostics$harm_threshold)
  expect_false(pc$harm_flag)
})

test_that("the harm channel weights errors by per-edge sensitivity", {
  # Per-pair audit stream present: harm_pred = |mean(m_e s_e)| * A with a
  # cluster-robust (per audited edge) resolution gate. Edge (0,1) is audited
  # twice to exercise the clustering.
  gauge = list(
    flip_rate = 0.001, noise_floor = 1e-4, se_mean = mean(c(0.05, 0.03, 0.05)),
    se_sd = 0.01, se_mcse = 0.002, n_ent = 100, n_ref = 3, n_capped = 0,
    pair_i = c(0L, 0L, 0L), pair_j = c(1L, 2L, 1L),
    pair_se = c(0.05, 0.03, 0.05), pair_mcse = c(0.002, 0.002, 0.002)
  )
  chains = list(list(zratio = list(gauge = gauge)))
  a = 9
  b = 1
  E = 190
  th = 0.9
  pip = rep(th, E)
  s = summarize_zratio_gauge(
    chains,
    verbose = FALSE,
    harm_inputs = list(pip = list(pip), a = a, b = b)
  )
  pc = s$per_chain

  m = th * (1 - th)
  theta_hat = (a + E * th) / (a + b + E)
  gain = E * m / (theta_hat * (1 - theta_hat) * (a + b + E))
  A = 1 / (1 - min(gain, 0.98))
  expect_equal(pc$kappa, m * A)
  x = m * c(0.05, 0.03, 0.05)
  expect_equal(pc$harm_pred, abs(mean(x)) * A)
  # Resolution gate: cluster-robust spread + reference noise.
  cr = ((sum(x[c(1, 3)]) - 2 * mean(x))^2 + (x[2] - mean(x))^2) / 9
  noise2 = sum((m * 0.002)^2 * 3) / 9
  expect_true(abs(mean(x)) > 2 * sqrt(cr + noise2))
  expect_true(pc$harm_flag)

  # Heterogeneous sensitivities: errors on pinned edges must not count.
  pip2 = rep(th, E)
  pip2[1] = 0.999 # edge (0,1) pinned: m_e ~ 0
  s2 = summarize_zratio_gauge(
    chains,
    verbose = FALSE,
    harm_inputs = list(pip = list(pip2), a = a, b = b)
  )
  m1 = 0.999 * (1 - 0.999)
  m2 = th * (1 - th)
  x2 = c(m1 * 0.05, m2 * 0.03, m1 * 0.05)
  expect_equal(
    s2$per_chain$harm_pred,
    abs(mean(x2)) * s2$per_chain$amplification
  )
  expect_lt(s2$per_chain$harm_pred, pc$harm_pred)
})

test_that("the extrapolation notice is graceful, gated, and back-compatible", {
  mk = function(nx, mx, np) {
    list(zratio = list(counters = c(
      n_hit = 0, n_miss = 0, n_pred = np,
      n_add = 0, cache_size = 0, n_extrap = nx, max_extrap_size = mx
    )))
  }
  # No extrapolation -> silent, returns FALSE.
  expect_silent(res0 <- bgms:::zratio_extrapolation_notice(list(mk(0, 0, 100))))
  expect_false(res0)
  # Extrapolation -> one graceful message reporting the largest block size.
  # These counters carry no warmup/retained split, so the notice reports the
  # whole-run share rather than claiming a phase for them.
  expect_message(
    bgms:::zratio_extrapolation_notice(list(mk(120, 73, 4000), mk(80, 66, 4000))),
    "beyond its anchored size range \\(largest 73"
  )
  # Old-format chains without the counter -> silent (back-compatible).
  old = list(zratio = list(counters = c(n_hit = 0, n_pred = 10)))
  expect_silent(res2 <- bgms:::zratio_extrapolation_notice(list(old)))
  expect_false(res2)
})

test_that("the harm channel maps audit records onto the correct edge", {
  # Regression for the %/% precedence bug in the edge-index formula: for the
  # 0-based pair (2, 3) at q = 7 the 1-based upper-triangle index is 12; the
  # unparenthesized form floor-divided (2q - i0 - 1) first and landed on 11.
  # pip isolates the weight on index 12, so a mis-mapped record zeroes the
  # sensitivity weight and the predictor.
  g = list(
    flip_rate = 0, noise_floor = 0, se_mean = 0.2, se_sd = 0.1,
    se_mcse = 1e-6, n_ref = 1L, n_ent = 1L, n_capped = 0L,
    pair_i = 2L, pair_j = 3L, pair_se = 0.2, pair_mcse = 1e-6
  )
  pip = numeric(21)
  pip[12] = 0.5
  res = summarize_zratio_gauge(
    list(list(zratio = list(gauge = g))),
    verbose = FALSE, harm_inputs = list(pip = list(pip))
  )
  # m = 0.5 * 0.5 on the audited edge, se = 0.2, amplification 1 (no a/b):
  # harm_pred = 0.25 * 0.2 = 0.05, resolved and above the 0.02 threshold.
  expect_equal(res$per_chain$harm_pred, 0.05, tolerance = 1e-12)
  expect_true(res$per_chain$harm_flag)
})

test_that("harm inputs stay aligned when a chain has no gauge output", {
  # Chain 1 carries no gauge block (e.g. an interrupt during its sweeps); the
  # summary must index harm_inputs$pip by the ORIGINAL chain position, not the
  # position after filtering, and report the original chain number.
  g = list(
    flip_rate = 0, noise_floor = 0, se_mean = 0.2, se_sd = 0.1,
    se_mcse = 1e-6, n_ref = 1L, n_ent = 1L, n_capped = 0L,
    pair_i = 2L, pair_j = 3L, pair_se = 0.2, pair_mcse = 1e-6
  )
  pip2 = numeric(21)
  pip2[12] = 0.5
  res = summarize_zratio_gauge(
    list(list(), list(zratio = list(gauge = g))),
    verbose = FALSE,
    harm_inputs = list(pip = list(numeric(21), pip2))
  )
  expect_equal(nrow(res$per_chain), 1L)
  expect_equal(res$per_chain$chain, 2L)
  expect_equal(res$per_chain$harm_pred, 0.05, tolerance = 1e-12)
})

test_that("a mixed hierarchical fit reports harm, not NA (F-022)", {
  skip_on_cran()
  # The mixed builder used to call summarize_zratio_gauge() without
  # harm_inputs, so harm_pred, amplification and kappa were permanently NA on
  # every mixed fit -- the gauge's second channel was dead on that path. The
  # fixture is dense-leaning and evidence-free on the continuous block so the
  # audit stream is non-empty in a couple of seconds.
  withr::local_options(bgms.zratio_gauge_sweeps = 2L, bgms.verbose = FALSE)
  set.seed(11)
  n = 50
  x = cbind(
    matrix(sample(0:2, n * 2, replace = TRUE), n, 2),
    matrix(rnorm(n * 8), n, 8)
  )
  fit = bgm(
    x, variable_type = c(rep("ordinal", 2), rep("continuous", 8)),
    interaction_prior = normal_prior(scale = 0.5),
    precision_scale_prior = gamma_prior(shape = 1, rate = 6),
    edge_prior = beta_bernoulli_prior(9, 1),
    precision_graph_prior = "hierarchical",
    iter = 100, warmup = 150, update_method = "adaptive-metropolis",
    chains = 1, cores = 1, seed = 5,
    display_progress = "none", verbose = FALSE
  )

  pc = fit$zratio_diag$per_chain
  expect_gt(pc$n_ref, 0L)
  expect_true(is.finite(pc$harm_pred))
  expect_true(is.finite(pc$amplification))
  expect_true(is.finite(pc$kappa))
  # The Beta-Bernoulli feedback is real on this fixture, so the wiring is
  # visibly doing something rather than passing a degenerate pool through.
  expect_gt(pc$amplification, 1)
  expect_false(is.na(pc$harm_flag))
})

test_that("the harm pool can be wider than the audited block", {
  # A mixed fit audits the continuous-continuous edges only, while a
  # Beta-Bernoulli theta is drawn from every edge of the graph. pool_pip
  # carries that wider pool; without it the feedback gain would be computed
  # from the audited block alone and understate the amplification.
  g = list(
    flip_rate = 0, noise_floor = 0, se_mean = 0.2, se_sd = 0.1,
    se_mcse = 1e-6, n_ref = 1L, n_ent = 1L, n_capped = 0L,
    pair_i = 0L, pair_j = 1L, pair_se = 0.2, pair_mcse = 1e-6
  )
  chains = list(list(zratio = list(gauge = g)))
  a = 9
  b = 1
  th = 0.9
  audited = rep(th, 10L) # 5 continuous variables
  pool = rep(th, 45L) # plus 5 discrete: the whole mixed graph

  narrow = summarize_zratio_gauge(
    chains,
    verbose = FALSE, harm_inputs = list(pip = list(audited), a = a, b = b)
  )
  wide = summarize_zratio_gauge(
    chains,
    verbose = FALSE,
    harm_inputs = list(
      pip = list(audited), a = a, b = b, pool_pip = list(pool)
    )
  )

  gain_of = function(E) {
    m = th * (1 - th)
    theta_hat = (a + E * th) / (a + b + E)
    E * m / (theta_hat * (1 - theta_hat) * (a + b + E))
  }
  amp_of = function(E) 1 / (1 - min(gain_of(E), 0.98))
  expect_equal(narrow$per_chain$amplification, amp_of(10L))
  expect_equal(wide$per_chain$amplification, amp_of(45L))
  expect_gt(wide$per_chain$amplification, narrow$per_chain$amplification)

  # The numerator is still weighted by the AUDITED edge's sensitivity, so the
  # two differ by the amplification alone.
  expect_equal(
    wide$per_chain$harm_pred / narrow$per_chain$harm_pred,
    amp_of(45L) / amp_of(10L)
  )
  # Omitting pool_pip must leave the GGM path exactly as it was.
  expect_equal(
    narrow$per_chain,
    summarize_zratio_gauge(
      chains,
      verbose = FALSE,
      harm_inputs = list(pip = list(audited), a = a, b = b, pool_pip = NULL)
    )$per_chain
  )
})

test_that("the prior sampler resolves gauge sweeps from the option (F-103)", {
  skip_on_cran()
  # sample_ggm_prior() hardwired 2 sweeps while the deployed path resolved
  # options(bgms.zratio_gauge_sweeps) through zratio_gauge_sweeps(). The two
  # now share one resolution, so an audit measured on the prior chain describes
  # the audit a deployed fit performs. zratio_diagnostics stays the on/off
  # switch: it is not the precision.
  prior_fit = function(...) {
    sample_ggm_prior(
      p = 10L, n_samples = 200L, n_warmup = 100L,
      interaction_prior = normal_prior(scale = 0.5),
      precision_scale_prior = gamma_prior(shape = 2, rate = 6),
      spec = "hierarchical", edge_prior = beta_bernoulli_prior(9, 1),
      update_method = "gibbs", seed = 3L, verbose = FALSE, ...
    )
  }
  n_ref_at = function(sweeps) {
    withr::local_options(bgms.zratio_gauge_sweeps = sweeps)
    prior_fit(zratio_diagnostics = TRUE)$zratio_diagnostics$per_chain$n_ref
  }

  # The gauge references at most a fixed number of pairs per sweep, so a
  # saturated audit grows one cap per sweep: the option, and nothing else,
  # sets the audit size.
  expect_equal(n_ref_at(2L), 2L * n_ref_at(1L))
  expect_equal(n_ref_at(6L), 6L * n_ref_at(1L))

  # The on/off switch still wins over the option.
  withr::local_options(bgms.zratio_gauge_sweeps = 6L)
  expect_null(prior_fit(zratio_diagnostics = FALSE)$zratio_diagnostics)
})