tests/testthat/test-within-cell-box-uniform.R

In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models

# The WITHIN-CELL construction a reported hyperparameter interval is read with
# (gcol33/tulpa#357).
#
# The outer grid's weights say how much mass each cell holds. They do not say
# how it is spread INSIDE the cell, and a quantile needs both. The shipped
# `chord` read puts the cumulative MID-mass at each cell COORDINATE and
# interpolates between coordinates; `box_uniform` puts the cumulative FULL mass
# at each cell EDGE and interpolates between edges -- the same masses over the
# same boxes with the knots moved half a cell, which is a whole order of
# convergence (measured 1.04 against 2.00, `dev_notes/issue353/RESULTS.md` 2.3).
#
# What is pinned here: that one setting holds the default; that the construction is
# selectable through `control` and reported back per axis; that its boxes tile
# the axis in the axis's OWN coordinate so a bounded axis cannot leave its
# support; that every decline is a recorded reason and a fall back to the chord
# read rather than an error; and that the shipped engine reproduces the arm
# `dev_notes/issue337/` scored.

PROBS_WC <- c(0.005, 0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975, 0.995)

skip_on_cran_budget()

# --------------------------------------------------------------- the boxes tile

test_that("the boxes tile the axis and carry the shipped masses exactly", {
  skip_on_cran_budget()
  v <- exp(seq(log(0.2), log(1.5), length.out = 6))
  w <- c(0.05, 0.15, 0.3, 0.25, 0.15, 0.10)
  e <- .nl_box_edges(v, "positive")
  expect_length(e, length(v) + 1L)
  expect_false(is.unsorted(e, strictly = TRUE))
  # Edge k is BOTH the upper edge of cell k and the lower edge of cell k+1 --
  # one number serving both, so no mass leaves a cell it was integrated for and
  # none is counted twice. Each coordinate sits inside its own box.
  expect_true(all(v > e[-length(e)] & v < e[-1L]))
  # And the outer two are the cell partition's own.
  expect_equal(c(e[1L], e[length(e)]), .nl_cell_edges(v, "positive"))
  # The CDF at edge k is the cumulative shipped mass through cell k.
  cw <- cumsum(w / sum(w))
  q <- .nl_summary_quantile(v, w, cw[-length(cw)], "positive", "density",
                            "box_uniform")
  expect_equal(unname(q), e[seq_len(length(v) - 1L) + 1L])
  # Interior edges are midpoints in the coordinate the partition lives in.
  expect_equal(e[3L], exp(mean(log(v[2:3]))))
})

test_that("a cumulative mass that rounds past 1 does not break the inversion", {
  skip_on_cran_budget()
  # gcol33/tulpa#388. `cumsum(m)` can round ABOVE 1, and a grid whose trailing
  # cells hold no mass carries that same value from the last positive cell
  # onward, so forcing the LAST entry of the CDF to 1 made the vector DECREASE
  # there and `findInterval()` refused it. The read errored on exactly the node
  # sets this construction exists to handle -- a grid whose outer cells'
  # softmax weight underflowed to zero -- and it was reachable from a default
  # fit as soon as this became the default read.
  #
  # The witness is found by search rather than by hand, because whether a
  # normalized sum rounds up is not something to assert by eye: over uniform
  # log-scale weights with two trailing zeros it takes a few dozen draws.
  set.seed(11L)
  hit <- NULL
  for (i in seq_len(5000L)) {
    k <- sample(3:8, 1L)
    cand <- c(exp(stats::rnorm(k, 0, 3)), 0, 0)
    cs <- cumsum(cand / sum(cand))
    if (cs[length(cs) - 2L] > 1) { hit <- cand; break }
  }
  expect_false(is.null(hit))
  expect_gt(cumsum(hit / sum(hit))[length(hit) - 2L], 1)
  v <- seq_along(hit)
  probs <- c(0.005, 0.025, 0.5, 0.975, 0.995)
  r <- .nl_box_quantile(v, hit, probs, NA_character_)
  expect_true(is.na(r$declined))
  expect_true(all(is.finite(r$q)))
  expect_false(is.unsorted(r$q))
  # It is the READ that was broken, not the masses: the interval still sits
  # inside the partition the coordinates define.
  e <- .nl_cell_edges(sort(unique(v)), NA_character_)
  expect_gte(r$q[1L], e[1L])
  expect_lte(r$q[length(r$q)], e[2L])
  # And the whole dispatcher survives it, which is the path a fit takes.
  expect_true(all(is.finite(.nl_summary_quantile(v, hit, probs, NA_character_,
                                                 "density"))))
})

test_that("the partition comes from the grid, the masses from the weights", {
  skip_on_cran_budget()
  # gcol33/tulpa#337's own rule is "keep the masses, tile the axis", and the two
  # halves come from different places. A cell whose integration weight underflows
  # to exactly 0 still SITS on the axis, and its coordinate is what fixes its
  # neighbour's box edge; filtering the coordinates by weight would shrink that
  # neighbour's box to nothing and collapse the read onto the chord one. Measured
  # on the coarsest rung of `dev_notes/issue357/coarse357b.R` -- 2 cells at 400
  # groups -- the softmax underflows one of the two on 43 of 150 seeds, so this
  # is the difference between measuring the construction and measuring nothing.
  #
  # The chord read filters both together, correctly: its knots ARE the
  # positive-weight coordinates. This one's knots are edges.
  v <- c(0.3309751, 0.9064126)
  w <- c(0, 1)
  h <- 0.5 * diff(log(v))
  bx <- c(exp(mean(log(v))), exp(log(v[2L]) + h))   # cell 2's own box
  p <- c(0.025, 0.5, 0.975)
  q <- .nl_summary_quantile(v, w, p, "positive", "density", "box_uniform")
  expect_equal(q, bx[1L] + p * diff(bx))
  expect_false(isTRUE(all.equal(
    q, .nl_summary_quantile(v, w, p, "positive", "density", "chord"))))
  # It is a box read, not a decline.
  expect_true(is.na(.nl_summary_quantile_read(v, w, p, "positive", "density",
                                              "box_uniform")$declined))

  # An interior empty box is a FLAT segment of the CDF, not a merged one: the
  # quantile is located on the cumulative mass and evaluated inside the box it
  # lands in, so an empty box is stepped over rather than interpolated across.
  v4 <- exp(seq(log(0.2), log(1.5), length.out = 4))
  e4 <- .nl_box_edges(v4, "positive")
  q4 <- .nl_summary_quantile(v4, c(0.5, 0, 0, 0.5), c(0.25, 0.5, 0.75),
                             "positive", "density", "box_uniform")
  expect_gte(q4[1L], e4[1L]); expect_lte(q4[1L], e4[2L])
  expect_equal(q4[2L], e4[4L])
  expect_gte(q4[3L], e4[4L]); expect_lte(q4[3L], e4[5L])
  # Interpolating across the two empty boxes would put q(0.75) below e4[4].
  expect_gt(q4[3L], e4[3L])
})

test_that("a bounded axis's boxes stay inside its support", {
  skip_on_cran_budget()
  # gcol33/tulpa#361 extended `bym2_rho` to six nodes including 0.999, so a
  # `unit` axis with a node that close to its boundary is live. Mirroring or
  # bisecting such an axis in `log` -- which is what a coordinate guessed from
  # `all(v > 0)` does -- puts an edge above 1, and a BYM2 mixing weight is
  # singular there (gcol33/tulpa#369).
  v <- c(0.05, 0.3, 0.6, 0.9, 0.9962, 0.99972)
  w <- c(0.05, 0.10, 0.20, 0.30, 0.20, 0.15)
  e <- .nl_box_edges(v, "unit")
  expect_true(all(e > 0 & e < 1))
  q <- .nl_summary_quantile(v, w, c(0.001, 0.5, 0.999), "unit", "density",
                            "box_uniform")
  expect_true(all(q > 0 & q < 1))
  # Without the domain the log guess is taken and does leave the support, which
  # is the fallback's own behaviour and not a claim about `unit`.
  expect_gt(.nl_box_edges(v, NA_character_)[length(v) + 1L], 1)
  r <- c(-0.9, -0.2, 0.3, 0.8)
  expect_true(all(abs(.nl_box_edges(r, "correlation")) < 1))
})

# ------------------------------------------------------ the per-axis provenance

test_that("the axis read says which construction produced each interval", {
  skip_on_cran_budget()
  v1 <- exp(seq(log(0.2), log(1.5), length.out = 4))
  v2 <- c(0.5, 1.5)
  tg <- as.matrix(expand.grid(sigma = v1, tau = v2))
  lm <- c(0, -0.4, -3, -6, -0.2, -0.6, -3.2, -6.2)
  qc <- .nl_axis_quantiles(tg, lm, within = "chord")
  expect_identical(unname(qc$within), c("chord", "chord"))
  expect_true(all(is.na(qc$within_declined)))
  qb <- .nl_axis_quantiles(tg, lm, within = "box_uniform")
  expect_identical(unname(qb$within), c("box_uniform", "box_uniform"))
  expect_false(isTRUE(all.equal(unname(qb$ci_lo), unname(qc$ci_lo))))
  # An unnamed read takes the engine default, whichever that is.
  q0 <- .nl_axis_quantiles(tg, lm)
  expect_identical(q0, .nl_axis_quantiles(tg, lm,
                                          within = .nl_diag("within_cell")))
  # A single-node axis falls back on its OWN, without taking the fit with it.
  tg1 <- cbind(sigma = tg[, "sigma"], fixed = rep(2, nrow(tg)))
  q1 <- .nl_axis_quantiles(tg1, lm, within = "box_uniform")
  expect_identical(unname(q1$within), c("box_uniform", "chord"))
  expect_identical(unname(q1$within_declined), c(NA_character_, "single_node"))
})

test_that("the per-axis resolution is reported in the axis's own coordinate", {
  skip_on_cran_budget()
  # `h / sd` is the regime variable for the position sensitivity a within-cell
  # reconstruction has, so a fit reports it rather than leaving the reader with
  # no way to tell a resolved axis from an unresolved one.
  v <- exp(seq(log(0.2), log(1.5), length.out = 9))
  u <- log(v)
  # A discretized Gaussian in u, so the answer is known.
  sd_true <- 0.35
  lm <- -0.5 * ((u - mean(u)) / sd_true)^2
  tg <- matrix(v, ncol = 1L, dimnames = list(NULL, "sigma"))
  rs <- .nl_axis_resolution(tg, lm, domains = list("positive"))
  expect_equal(unname(rs$h[["sigma"]]), diff(u)[1L])
  expect_equal(unname(rs$sd[["sigma"]]), sd_true, tolerance = 1e-8)
  expect_equal(unname(rs$h_over_sd[["sigma"]]), diff(u)[1L] / sd_true,
               tolerance = 1e-8)
  # The SD is the one `.nl_laplace_at_mode_sd_axis()` reports in that same
  # coordinate; on a positive axis the linear-scale value is the delta map of
  # it, and the ratio would be a different (and meaningless) number.
  expect_false(isTRUE(all.equal(
    unname(rs$sd[["sigma"]]),
    .nl_laplace_at_mode_sd_axis(v, lm))))
  # An axis whose mode sits on its own boundary has no 3-point parabola, so the
  # ratio is withheld rather than guessed.
  lm_rail <- seq(0, 5, length.out = 9)
  expect_true(is.na(.nl_axis_resolution(tg, lm_rail,
                                        domains = list("positive"))$sd[[1L]]))
})

skip_on_cran()

# --------------------------------------------------------- the default is fixed

test_that("the default construction is log-quadratic, and one setting decides it", {
  skip_on_cran()
  # The default moved on the measurements recorded in `R/settings.R`: chord to
  # box-uniform at 0.0.188 (gcol33/tulpa#357), box-uniform to log-quadratic at
  # gcol33/tulpa#932. What this has to keep saying is that there is ONE default
  # and every place that could name one defers to it -- the readers' `within`
  # arguments default to NULL and resolve through `.nl_within_cell_mode()` --
  # since disagreeing defaults are what make a read depend on which door it
  # came through.
  v <- exp(seq(log(0.2), log(1.5), length.out = 5))
  w <- c(0.30, 0.25, 0.20, 0.15, 0.10)
  dflt <- .nl_diag("within_cell")
  expect_identical(dflt, "log_quadratic")
  expect_identical(.nl_within_cell_mode(NULL), dflt)
  for (fn in list(.nl_posterior_moments, .nl_summary_quantile,
                  .nl_summary_quantile_read, .nl_read_sd, .nl_axis_quantiles,
                  .nl_posterior_moments_multi, .joint_posterior_moments_multi)) {
    expect_null(formals(fn)$within)
  }
  expect_true(dflt %in% .NL_SUPPORT[["density"]]$within)
  for (sup in .NL_SUPPORT_KINDS) {
    for (dm in c("positive", NA_character_)) {
      # An unnamed call takes the default, and a kind that does not admit it
      # falls back to the chord read -- which is the same number either way.
      want <- if (dflt %in% .NL_SUPPORT[[sup]]$within) dflt else "chord"
      expect_identical(.nl_summary_quantile(v, w, PROBS_WC, dm, sup),
                       .nl_summary_quantile(v, w, PROBS_WC, dm, sup, want))
    }
  }
  # `.nl_summary_quantile()` is `.nl_summary_quantile_read()`'s `$q` -- one
  # dispatcher, so the vector caller and the provenance caller cannot drift.
  r <- .nl_summary_quantile_read(v, w, PROBS_WC, "positive", "density")
  expect_identical(r$q, .nl_summary_quantile(v, w, PROBS_WC, "positive",
                                             "density"))
  # A quantity read without its grid rows has none to reconstruct, so the
  # log-quadratic read declines to the boxes and says why.
  expect_identical(r$within, "box_uniform")
  expect_identical(r$declined, "not_a_grid_axis")
  # The paired assertion the flip needs: the default moved the read on the ONE
  # kind that admits box-uniform and on nothing else, and `"chord"` is still
  # reachable and still different there.
  for (sup in .NL_SUPPORT_KINDS) {
    for (dm in c("positive", NA_character_)) {
      d <- .nl_summary_quantile(v, w, PROBS_WC, dm, sup)
      c0 <- .nl_summary_quantile(v, w, PROBS_WC, dm, sup, "chord")
      if (identical(sup, "density")) expect_false(isTRUE(all.equal(d, c0)))
      else expect_identical(d, c0)
    }
  }
})

test_that("`.nl_cell_edges()` is unchanged by the partition extraction", {
  skip_on_cran()
  # `.nl_cell_partition()` is what the box read needs: the coordinate the outer
  # edges were mirrored in is not recoverable from the edges, because the guard
  # falls back when a mapped edge leaves the support. `.nl_cell_edges()` is now
  # its `$edges` and must be the same function it was.
  cases <- list(exp(seq(log(0.2), log(1.5), length.out = 5)),
                c(-2, 0, 2), c(0.05, 0.3, 0.6, 0.9, 0.9962, 0.99972),
                c(-0.9, -0.2, 0.3, 0.8), c(0.5, 1.5, 2.5), 3)
  for (v in cases) {
    for (dm in c(NA_character_, "positive", "unit", "correlation",
                 "unbounded")) {
      p <- .nl_cell_partition(sort(v), dm)
      expect_identical(.nl_cell_edges(sort(v), dm), p$edges)
      # The reported coordinate really is the one the edges came out of.
      if (length(v) > 1L) {
        expect_equal(p$tr$to(p$edges),
                     local({
                       u <- p$tr$to(sort(v)); n <- length(u)
                       c(u[1L] - 0.5 * (u[2L] - u[1L]),
                         u[n] + 0.5 * (u[n] - u[n - 1L]))
                     }))
      }
    }
  }
})

test_that("the interior bisector is byte-identical to the sum form it replaced", {
  skip_on_cran()
  # gcol33/tulpa#378 changed HOW the interior midpoint is formed, and the whole
  # claim is that it changes no number any axis actually carries. `a / 2 + b / 2`
  # and `(a + b) / 2` are not the same expression -- the first rounds twice --
  # but halving a NORMAL double is exact, so both carry a single rounding of the
  # same real and agree bit for bit. Asserted with `identical()` against a
  # reference built here from the sum form, over randomized node sets on every
  # domain, not on a handful of values.
  set.seed(378L)
  ref_edges <- function(v, dm) {
    pt <- .nl_cell_partition(v, dm)
    n <- length(v)
    u <- pt$tr$to(v)
    ue <- pt$tr$to(pt$edges)
    e <- pt$tr$from(c(ue[1L], (u[-1L] + u[-n]) / 2, ue[2L]))
    if (!all(is.finite(e)) || is.unsorted(e, strictly = TRUE)) return(NULL)
    e
  }
  gen <- list(
    positive    = function(k) exp(runif(k, -300, 300)),
    unit        = function(k) stats::plogis(runif(k, -300, 300)),
    correlation = function(k) tanh(runif(k, -15, 15)),
    unbounded   = function(k) ifelse(runif(k) < 0.5, -1, 1) * 10^runif(k, -8, 8))
  n_case <- 0L
  for (i in seq_len(400L)) {
    for (dm in names(gen)) {
      v <- sort(unique(gen[[dm]](sample(2:9, 1))))
      if (length(v) < 2L) next
      n_case <- n_case + 1L
      # The edges are the whole difference: the read is a pure function of them
      # and the masses, so identical edges are identical intervals.
      expect_identical(.nl_box_edges(v, dm), ref_edges(v, dm))
    }
  }
  expect_gt(n_case, 1500L)
  # The one place the two forms part is the one #378 is about: a sum that leaves
  # the double range while the midpoint does not.
  expect_true(is.infinite((1e300 + .Machine$double.xmax) / 2))
  expect_true(is.finite(1e300 / 2 + .Machine$double.xmax / 2))
})

# ------------------------------------------------------ declines, never errors

test_that("every decline records a reason and falls back to the chord read", {
  skip_on_cran()
  v <- exp(seq(log(0.2), log(1.5), length.out = 5))
  w <- c(0.30, 0.25, 0.20, 0.15, 0.10)
  # A support whose node set is not a tiling partition declines by NAME, so a
  # caller can tell "this kind does not admit it" from "this node set failed".
  for (sup in c("moment_rule", "mixed", "sample")) {
    r <- .nl_summary_quantile_read(v, w, PROBS_WC, "positive", sup,
                                   "box_uniform")
    expect_identical(r$within, "chord")
    expect_identical(r$declined, paste0("support_", sup))
    expect_identical(r$q, .nl_summary_quantile(v, w, PROBS_WC, "positive", sup))
  }
  # One node has no box; the chord read's answer (the value itself) stands.
  r1 <- .nl_summary_quantile_read(3, 1, PROBS_WC, "positive", "density",
                                  "box_uniform")
  expect_identical(r1$declined, "single_node")
  expect_identical(r1$q, rep(3, length(PROBS_WC)))
  # No usable node at all is NA, as it is under either construction.
  r0 <- .nl_summary_quantile_read(c(1, 2), c(0, 0), PROBS_WC, "positive",
                                  "density", "box_uniform")
  expect_identical(r0$declined, "no_usable_node")
  expect_true(all(is.na(r0$q)))
  # `slice_masses()` in `dev_notes/issue337/recon.R` STOPS when the boxes do not
  # tile, because mass would leave the cell it was integrated for. A reported
  # interval cannot take that behaviour, so the engine declines instead
  # (gcol33/tulpa#293).
  #
  # The interior midpoint is no longer one of the ways that happens
  # (gcol33/tulpa#378). It is formed as `a / 2 + b / 2`, so two coordinates near
  # the top of the double range have a representable bisector and the partition
  # they define is built; `(a + b) / 2` reached `Inf` on the sum and declined a
  # partition every one of whose edges is a double.
  vd <- c(1, 1e300, .Machine$double.xmax)
  expect_true(is.infinite((1e300 + .Machine$double.xmax) / 2))
  eb <- .nl_box_edges(vd, "unbounded")
  expect_length(eb, length(vd) + 1L)
  expect_false(is.unsorted(eb, strictly = TRUE))
  expect_true(all(vd >= eb[-length(eb)] & vd <= eb[-1L]))
  expect_identical(eb[3L], 1e300 / 2 + .Machine$double.xmax / 2)
  expect_identical(
    .nl_summary_quantile_read(vd, c(1, 1, 1), PROBS_WC, "unbounded", "density",
                              "box_uniform")$within, "box_uniform")
  # UNDECLARED, that node set no longer declines either (gcol33/tulpa#379). Its
  # outer mirror -- `xmax + (xmax - 1e300) / 2` -- is past the double range in
  # any form, and no midpoint arithmetic reaches it, which is what made it the
  # remaining fixture after gcol33/tulpa#378. What it needed was the guard the
  # DECLARED branch already had: an unusable mirror falls back to the extreme
  # coordinates. So the partition is finite, the boxes tile, and the read the
  # caller asked for runs.
  expect_true(is.infinite(vd[3L] + 0.5 * (vd[3L] - vd[2L])))
  expect_identical(.nl_cell_edges(vd, NA_character_), c(vd[1L], vd[3L]))
  expect_false(is.null(.nl_box_edges(vd, NA_character_)))
  rd <- .nl_summary_quantile_read(vd, c(1, 1, 1), PROBS_WC,
                                  NA_character_, "density", "box_uniform")
  expect_true(is.na(rd$declined))
  expect_identical(rd$within, "box_uniform")
  expect_identical(rd$edge_declined, "mirrored_edge_not_representable")
  expect_true(all(is.finite(rd$q)))
  # The decline that REMAINS is a node set the coordinate map cannot separate. A
  # `unit` axis reaching the subnormal floor has `plogis(qlogis(1e-320))`
  # underflow to 0, so two adjacent edges come back as the same double. That is
  # a property of the node set against the coordinate, which no arithmetic and
  # no edge guard changes.
  vu <- c(1e-320, 1e-310, 0.5)
  expect_null(.nl_box_edges(vu, "unit"))
  ru <- .nl_summary_quantile_read(vu, c(1, 1, 1), PROBS_WC, "unit", "density",
                                  "box_uniform")
  expect_identical(ru$declined, "boxes_do_not_tile")
  expect_identical(ru$within, "chord")
  expect_identical(ru$q, .nl_summary_quantile(vu, c(1, 1, 1), PROBS_WC, "unit",
                                              "density"))
  # A spacing below the coordinate's own resolution collapses a box to zero
  # width and is refused at the edge builder, but it does not reach the read:
  # `.nl_axis_atoms()` has already merged coordinates that are equal as doubles,
  # so the partition the read sees never has a zero-width interior box.
  expect_null(.nl_box_edges(c(1, 1, 2), "positive"))
  expect_false(is.null(.nl_box_edges(
    .nl_axis_atoms(c(1, 1, 2), c(1, 1, 1))$v, "positive")))
  # Declared `positive`, the SAME node set does not reach that decline at all
  # any more (gcol33/tulpa#377). Its mirrored upper edge overflows, and the
  # partition then holds the declaration instead of falling through to the value
  # guess: the outer edges are the extreme coordinates, which the containment
  # test has already put inside the support, so both reads have a finite
  # in-support partition to work with.
  expect_identical(.nl_cell_edges(vd, "positive"), c(vd[1L], vd[3L]))
  expect_false(is.null(.nl_box_edges(vd, "positive")))
  expect_identical(
    .nl_summary_quantile_read(vd, c(1, 1, 1), PROBS_WC, "positive", "density",
                              "box_uniform")$within, "box_uniform")
  # An unknown construction is refused rather than translated.
  expect_error(.nl_summary_quantile(v, w, 0.5, "positive", "density", "boxes"))
  expect_error(.nl_within_cell_mode("uniform"))
})

# ----------------------------------------------------------- through a real fit

test_that("a fit selects the construction, reports it, and stays byte-identical by default", {
  skip_on_cran()
  skip_if_fast()
  set.seed(357L)
  G <- 40L; npg <- 30L; N <- G * npg
  g <- sample.int(G, N, replace = TRUE)
  x <- rnorm(N); X <- cbind(1, x)
  y <- as.numeric(X %*% c(-0.2, 0.7)) + rnorm(G, 0, 0.7)[g] +
       rnorm(N, 0, sqrt(0.5))
  sg <- exp(seq(log(0.2), log(1.5), length.out = 9L))
  # The read is rebuilt below from the nine nodes written down, so no node is
  # added to the axis.
  ctrl <- list(max_iter = 100L, tol = 1e-8, n_threads = 1L, diagnose_k = FALSE,
               diagnose_skew = FALSE, integration = "grid", local_ccd = NULL,
               auto_recenter = FALSE, axis_refine = "none")
  fit <- function(wc) suppressWarnings(tulpa_nested_laplace_joint(
    responses = list(a = list(y = y, n_trials = rep(1L, N), X = X,
                              family = "gaussian", phi = 0.5)),
    prior = list(list(type = "iid", obs_idx = list(g), n_units = G,
                      sigma_grid = sg)),
    control = if (is.null(wc)) ctrl else c(ctrl, list(within_cell = wc))))

  f0 <- fit(NULL)
  fc <- fit("chord")
  fb <- fit("box_uniform")
  fl <- fit("log_quadratic")

  # Asking for the default explicitly changes nothing. The default is
  # `log_quadratic` (gcol33/tulpa#932), so the arm `f0` matches is `fl`; `fb`
  # and `fc` are the alternatives that must still be reachable.
  expect_identical(f0$theta_median, fl$theta_median)
  expect_identical(f0$theta_ci_lo, fl$theta_ci_lo)
  expect_identical(f0$theta_ci_hi, fl$theta_ci_hi)
  expect_false(isTRUE(all.equal(unname(fb$theta_ci_lo),
                                unname(fc$theta_ci_lo))))
  # The construction is reported, and every path stamps the node-set kind now,
  # not only the multi-block driver (gcol33/tulpa#357).
  expect_identical(f0$within_cell_requested, "log_quadratic")
  expect_identical(fc$within_cell_requested, "chord")
  expect_identical(f0$theta_within_cell, fl$theta_within_cell)
  expect_identical(unname(fb$theta_within_cell), "box_uniform")
  expect_identical(unname(fc$theta_within_cell), "chord")
  expect_identical(f0$theta_interval_read, "density")
  expect_equal(f0$theta_interval_design_mass, 0)
  # The moments are the SAME under either construction: they are sums over the
  # cells, and this changes only where inside a cell the mass sits.
  expect_identical(fc$theta_mean, fb$theta_mean)
  expect_identical(fc$theta_sd, fb$theta_sd)
  # And the interval is not: the box read is the narrower of the two.
  expect_lt(unname(fb$theta_ci_hi - fb$theta_ci_lo),
            unname(fc$theta_ci_hi - fc$theta_ci_lo))
  # The fit knows how coarse its own grid is.
  expect_true(is.finite(f0$outer_grid_h_over_sd[[1L]]))
  expect_equal(unname(f0$outer_grid_cell_width[[1L]]), diff(log(sg))[1L])

  # The shipped read IS the arm `dev_notes/issue337/recon.R` scored: cell masses
  # over the cells' own boxes, the CDF read by the chord between edges. Rebuilt
  # here from the fit's own grid, independently of the engine's code path.
  w <- exp(fb$log_marginal - max(fb$log_marginal)); w <- w / sum(w)
  v <- as.numeric(fb$theta_grid); o <- order(v); v <- v[o]; w <- w[o]
  u <- log(v); n <- length(u)
  mid <- (u[-1L] + u[-n]) / 2
  e <- exp(c(u[1L] - (mid[1L] - u[1L]), mid, u[n] + (u[n] - mid[n - 1L])))
  Fc <- c(0, cumsum(w))
  ref <- stats::approx(Fc, e, xout = c(0.025, 0.5, 0.975), ties = "ordered",
                       rule = 2)$y
  expect_equal(unname(c(fb$theta_ci_lo, fb$theta_median, fb$theta_ci_hi)),
               ref, tolerance = 1e-12)

  # `diagnostics()` reads the construction back and says what it means. The
  # note fires on a read that is NOT the engine default, so it is the `chord`
  # fit that carries one and the default fit that does not (gcol33/tulpa#357)
  # -- naming the default as a literal in the note is what made it say "rather
  # than the default 'chord'" on a fit read with the default.
  d <- diagnostics(f0)
  expect_identical(attr(d, "within_cell_requested"), "log_quadratic")
  expect_identical(attr(d, "interval_read"), "density")
  expect_false(any(grepl("within-cell", attr(d, "interval_read_note"))))
  dc <- diagnostics(fc)
  expect_identical(attr(dc, "within_cell_requested"), "chord")
  expect_true(any(grepl("'chord' within-cell", attr(dc, "interval_read_note"),
                        fixed = TRUE)))
  expect_true(is.finite(attr(d, "grid_h_over_sd_max")))
  # This fixture is in the unresolved regime -- as every configuration of the
  # 34-row census is -- so the resolution note fires and names the axis.
  expect_gt(attr(d, "grid_h_over_sd_max"), .nl_diag("grid_resolved"))
  expect_true(grepl(attr(d, "grid_coarsest_axis"),
                    attr(d, "grid_resolution_note"), fixed = TRUE))
})

test_that("an unknown within_cell is refused at the front door", {
  skip_on_cran()
  skip_if_fast()
  set.seed(358L)
  N <- 200L
  X <- cbind(1, rnorm(N))
  g <- sample.int(10L, N, replace = TRUE)
  y <- as.numeric(X %*% c(0, 1)) + rnorm(N)
  expect_error(suppressWarnings(tulpa_nested_laplace_joint(
    responses = list(a = list(y = y, n_trials = rep(1L, N), X = X,
                              family = "gaussian", phi = 1)),
    prior = list(list(type = "iid", obs_idx = list(g), n_units = 10L,
                      sigma_grid = c(0.5, 1, 2))),
    control = list(diagnose_k = FALSE, diagnose_skew = FALSE,
                   within_cell = "quadratic"))))
})