# Tests for the hierarchical-spec per-edge Z-ratio engine: parity against a
# pre-generated exact-reference fixture (zratio_reference.rds), the fit-time
# constant builders, and engine invariants.

skip_heavy_guard_on_cran() # CI-only on CRAN's clock; see helper-tiers.R

fixture_path = testthat::test_path("fixtures", "zratio_reference.rds")

test_that("engine reproduces the reference log Z-ratios (all cells/variants)", {
  skip_if_not(file.exists(fixture_path), "fixtures/zratio_reference.rds not in this build (the tarball excludes tests/testthat/fixtures/) - regenerate with certification/regenerate_zratio_reference.R")
  fx = readRDS(fixture_path)
  for(cell in fx) {
    addc13 = c(cell$addc6, cell$fc_coef, 1)
    addc23 = c(addc13, cell$hull)
    for(gn in names(cell$graphs)) {
      g = cell$graphs[[gn]]
      eu = which(upper.tri(g), arr.ind = TRUE)
      edges = cbind(eu[, 1], eu[, 2])
      for(variant in c("base", "direct", "clamp")) {
        addc = switch(variant,
          base = cell$addc6,
          direct = addc13,
          clamp = addc23
        )
        res = zratio_test_eval(
          g, edges, addc, cell$tg, cell$ihat, cell$ghat, cell$wt, cell$psi0
        )
        # The OLS correction slots (addc[6..22]) are inert: the engine reads
        # only the 6 additive constants, so passing the extended blocks
        # (direct/clamp) yields the same additive saddle as the base block.
        ref = cell$evals[[paste(gn, "base", sep = "_")]]
        expect_equal(
          as.numeric(res$log_zratio), unname(ref),
          tolerance = 1e-12,
          label = sprintf(
            "delta=%.3f sigma=%g beta=%g %s %s",
            cell$delta, cell$sigma, cell$beta, gn, variant
          )
        )
      }
    }
  }
})

test_that("fit-time constant builders match the reference builders", {
  skip_if_not(file.exists(fixture_path), "fixtures/zratio_reference.rds not in this build (the tarball excludes tests/testthat/fixtures/) - regenerate with certification/regenerate_zratio_reference.R")
  fx = readRDS(fixture_path)
  # The builders run in the standardized cell (sigma = 1), so match against
  # the two sigma = 1 reference cells; eta = sigma * beta = beta there.
  for(cell in fx[c(1, 4)]) {
    zc = bgms:::zratio_constants(cell$delta, cell$sigma * cell$beta)
    expect_equal(zc$addc, cell$addc6, tolerance = 1e-8)
    expect_equal(zc$psi0, cell$psi0, tolerance = 1e-8)
    expect_equal(zc$tg, cell$tg, tolerance = 1e-12)
    expect_equal(zc$wt, cell$wt, tolerance = 1e-12)
    expect_equal(zc$ihat, cell$ihat, tolerance = 1e-8)
    expect_equal(zc$ghat, cell$ghat, tolerance = 1e-8)
  }
})

test_that("the pair table and the saddle integration assert a common domain", {
  # The pair integrals are interpolated, and a read past the table's end clamps
  # to its edge value: a plateau where the integrand should decay. The saddle
  # grid therefore refuses a table it would have to read past.
  short = bgms:::zratio_pair_integrals(delta = 1, sigma = 1, beta = 1, cmax = 8)
  expect_error(
    bgms:::zratio_saddle_grid(short, Cmax = 40),
    "pair table covers"
  )
  # The shipped pairing satisfies it, checked on the domains rather than by
  # running the transform.
  wide = bgms:::zratio_pair_integrals(delta = 1, sigma = 1, beta = 1)
  expect_gte(max(wide$cg), eval(formals(bgms:::zratio_saddle_grid)$Cmax))
})

test_that("Gauss quadrature nodes integrate known moments exactly", {
  gl = bgms:::zratio_gauss_quad(24, "laguerre")
  expect_equal(sum(gl$weights), 1, tolerance = 1e-12) # int e^-x
  expect_equal(sum(gl$weights * gl$nodes^3), 6, tolerance = 1e-9) # Gamma(4)
  gh = bgms:::zratio_gauss_quad(24, "hermite")
  expect_equal(sum(gh$weights), sqrt(pi), tolerance = 1e-12)
  expect_equal(sum(gh$weights * gh$nodes^2), sqrt(pi) / 2, tolerance = 1e-9)
  lg = bgms:::zratio_gauss_quad(24, "legendre")
  expect_equal(sum(lg$weights), 2, tolerance = 1e-12)
  expect_equal(sum(lg$weights * lg$nodes^2), 2 / 3, tolerance = 1e-12)
})

test_that("I_spike matches its delta = 0 Bessel identity", {
  cg = seq(0.2, 8, length.out = 12)
  # At delta = 0, beta = 0.5: I_spike(c) = 4 c K_1(c).
  expect_equal(
    bgms:::zratio_ispike(cg, delta = 0, beta = 0.5),
    4 * cg * besselK(cg, 1),
    tolerance = 1e-10
  )
})

test_that("engine invariants: state-invariance, isolated edge, precompute", {
  skip_if_not(file.exists(fixture_path), "fixtures/zratio_reference.rds not in this build (the tarball excludes tests/testthat/fixtures/) - regenerate with certification/regenerate_zratio_reference.R")
  cell = readRDS(fixture_path)[[2]]
  g = cell$graphs$er3
  edges = cbind(1L, 2L)

  # The ratio is invariant to the toggled edge's own state.
  g_on = g
  g_on[1, 2] = g_on[2, 1] = 1L
  g_off = g
  g_off[1, 2] = g_off[2, 1] = 0L
  r_on = zratio_test_eval(
    g_on, edges, cell$addc6, cell$tg, cell$ihat, cell$ghat, cell$wt, cell$psi0
  )
  r_off = zratio_test_eval(
    g_off, edges, cell$addc6, cell$tg, cell$ihat, cell$ghat, cell$wt, cell$psi0
  )
  expect_identical(r_on$log_zratio, r_off$log_zratio)

  # No mediating structure: the isolated-edge ratio psi0.
  g_empty = diag(1L, 10)
  r_iso = zratio_test_eval(
    g_empty, edges, cell$addc6, cell$tg, cell$ihat, cell$ghat, cell$wt,
    cell$psi0
  )
  expect_equal(r_iso$log_zratio[1], log(cell$psi0), tolerance = 1e-14)

  # Preloading the count table removes every cache miss and preserves values.
  eu = which(upper.tri(g), arr.ind = TRUE)
  all_edges = cbind(eu[, 1], eu[, 2])
  lazy = zratio_test_eval(
    g, all_edges, cell$addc6, cell$tg, cell$ihat, cell$ghat, cell$wt,
    cell$psi0
  )
  pre = zratio_test_precompute(
    g, all_edges, cell$addc6, cell$tg, cell$ihat, cell$ghat, cell$wt,
    cell$psi0, 10L, 30L
  )
  expect_equal(as.numeric(pre$log_zratio), as.numeric(lazy$log_zratio), tolerance = 1e-14)
  expect_true(pre$n_miss == 0)
})