R/nested_laplace_joint.R

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

#' Joint multi-likelihood nested Laplace approximation
#'
#' @description
#' Outer-grid nested Laplace driver for *joint* models -- multiple response
#' arms sharing one latent prior block, parameterized as a per-arm field
#' amplitude (sigma) on a unit-precision latent.
#'
#' Supported priors:
#'  * `"bym2"`       -- outer grid over `(sigma, rho [, alpha])`.
#'                     Latent: `phi (n_s) | theta (n_s)` with unit-precision
#'                     ICAR + iid components.
#'  * `"icar"`       -- outer grid over `(sigma [, alpha])`.
#'                     Latent: `phi (n_s)` with unit-precision ICAR.
#'  * `"car_proper"` -- outer grid over `(sigma, rho_car [, alpha])`.
#'                     Latent: `phi (n_s)` with `Q = D - rho_car * W`.
#'
#' Other backends (NNGP, HSGP, RW1/2, AR1) follow the same interface and
#' land under Phase 3.
#'
#' @section (sigma, alpha) parameterization:
#' Each arm's linear predictor reads
#' \deqn{\eta_{arm} = X_{arm} \beta_{arm} + \sigma_{arm} \cdot z_s,}
#' where \eqn{z_s} is a unit-precision latent (ICAR(tau=1) for ICAR/BYM2,
#' or the BYM2 mix; CAR_proper uses the structure of
#' \eqn{D - \rho_{car} W}). All arms share the donor amplitude
#' \eqn{\sigma} from the outer-grid `sigma_grid` axis. The copy arm's
#' amplitude is \eqn{\sigma_{arm} = \alpha \cdot \sigma}, where \eqn{\alpha}
#' is a direct outer-grid axis taken from `copy$alpha_grid`. The Cartesian
#' product is over `(sigma, [rho/rho_car,] alpha[, phi])`. Direct \eqn{\alpha}
#' as an outer-grid axis (rather than a post-hoc ratio
#' \eqn{\alpha = \sigma_{pos} / \sigma_{occ}}) avoids plug-in bias on the
#' weakly-identified ratio at small `n_pos` and lets a regularizing
#' hyperprior land on \eqn{\alpha} directly.
#'
#' @param responses A named list of arm specs (length >= 1). Each arm:
#'   * `y`           -- numeric `[N_arm]` response.
#'   * `n_trials`    -- integer `[N_arm]` (use `rep(1L, N_arm)` for non-binomial).
#'   * `X`           -- numeric matrix `[N_arm x p_arm]` fixed-effects design.
#'   * `spatial_idx` -- integer `[N_arm]`, 1-based map obs -> spatial unit.
#'   * `re_idx`      -- optional numeric `[N_arm]` 1-based RE group index;
#'                     defaults to `rep(0, N_arm)` (no RE).
#'   * `n_re_groups` -- optional integer (default `0L`).
#'   * `sigma_re`    -- optional numeric (default `1`); ignored when
#'                     `n_re_groups == 0`.
#'   * `family`      -- one of `"binomial"`, `"gaussian"`, `"poisson"`,
#'                     `"neg_binomial_2"`, `"beta"`, `"lognormal"`,
#'                     `"gamma"`, `"inverse_gaussian"`. For `"lognormal"`,
#'                     `y` is on the natural scale and the linear
#'                     predictor `eta = E[log y]` (identity link on the
#'                     log scale); the `-log(y)` Jacobian is included in
#'                     the kernel's `log_lik`.
#'   * `phi`         -- numeric dispersion (gaussian/lognormal residual
#'                     SD, negbin size, beta precision); default `1`.
#'   * `field_coef`  -- optional per-arm field coefficient controlling
#'                     this arm's multiplier on the shared latent field's
#'                     amplitude. One of:
#'                     * numeric scalar (default `1`) -- constant
#'                       multiplier. `0` means the arm carries NO field
#'                       at all (the per-row scatter `eta += field_coef *
#'                       sigma * z` is skipped for that arm).
#'                     * character of length 1 -- names an outer-grid
#'                       hyperparam axis (currently `"alpha"`); the
#'                       coefficient varies across the grid.
#'                     * `list(name = , grid = )` -- embedded axis
#'                       declaration, equivalent to declaring the axis
#'                       and naming it on this arm.
#'                     At most one arm may declare a hyperparam-driven
#'                     axis (the cover hurdle and
#'                     occu_cover both need only one). Shared axes
#' across multiple arms are deferred. A single-block
#' copy coefficient is declared here, on the arm --
#' not through a separate `copy` argument.
#'
#' @param prior A list describing the shared latent prior block. Required
#'   field `type`. Backend-specific fields:
#'   * **bym2**: `n_spatial_units`, `adj_row_ptr`, `adj_col_idx`,
#'     `n_neighbors`, `scale_factor` (default `1`); optional `sigma_grid`
#'     (donor-arm field amplitude, default 5 log-spaced values in
#'     `[0.1, 3]`), `rho_grid` (default
#'     `c(0.2, 0.5, 0.8, 0.95, 0.99, 0.999)`).
#'   * **icar**: `n_spatial_units`, `adj_row_ptr`, `adj_col_idx`,
#'     `n_neighbors`; optional `sigma_grid` (default 5 log-spaced values
#'     in `[0.1, 3]`).
#'   * **car_proper**: same as icar plus `rho_car_grid`
#'     (default `c(0.5, 0.8, 0.95, 0.99)`).
#'
#'   `sigma_grid`'s default is a starting axis, not a hard ceiling: when the
#'   fitted field-SD posterior mode rails the top node (the axis's own
#'   marginal is maximal there, or the whole grid collapsed onto it:
#'   `pareto_k_regime = "collapsed_edge"`, see below), the driver finds the
#'   outer posterior mode -- a damped Newton search over every outer axis,
#'   started at the grid's heaviest cell and read through full inner solves,
#'   the same mode-find the CCD integrator places its design with -- lays the
#'   axis at that mode +/- 2.5 of the SDs its curvature there gives, and refits
#'   (up to two attempts, the second adding a light
#'   default PC(U=3, alpha=0.01) prior on sigma unless `prior_sigma` was pinned
#'   -- see there), so a sparse or strongly-identified species is not silently
#'   truncated at 3.0. This engages whether or not `control$diagnose_k`
#'   computed the full outer Pareto-k diagnostic, and only when an axis
#'   actually railed. A field SD still on a boundary of its grid after the
#'   attempts is named in a warning, with the reason it stayed, in
#'   `outer_grid_axis_declined` (gcol33/tulpa#919). The mode the placement
#'   was laid around is on the fit as `outer_mode_u` / `outer_mode_cov_u` (in
#'   each axis's `outer_mode_axis_tags` coordinate), with `outer_mode_status`,
#'   `outer_mode_rounds` and `outer_mode_evals`. A `sigma_grid` the caller
#'   PINNED always wins: auto-recenter
#'   engages when the field is left `NULL`, when it is marked with
#'   [auto_grid()] (how a wrapper package declares an axis it defaulted rather
#'   than one the user chose), or when its nodes are exactly
#'   the engine's own default axis. Declines gracefully (keeps the fixed-grid
#'   fit) when another axis in the same grid has unguessable support
#'   (car_proper's `rho_car`); whichever way it declines,
#'   `outer_grid_recenter_declined` says which (see below).
#'
#'   On the list-of-blocks path an `icar`, `bym2`, `car_proper` or `nngp`
#'   block may carry `projector` in place of `spatial_idx`: a list of `n_arms`
#'   matrices (dense or a \pkg{Matrix} sparse matrix), arm `k`'s of size
#'   `N_k x n_spatial_units` (`N_k x n_spatial` for `nngp`), or one matrix
#'   shared by every arm; a `NULL` entry is an arm the field does not reach. Arm `k`'s predictor then receives
#'   `A_k z` with `z` on the block's own prior, so restricted spatial
#'   regression is `A_k = P_k S_k`, with `P_k = I - X_k (X_k' X_k)^-1 X_k'` and
#'   `S_k` the observation-to-unit incidence. An areal field's level is
#'   placed from `A_k 1`: where it is zero on every arm the level never
#'   reaches the predictor and is removed with no fold; where it is the
#'   intercept column it folds into the intercept. The `nngp` and
#'   `car_proper` priors are proper and carry their level themselves. Restricting the field changes the estimand
#'   to the marginal association (see [spatial_rsr()]).
#'
#'   An empty `prior = list()` is a list of zero blocks: the model carries no
#'   latent field, each arm's predictor is its fixed effects, and the outer
#'   grid is spanned by the per-arm dispersion axes of `phi_grid` alone (one
#'   cell when there are none).
#'
#' @param copy Multi-block copy specification (multi-block `prior` only). For a
#' single-block fit there is no `copy` argument: declare the copy coefficient
#' on the arm via `responses[[X]]$field_coef = list(name = "alpha", grid = G)`.
#' On the multi-block path `copy` is an unnamed list of specs --
#' `list(list(arm, block, alpha_grid),...)` -- coupling N distinct shared
#' latent fields, each onto its own arm with its own \eqn{\alpha} axis,
#' integrated over the product outer grid. Each spec must name a distinct
#' block.
#'
#' The \eqn{\alpha} axis takes EITHER `alpha_grid`, which states its nodes, OR
#' `alpha_n` (`field_coef$n` on the single-block path), which re-reads the
#' engine's own axis at a higher RESOLUTION -- more nodes between the same
#' bounds, keeping the atom at 0 that gives the "no copy" base model posterior
#' mass. They are different requests and supplying both is an error. Use
#' `alpha_n` to integrate the same axis more accurately: it is the only outer
#' axis a copy fit could not otherwise raise, since stating nodes for it also
#' restates the prior structure it carries (gcol33/tulpa#633).
#' The copy block may be any of `icar` / `bym2` / `car_proper` / `rw1`
#' / `rw2` / `ar1` / `iid`; blocks with their own per-arm scaling (`lf`,
#'   `hsgp_mo`) or a precomputed precision (`tgmrf`) do not take a copy. A
#'   copy block's own `sigma_grid` (the donor field amplitude, same default
#'   ceiling as the single-block `prior$sigma_grid` above) auto-recenters on
#'   `collapsed_edge` the same way, one block per attempt.
#'
#' @param phi_grid Optional list specifying per-arm dispersion axes on the
#'   outer grid. Accepts either a named list (keys = arm names) or a
#'   positional list of length `n_arms`. Each entry is one of:
#'   * `NULL` or scalar -- no axis for that arm; the kernel uses the
#'     parse-time scalar `responses[[k]]$phi`.
#'   * numeric vector of length > 1 -- adds a new outer-grid axis
#'     `phi_<arm>` taking those values; the kernel rewrites
#'     `arms[k].phi` at each grid point before the inner Newton solve.
#'
#'   Family-specific interpretation of `arm$phi` (the parse-time scalar
#'   and the grid values):
#'   * `gaussian` -- residual VARIANCE, the one engine convention (the SD is
#'     `sqrt(phi)`). Use `phi_grid` to estimate it as a hyperparameter
#'     instead of pinning it pre-fit.
#'   * `lognormal` -- residual variance on the log scale; identical kernel
#'     parameterization as `gaussian` plus the `-log(y)` Jacobian.
#'   * `neg_binomial_2` -- dispersion (variance is `mu + mu^2/phi`).
#'   * `beta` -- precision (variance is `mu(1-mu)/(1+phi)`).
#'   * `gamma`, `inverse_gaussian` -- shape / dispersion.
#'   * `binomial`, `poisson` -- ignored.
#'
#'   A dispersion axis is PLACED like a prior block's scale axis when the caller
#'   declares it a default with [auto_grid()] -- `phi_grid = list(pos =
#'   auto_grid(nodes))`. The engine has no default dispersion axis of its own,
#'   so an unmarked vector is read as a pin and integrated exactly as written;
#'   the fit records that per axis in `outer_grid_axis_declined`. Marked, the
#'   axis is re-laid on `mode +/- span * sd` from the same mode/Hessian stencil
#'   the field-SD axes are placed from, which on a span wide enough to hold no
#'   node near the posterior is the difference between an estimate and an
#'   endpoint.
#'
#'   Each `phi_<arm>` axis is appended to the Cartesian product and
#'   varies slowest (within-spatial warm starts hold). The axis appears
#'   as a regular hyperparameter in `theta_grid`, `theta_mean`, and
#'   `theta_sd`, and participates in adaptive-grid refinement.
#'
#' @param prior_sigma,prior_alpha Optional regularizing hyperpriors on the
#'   donor field amplitude \eqn{\sigma} and on the copy coefficient
#'   \eqn{\alpha}. Each is `NULL` (default: the axis carries what `hyperprior`
#'   gives it) or a list of the form `list(family, params)`:
#'   * `list("pc.prec", c(U, alpha))` -- Penalized Complexity prior,
#'     calibrated by `P(theta > U) = alpha`. Closed-form density
#'     `lambda * exp(-lambda * theta)` with `lambda = -log(alpha)/U`.
#'     Drop-in for the weakly-identified small-`n_pos` regime
#'. Pick `U` at the upper end of plausible values
#'     so the prior shrinks the tail without biasing the modal cell when
#'     the data identifies it: default-friendly choice on \eqn{\sigma} is
#'     `c(U = 1.0, alpha = 0.01)` (donor amplitude); on the dimensionless
#'     copy coefficient \eqn{\alpha} the recommended choice is
#'     `c(U = 4.0, alpha = 0.01)`. What counts as strong is set by the
#'     measure the axis already carries, and on the copy scale that is not
#'     flat: with no `prior_alpha` it takes the engine's own exponential
#'     slab, whose rate puts 5 % of the prior above the largest declared
#'     node (`control$copy_slab`). A `prior_alpha` REPLACES that slab, so it
#'     shrinks harder than the default only once its own
#'     `lambda = -log(alpha)/U` beats the slab's -- on a copy grid reaching
#'     `alpha = 2`, `lambda = 1.50`, the strength of `U = 3.07`. Too small a
#'     `U` over-shrinks the copy coefficient past the modal cell and, through
#'     the `alpha * sigma` copy axis, inflates the coupled donor amplitude
#'     `sigma`: on a fixture with truth \eqn{\alpha = 1} and `sigma = 0.6`,
#'     `c(U = 2.0, alpha = 0.01)` pulls the \eqn{\alpha} geometric mean to
#'     0.93 and lifts `sigma` to 0.68, and `c(U = 0.25, alpha = 0.01)` to
#'     0.51 and 0.94.
#'   * `list("half_normal", scale)` -- half-normal with scale `scale > 0`.
#'     Sharper tail decay than PC; use when stronger regularization is
#'     desired and the truth is well inside the prior.
#'   The contribution is added to `log_marginal` cell-by-cell at the
#'   kernel-call boundary, so refinement passes (adaptive grid,
#'   var-of-means consistency) see the regularized posterior. When the
#'   data identifies the parameter (e.g. `n_pos >= ~200`) the prior is
#'   essentially harmless -- the lever is tail-shrinkage at small
#'   `n_pos`. `prior_alpha` only applies when `copy` is active;
#'   `prior_sigma` applies on any `sigma`-named axis.
#'
#'   WHICH AXES ONE SPEC REACHES. A list-of-blocks `prior` can carry the same
#'   axis on several blocks at once -- two copied fields carry `b1.alpha` and
#'   `b2.alpha`, and each copy block carries its own donor `b<k>.sigma`. One
#'   spec applies to EVERY block carrying the axis it names, each block taking
#'   its own independent copy of the density (its atom / continuum split is
#'   normalised against its own axis, so nothing is shared between blocks).
#'   To regularize the blocks differently, or only some of them, key the prior
#'   by block the way a multi-block `copy` spec does:
#'   ```
#'   prior_alpha = list(
#'     list(block = 1, prior = list("pc.prec", c(4.0, 0.01))),
#'     list(block = 2, prior = list("half_normal", 2.0))
#'   )
#'   ```
#'   A block the list does not name carries no prior on that axis. An entry
#'   naming a block that carries no such axis is an error, and the message
#'   names the blocks that do. The per-block shape is a multi-block one: on a
#'   single-block `prior`, where the axis is unique, pass the spec itself.
#'   A fit reports what it applied in `hyperprior_axes`: one record per role
#'   with the scope the caller stated it at (`"every_block"` / `"per_block"`),
#'   the block indices reached, and the grid columns those are.
#'
#'   `prior_sigma` also interacts with the auto-recenter above: its second
#'   attempt engages the engine's own weakly-informative PC(U = 3,
#'   alpha = 0.01) prior, and a `prior_sigma` the caller PINNED suppresses
#'   that (the caller's prior stands). Pinning is decided by provenance, not
#'   presence: a spec marked with [auto_grid()], or one
#'   equal by value to the engine's own default, is a default and does not
#'   suppress the escalation. When it does, the fit carries
#'   `outer_grid_prior_declined = "prior_pinned"`, so a second attempt that
#'   changed only the grid geometry is legible rather than looking like the
#'   full escalation.
#'
#' @param prior_phi Optional regularizing hyperprior on the per-arm
#'   dispersion axes declared through `phi_grid` (e.g. a Beta precision on a
#'   cover arm, a negbin dispersion, a Gaussian residual SD). Same families
#'   as `prior_sigma` -- `NULL` (default: the axis carries what `hyperprior`
#'   gives it), `list("pc.prec", c(U, alpha))`, or `list("half_normal", scale)`.
#'   A single spec re-weights every `phi_<arm>` axis on the grid, the way
#'   `prior_sigma` re-weights any sigma-named axis; with no `phi_grid` it is
#'   a no-op. The PC scale is the dispersion's own units (a precision for
#'   `beta`, a size for `neg_binomial_2`), so pick `U` at the upper end of
#'   plausible values.
#'
#' @template hyperprior
#'
#' @param cell_coupling Character scalar naming a per-cell coupled likelihood
#'   registered against tulpa's process-global registry (see
#'   `src/cell_coupling_registry.h`). Defaults to `"separable"`, the arm-
#'   separable per-obs sum every existing joint fit uses. Consumer packages
#'   (e.g. tulpaObs) compile a `tulpa::CellCouplingSpec` subclass in their
#'   own `src/` and register it from `R_init_<pkg>` via the
#'   `tulpa_register_cell_coupling` C callable; the R driver validates the
#'   name against the registry and the inner Newton routes the per-cell
#'   contribution through `evaluate_cell()` when the spec couples at least
#' one arm.
#'
#' @param control Optional list of perf/numerical tuning knobs (statistical
#'   arguments stay top-level), following the `control` convention of
#'   [tulpa()]. Recognised elements (defaults in parentheses):
#'   * `max_iter` (`50L`), `tol` (`1e-6`) -- inner Newton iteration budget and
#'     tolerance.
#'   * `n_threads` (`1L`) -- inner-loop OpenMP threads (per-observation
#'     scatter, compute_eta, log-likelihood reduction). For typical joint
#'     workloads (`N` in the hundreds to a few thousand) inner parallelism is
#'     overhead-dominated; prefer `n_threads_outer`, which stacks better on
#'     many-core hardware.
#'     Capped at the physical performance-core count by default (see
#'     `n_threads_scatter`), as the inner per-observation loops oversubscribe a
#'     hybrid CPU's efficiency cores past that point.
#'     A fit is reproducible bit for bit at a GIVEN `n_threads`: the
#'     per-observation sum cuts its range into that many contiguous chunks and
#'     adds the chunk sums in chunk order, so nothing about the answer is left
#'     to the OpenMP runtime. Chunking imposes its own association, so two
#'     different `n_threads` agree only to floating-point tolerance (measured
#'     at 6e-14 on `log_marginal` over four families).
#'   * `n_threads_scatter` (performance-core count) -- cap on the inner
#'     per-observation threads. Overrides the default performance-core cap on
#'     `n_threads`; raise it to use all logical cores or lower it to leave
#'     headroom. No effect where the core topology cannot be resolved (off
#'     Windows), where `n_threads` is used as requested.
#'   * `n_threads_outer` (`1L`) -- outer-grid OpenMP threads. When `> 1`, a
#'     pilot Laplace at the centre cell warm-starts the remaining cells, each
#'     dispatched across `n_threads_outer` threads with its own CHOLMOD solver
#'     and NewtonScratch (inner OpenMP auto-disabled). `1L` is serial, chained
#'     warm-starts -- bitwise identical to the pre-speedup driver. Recommended
#'     on multi-core workstations: `parallel::detectCores() - 1L`.
#'     A solve that runs alone holds every thread the fit was granted, the
#'     larger of `n_threads` and `n_threads_outer` under the `n_threads_scatter`
#'     cap, on its inner loop: the pilot, every cell of a refinement or
#'     var-of-means consistency round, and the skew probe. A fit is therefore
#'     reproducible bit for bit at a given pair of `n_threads` and
#'     `n_threads_outer`.
#'     It is a REQUEST: the driver clamps it to the team the OpenMP environment
#'     hands out (`omp_get_max_threads()`, `OMP_THREAD_LIMIT`, and the check
#'     farm's core cap), because the per-cell block cache sizes its slot array
#'     from that same number. A job script exporting `OMP_NUM_THREADS=1` to pin
#'     the INNER threads therefore runs the outer grid serially whatever is
#'     asked for here; the fit says so once, in a message, and records the pair
#'     in `n_threads_outer_requested` / `n_threads_outer_realised`.
#'   * `tile_warm` (`TRUE`) -- when `n_threads_outer > 1` and a copy block is
#'     present, group outer cells into tiles sharing every axis except the copy
#'     coefficient `alpha`, solve one warm Tier-2 per tile from the centre
#'     pilot, and warm-start the rest from their tile pilot. Falls back to the
#'     single-tier path when no copy block / single tile / `n_threads_outer
#'     <= 1L`. `FALSE` recovers the pre-tiling behaviour (e.g. regression
#'     testing).
#'   * `prune` (`TRUE`), `prune_tol` (`1e-3`) -- cheap-pass screening. The
#'     driver sweeps the outer lattice running a short inner Newton per cell,
#'     each warm-started from the previous screened cell's quasi-mode
#'     (lattice-adjacent), computes a screening Laplace log-marginal,
#'     softmax-normalises, and skips the full inner Newton on cells whose
#'     normalised weight is `< prune_tol`. The neighbour-warm-start sweep keeps
#'     every cheap mode near its cell's true mode, so the cheap ranking is
#'     faithful to the full-solve ranking even when the inner latent mode moves
#'     substantially across the grid. Pruned cells get `log_marginal = -Inf`,
#'     `n_iter = 0`, and inherit the pilot mode; the pilot cell is never
#'     pruned, and at least `prune_min_keep` cells (5) are solved in full
#'     whatever the tolerance says -- the highest-ranked dropped cells are
#'     restored up to that floor, because the outer grid placement pass reads a
#'     finite-difference curvature off the cells that were SOLVED and a kept set
#'     of one carries none. After the full pass the kept cells say how far off
#'     the screen was, and the posterior mass the dropped cells may carry is
#'     bounded by lifting each by the worst error measured on a kept cell
#'     (`prune_dropped_mass_bound`). Past 1% the driver solves the dropped cells
#'     that contribute most to the bound until it is under 1%
#'     (`prune_n_repaired`), so what the fit returns holds that bound. Two
#'     fallbacks to the full grid remain, each with a warning: a bound the
#'     repair could not bring under 1%, and a screened fit whose grid placement
#'     declined for want of curvature. On by default: a grid placed around its
#'     mode crosses axes the posterior does not fill (gcol33/tulpa#925), and
#'     solving every one of those cells is most of a large fit's cost. `FALSE`
#'     solves every cell, and is the default when the fit asks for local CCD
#'     refinement (`local_ccd`, or `k_refine = "ccd"`), whose stencil reads each
#'     refined cell's axis neighbours. Stacks with `n_threads_outer`.
#'     `prune_tol` must be in `[0, 1)`.
#'   * `prune_log_gap` (unset) -- the screening cut in nats: keep every cell
#'     within this many nats of the best screened cell. `prune_tol` is a
#'     normalised weight, so what it cuts at is the gap
#'     `-(log(prune_tol) + log(Z))`, and the whole range a caller would type is
#'     a few tens of nats -- on an outer surface spanning thousands every
#'     setting returns the same kept set. It reaches the kernel as
#'     `prune_tol = exp(-gap)`, so the realised cut is the requested gap less
#'     `log(Z)`; the fit reports the realised value in `prune_log_gap_cut`.
#'     Setting both this and `prune_tol` is an error.
#'   * `screen_iters` (`5L`) -- inner Newton steps the cheap screen takes per
#'     cell, read only when `prune = TRUE`. The screen only has to RANK the
#'     cells, and each one is warm-started from its already-screened lattice
#'     neighbour, so its quasi-mode is near the cell's own mode after very few
#'     steps. Every step is paid on every cell of the grid, including the cells
#'     the screen keeps and then solves in full, so a depth above what the
#'     ranking needs makes screening cost more than the solves it avoids. Must
#'     be a single integer `>= 1`.
#'   * `fitted_var` (`TRUE`) -- also fill `fitted_eta_var`, the per-cell,
#'     per-observation predictive variance of the linear predictor. It is the
#'     WITHIN-cell spread a grid-mixture replicate is drawn with
#'     ([posterior_predict()]); the per-cell predictor `fitted_eta` itself is
#'     stored either way. Computed by a per-cell solve sweep on top of the inner
#'     Newton, so a fit that will not be predicted from can decline it; the
#'     replicates of a fit run without it then carry the across-cell spread
#'     only. Read on a single-arm fit, which is the one that carries a single
#'     linear predictor to report.
#'   * `x_init` (`NULL`) -- warm-start for the first grid point's inner solve.
#'   * `verbose` (`FALSE`) -- when `TRUE`, announce the engaged outer integrator
#'     for a multi-block prior in one line at selection time (see `integration`),
#'     and report each CCD decline reason.
#'   * `store_Q` (`FALSE`) -- also return the per-grid joint precision Q (lower
#'     triangle, CSC) as `Q_csc_p_per_grid`, `Q_csc_i_per_grid`,
#'     `Q_csc_x_per_grid`, `Q_csc_n`, letting callers compute INLA-style
#'     total-variance posterior moments (`Var-of-means + Mean-of-Var`) on inner
#'     latent coordinates such as fixed-effect betas.
#'   * `keep_grid_hessians` (`TRUE`) -- retain the per-grid-point fixed-effect
#'     mode and marginal precision as `$grid_modes` / `$grid_hessians`, which is
#'     what lets [summary.tulpa_fit()], [confint.tulpa_fit()] and
#'     [vcov.tulpa_fit()] report the grid-marginalized fixed-effect covariance
#'     instead of `NA`. Memory is `O(n_fixed^2)` per cell. When the retention is
#'     not available the reason is recorded on `$grid_fixed_declined`.
#'   * `adaptive_grid` (`FALSE`), `adaptive_grid_edge_thresh` (`0.02`),
#'     `adaptive_grid_max_passes` (`1L`) -- when `adaptive_grid = TRUE`, a
#'     mode-tracked 1D refinement pass triggers on any axis whose marginal
#'     boundary weight exceeds `adaptive_grid_edge_thresh`. New points are
#'     added to that axis as levels (interior densification + outward
#'     log-spaced extension), each laid in every row of the other axes that
#'     holds the posterior -- at most `n_new_points * n_rows` kernel solves. The
#'     edge score is `max(marginal_weight_at_boundary, exp(max_log_marginal_at
#'     _boundary - max_log_marginal_overall))`, catching both boundary pile-up
#'     and integrand truncation; `0.02` is ~4 log units of decay.
#'     `adaptive_grid_max_passes` caps the passes (one usually suffices). Fixes
#'     posterior CI under-coverage when truth sits near a grid edge. Both the
#'     single-block and the multi-block driver run the passes on a tensor grid;
#'     a CCD design, the adaptive lattice and a locally CCD-refined grid keep the
#'     nodes they were laid on.
#'
#'     How far a given axis may be moved depends on where its nodes came from.
#'     An axis whose nodes the CALLER wrote down -- a block's `sigma_grid`,
#'     `copy$alpha_grid` / `field_coef$grid`, an entry of `phi_grid` -- states
#'     where the fit integrates, so refinement densifies it and never places a
#'     node past either end node; a mode outside that range shows up as mass at
#'     the edge.
#'     An axis the ENGINE placed (`alpha_n` / `field_coef$n`, or no nodes given)
#'     carries no such statement, so refinement may follow the posterior out
#'     past its ends. Nodes a wrapper package computed as a default of its own
#'     are declared with [auto_grid()] and count as engine-placed.
#'   * `axis_refine` (`NULL`) -- per-axis override of the rule above, named by
#'     outer-grid axis (`"sigma"`, `"alpha"`, `"phi_<arm>"`; `"b<k>.sigma"`,
#'     `"b<k>.tau"` and `"b<k>.alpha"` on a multi-block grid): `"none"` (the
#'     axis takes no refinement nodes), `"densify"` (nodes inside the declared
#'     span only) or
#'     `"extend"` (nodes past the end nodes too). A single unnamed value applies
#'     to every refinable axis. An unknown axis name, or a request for nodes on
#'     an axis this driver does not refine, is an error rather than a silent
#'     no-op. Read only when `adaptive_grid = TRUE` or
#'     `var_of_means_consistency = TRUE`.
#'   * `var_of_means_consistency` (`TRUE`) -- run a post-integration
#'     consistency pass over the refinable axes (a block's field SD, the copy
#'     scale and the per-arm dispersions). An axis whose marginal has collapsed
#'     onto too few nodes to carry a spread, such as a field SD narrower than
#'     its placed cell, has levels added in every row of the other axes that
#'     holds the posterior until its marginal resolves, so the reported interval is read off the
#'     posterior and not off a cell box. Attaches
#'     `var_of_means_consistency_info`.
#'   * `inner_factorization` (`"auto"`) -- which factorization the dense inner
#'     Newton applies to the Hessian it assembled: `"auto"` by the latent
#'     dimension, `"sparse"` for CHOLMOD, `"dense"` for the dense Cholesky.
#'     Independent of `force_sparse`, which selects the assembly path.
#'   * `force_sparse` (`FALSE`) -- linear-algebra backend for the inner joint
#'     solve. `TRUE` / `FALSE` select the sparse or dense path outright,
#'     regardless of the dense/sparse heuristic. `"auto"` selects by the latent
#'     dimension the fit will actually build, taking the sparse path above
#'     `n_x > 1000` and the dense one at or below it, where the sparse
#'     symbolic-analysis and indirection overhead outweighs the fill-in saving.
#'     A model whose latent dimension cannot be determined resolves to dense.
#'   * `integration` (`"auto"`) -- outer-grid node layout for a multi-block
#'     prior. A central composite design (CCD) integrates the hyperparameter
#'     posterior on `1 + 2d + 2^d` nodes oriented by the Cholesky of the
#'     posterior covariance at the joint hyperparameter mode -- far fewer inner
#'     solves than the `d`-dimensional tensor product (`25` vs `81` at `d = 4`).
#'     `"auto"` (the default) uses the CCD only at `>= 4` transformable axes,
#'     where the tensor product's `k^d` blow-up bites hardest, and keeps the
#'     cheaper, more ridge-robust tensor grid at `<= 3` axes; `"ccd"` lowers the
#'     CCD threshold to `>= 3` axes; `"grid"` always forces the full tensor
#'     product. `"grid_adaptive"` is the low-dimensional companion to the CCD:
#'     it seeds a coarse subsample of the SAME tensor lattice (latent block axes
#'     and phi axes together), floods outward from the posterior mode on the fine
#'     lattice, and evaluates only the cells within a log-density cutoff of the
#'     peak -- a strict, uniform-weight subset of the dense tensor, so its
#'     posterior matches the dense grid to that cutoff at fewer inner solves when
#'     the hyperparameter posterior concentrates (a sharply-identified field SD /
#'     precision). It declines back to the dense tensor on a diffuse posterior
#'     (kept region would rival the tensor) or a degenerate lattice, so it never
#'     costs accuracy; tune it with `adaptive_grid_cutoff` / `adaptive_grid_stride`
#'     / `adaptive_grid_max_frac`. The CCD auto-falls back to the tensor grid for
#'     an axis whose support is not safely transformable (a CAR_proper `rho_car`
#'     or a non-BYM2 `rho`), or a flat / ridged / degenerate outer mode-find or
#'     Hessian; an active `phi_grid` rides as a tensor axis crossed on top of the
#'     CCD.
#'     Single-block joint priors always use the tensor grid. The CCD mode-find
#'     runs cheap warm-started inner solves and does not write to the checkpoint
#'     file. Under `verbose = TRUE` the engaged integrator is announced in one
#'     line at selection time (e.g. `outer integration: CCD (4 latent axes, 25
#'     nodes)` or `tensor grid (72 cells)`, or `CCD declined -> tensor grid`),
#'     so the auto switch to the CCD at `>= 4` axes is never silent. The
#'     resolved integrator is also returned on the joint result as
#'     `$integration`, alongside `$integration_requested` and
#'     `$integration_declined` (see Value), so a fallback is on the fit itself
#'     and not only in a verbose message.
#'   * `ccd_budget` (`1`) -- ceiling on what PLACING a CCD may spend, as a
#'     multiple of the tensor grid the design replaces. Every mode-find probe is
#'     one full inner Newton solve, the same solve an outer grid cell costs, and
#'     a placement round is `1 + 2d + 4 C(d, 2)` of them (33 at `d = 4`), so at
#'     the round caps alone a placement can spend far more than the integration
#'     it avoids. At `1` a placement plus its design costs no more than that
#'     tensor grid; the fit declines to the tensor grid with
#'     `integration_declined = "placement_budget"` when the ceiling is reached,
#'     and reports what it spent as `ccd_modefind_evals` (see Value). `Inf`
#'     places without a ceiling. `ccd_budget_floor` (`TRUE`) exempts the
#'     designed first attempt -- the eight rounds from the grid-median seed,
#'     whose cost is set by the axis count and not by the grid -- so the ceiling
#'     bounds the escalation (grid seed, step calibration, the 30-round rescue
#'     mode-find) rather than the path the CCD was built for; `FALSE` makes
#'     `ccd_budget` a hard cap on the whole placement.
#'   * `local_ccd` (`NULL`) -- local CCD refinement of a multi-block tensor grid.
#'     `TRUE` (defaults) or a `list(max_cells =, f0 =, skew_max =, rank =)` refines a few high-weight,
#'     mutually non-adjacent interior cells, replacing each with a small
#'     curvature-aware CCD node cloud so a coarse base grid resolves the
#'     sharply-peaked directions without the `k^d` tensor blow-up. The local
#'     curvature is a diagonal finite difference of the outer log-marginal over
#'     the cell's own grid neighbours (no mode-find; only the off-centre nodes are
#'     new solves), warm-started from the cell's inner mode; each refined cell's
#'     sub-nodes carry partition-of-unity design weights so the total integration
#'     weight is conserved (no double-count). `max_cells` (`8L`) caps the refined
#'     cells, chosen in order of `rank`: `"weight"` (default), the cell's
#'     integration weight, or `"mass_moved"`, the weight times
#'     `|exp(log_box_ratio) - 1|` the box-mass rule predicts the cell's midpoint
#'     atom misplaces; `f0` (`1.1`) is the CCD radius. The design scale is shrunk per cell
#'     so the cloud fits the cell's Voronoi box (the local-Gaussian mass beyond it
#'     belongs to the neighbouring cells). A cell keeps its cloud only while the
#'     nodes' own log-marginals stay within `skew_max` (the `gamma3_ok` band,
#'     `0.5`) of the quadratic the cloud was placed from, measured as a
#'     standardized cubic magnitude; above it the cell is put back as its own
#'     mass atom, which on a skewed outer target is measurably closer than the
#'     design. Engages only on the tensor path (the
#'     curvature stencil needs axis neighbours), at `>= 4` transformable latent
#'     axes, with no active `phi_grid`; otherwise it is a no-op. The applied
#'     refinement is summarised on the result as `$local_ccd_info`. Also driven
#'     automatically by `k_refine = "ccd"`.
#'   * `adaptive_grid_cutoff` (`10`), `adaptive_grid_stride` (`2L`),
#'     `adaptive_grid_max_frac` (`0.75`), `adaptive_grid_min_cells` (`48`) --
#'     tuning for `integration = "grid_adaptive"`. `adaptive_grid_cutoff` is the
#'     log-density keep / expand radius from the peak (larger keeps more cells,
#'     closer to the dense tensor); `adaptive_grid_stride` the coarse-seed
#'     subsample stride per axis; `adaptive_grid_max_frac` the kept-fraction
#'     ceiling past which the builder declines back to the dense tensor;
#'     `adaptive_grid_min_cells` the smallest dense tensor worth the adaptive
#'     machinery -- below it the builder declines BEFORE any inner solve (on a
#'     small tensor the coarse seed is already most of the grid, so there is no
#'     tail to skip). Ignored by the other integrators. The kept-cell / dense /
#'     solve counts are returned as `$adaptive_grid_info`.
#'   * `hessian` (`"auto"`) -- the curvature each inner Newton step is built
#'     from, where a cell-coupled arm makes the joint Hessian indefinite away
#'     from the mode. `"auto"` takes the observed Hessian wherever it factors as
#'     it stands, the Newton step, and the coupled cells' complete-data expected
#'     information (Fisher scoring, positive semi-definite) where it does not;
#'     `"lm"` always takes the observed Hessian, loading its diagonal until it
#'     factors; `"fisher"` always takes the expected information, which
#'     converges linearly; `"psd"` eigen-clamps the observed Hessian (small
#'     `n_x` only). The final mode-pass behind `log_det`, the log marginal and
#'     the SEs uses the observed Hessian under every choice, and a fit with no
#'     coupled arm runs `"auto"` and `"fisher"` as `"lm"`.
#'   * `inner_refresh` (`1L`) -- inner-Newton Cholesky factor reuse interval
#'     (Shamanskii / chord method). For a non-quadratic positive arm (e.g. a
#'     beta cover arm) the latent Hessian changes every inner iteration, so the
#'     default re-factorizes the sparse Cholesky on each step -- the dominant
#'     per-grid-cell cost. `inner_refresh = m > 1` re-factorizes only every
#'     `m`-th inner step and reuses the cached factor in between (refreshing
#'     early whenever a reused solve fails). The gradient is exact on every step
#'     and each step is line-search safeguarded, so the converged mode is
#'     unchanged and the final mode-pass Hessian (`log_det`, SEs) is always
#'     fresh; only the path to the mode uses a stale curvature, which may cost a
#'     few extra inner iterations. Applies to the sparse joint path with an
#'     LM-guarded step (`hessian` `"auto"`, `"lm"` or `"fisher"`); the dense small-`n_x` path
#'     re-factorizes a cheap Hessian and ignores it. `2L`-`4L` is a good range
#'     for a slow beta arm.
#'   * `k_quality` (`"report"`) -- the reliability intent for the outer
#'     Pareto-\eqn{\hat{k}}, a single statement of how reliable the fit should be.
#'     `"report"` (default) computes the diagnostic and reports the achieved band.
#'     `"ok"` / `"good"` additionally name a TARGET band (the \eqn{\hat{k}}
#'     confidently usable, resp. good) and raise the default `k_samples`
#'     (to `800L` / `2000L`, unless you set it) so the bootstrap CI can resolve it.
#'     `"none"` disables the diagnostic. The fit carries an honest verdict --
#'     `k_quality_requested`, `k_quality_reached`, `k_quality_best`,
#'     `k_quality_miss`, `k_quality_reason`, `k_quality_rounds`,
#'     `k_quality_k_trace` -- and never silently downgrades: if the requested band
#'     is not confidently met it reports the band actually reached and why.
#'     For `"ok"` / `"good"`, when the first fit does not reach the band the
#'     engine escalates by REFINING THE INTEGRATION GRID (see `k_refine`),
#'     widening / densifying it where the posterior mass escapes its current
#'     bounds and re-diagnosing, up to `k_max_rounds` rounds. Refinement is the
#'     first lever on either kind of miss because it lowers \eqn{\hat{k}}
#'     itself rather than the uncertainty around it. Once refinement is
#'     exhausted -- no boundary mass left to chase -- what remains depends on
#'     the miss (`k_quality_miss`): a `"resolution"` miss (\eqn{\hat{k}}
#'     confidently outside the band) ends the chase, while a `"precision"` miss
#'     (the bootstrap CI straddling a boundary, so \eqn{\hat{k}} may already be
#'     inside the band) spends its remaining rounds doubling the
#'     importance-draw budget, the one lever refinement does not touch.
#'   * `k_refine` (`"grid"`) -- the integration-refinement rung for `k_quality`
#'     `"ok"` / `"good"`. `"grid"` (default) re-fits with adaptive grid refinement
#'     (`adaptive_grid`) each escalation round, driven by the bad \eqn{\hat{k}}, so
#'     a too-coarse / too-narrow grid is widened / densified where the importance
#'     weight concentrates until the band is reached or the budget is spent.
#'     `"ccd"` instead refines a few high-weight, mutually non-adjacent interior
#'     cells with local curvature-aware CCD node clouds (see `local_ccd`), the
#'     right rung when the grid is too coarse to resolve a sharply-peaked
#'     direction rather than too narrow; it forces a tensor base grid (the
#'     curvature stencil needs axis neighbours) and engages only on the
#'     multi-block path at >= 4 transformable latent axes. `"none"` disables
#'     refinement: the band is reported but not chased.
#'   * `k_max_rounds` (`2L`) -- the escalation round budget for `k_quality`
#'     `"ok"` / `"good"`: the maximum number of re-fit rounds after the first fit.
#'     A refinement round allows one more pass than the last refinement round
#'     did; a round taken after refinement is exhausted leaves the grid where it
#'     stands and doubles the importance-draw budget instead. `0L` disables
#'     escalation (single-shot, the band is reported but not chased).
#'   * `diagnose_k` (`TRUE`), `k_samples` (`500L`) -- compute the outer
#'     Pareto-\eqn{\hat{k}} accuracy diagnostic by importance-sampling the joint
#'     hyperparameter posterior against the proposal the integrator fits (mixed
#'     per-axis transforms: `log` for positive scales, logit for the BYM2 mixing
#'     weight, identity for the copy coefficient \eqn{\alpha}). `k_samples`
#'     is the number of importance draws, each one an extra inner joint solve, and
#'     is where more tail information comes from: a tighter k-hat needs MORE
#'     actual tail ratios, so increase `k_samples` (not `k_bootstrap`). It is not
#'     a pure precision knob, though. Under the automatic PSIS tail rule
#'     (`min(S/5, 3 sqrt(S))`) the fitted tail FRACTION shrinks as `3 / sqrt(S)`,
#'     so a larger budget describes a DEEPER quantile of the weight distribution
#'     and the reported k-hat can move materially -- measured on a synthetic
#'     heavy-tailed outer target, 0.57 at 500 draws and 7.94 at 50000
#'     (gcol33/tulpa#631). Set `k_tail_points` in proportion to `k_samples` to
#'     hold the fraction and get a strictly sharper estimate of the same number;
#'     the `k_quality` escalation's draw rung does exactly that. The draws are
#'     RNG-restored so the fit's modes / draws are unchanged. A fit carrying an
#'     axis whose support is not safely known (CAR_proper's `rho_car`) declines to
#'     the quadrature-ESS fallback (`pareto_k = NA`). `FALSE` skips the diagnostic.
#'     `"by_arm"` additionally computes a k-hat restricted to each arm's
#'     hyperparameter axes (the other arms held at their posterior mean), reported
#'     in `pareto_k_by_arm`, to localise which arm drives a tail-heavy joint k; the
#'     joint k itself is unchanged. Per-arm k is defined for the multi-block layout
#'     with two or more arms and declines for the single-block shared-field layout.
#'   * `k_threads` (`NULL`) -- outer-thread width for the diagnostic's importance
#'     batch. The `k_samples` re-solves are independent and run after the grid
#'     (every core free), each solved single-threaded once the batch saturates the
#'     pool, so widening it is a bit-identical wall-clock speedup (the k-hat is
#'     unchanged). `NULL` follows the fit's own thread grant -- the larger of
#'     `n_threads_outer` and the inner `n_threads` -- so a serial fit keeps a serial
#'     diagnostic while a threaded fit gets a free parallel one. `"auto"` uses the
#'     physical performance-core count (capped at 2 under R CMD check); an integer
#'     pins the width (`1L` forces serial). Always capped at `k_samples`.
#'   * `k_bootstrap` (`1000L`) -- bootstrap replicates for the outer
#'     Pareto-\eqn{\hat{k}} uncertainty. The k-hat is a single fixed number for a
#'     fit + proposal; its sampling uncertainty GIVEN the proposal is estimated by
#'     resampling the diagnostic's raw importance log-ratios with replacement and
#'     re-fitting the GPD tail `k_bootstrap` times (no new inner solves). Reports
#'     `pareto_k_se_boot` (bootstrap SE), `pareto_k_ci_low` / `pareto_k_ci_high`
#'     (the 2.5\% / 97.5\% bootstrap quantiles), `pareto_k_se_formula` (the
#'     closed-form GPD-shape MLE asymptotic SE \eqn{(1 + k)/\sqrt{M}}, a
#'     cross-check), and `pareto_k_band_confident` (TRUE iff the bootstrap CI lies
#'     within one reliability band). The bootstrap measures how UNSTABLE the
#'     current tail estimate is; it cannot create tail information. Increase
#'     `k_samples`, not `k_bootstrap`, to obtain more tail information. `0L`
#'     skips it (point k-hat only). The per-arm k carries the same fields.
#'   * `k_tail_points` (`NULL`) -- number of upper-tail order statistics for the
#'     GPD fit. `NULL` uses the automatic PSIS rule
#'     \eqn{\lceil\min(0.2 N, 3\sqrt{N})\rceil}. An explicit value is an EXPERT
#'     tail-threshold control, capped at the 20\%-of-draws ceiling (with a warning)
#'     so the fit stays an extreme tail; it is NOT a precision knob (a request that
#'     drags body ratios into the tail lowers variance but biases the k-hat). The
#'     used and requested counts are reported in `pareto_k_tail_points` /
#'     `pareto_k_tail_points_requested`.
#'   * `k_conf_bands` (`NULL`) -- the reliability-band boundaries the bootstrap CI
#'     is tested against for `pareto_k_band_confident`. `NULL` (default) uses the
#'     sample-size-dependent boundaries \eqn{c(0.5, \min(1 - 1/\log_{10} S, 0.7))}
#'     at the realised draw count \eqn{S} (Vehtari et al. 2024): the good cut is
#'     0.5 and the usable cut tightens below 0.7 for small \eqn{S} (about 0.565 at
#'     \eqn{S = 200}, reaching 0.7 only past \eqn{S \approx 2154}). Supply a
#'     strictly-increasing numeric vector to fix the boundaries instead, e.g.
#'     `c(0.5, 0.7)` for the size-independent good / ok / unreliable split.
#'   * `diagnose_skew` (`TRUE`), `skew_idx` (`NULL`) -- compute the
#'     inner-Laplace skewness diagnostic (`$inner_skew`, gamma_3, Rue Martino
#'     & Chopin 2009 Sec 3.2.3) at the fitted MAP grid cell: one extra Newton
#'     solve via the same `kernel_fn` the outer diagnostic reuses, scoring
#'     every arm's fixed-effects coefficients by default (pass `skew_idx`,
#'     1-based indices in the joint `[arm1_beta | arm1_re | arm2_beta | ... |
#'     blocks]` latent layout -- see `$arm_layout` -- to probe additional
#'     indices). This is the complementary layer to `diagnose_k`: that scores
#'     the outer hyperparameter-grid integration around a FIXED inner
#'     Laplace, this scores whether that inner Gaussian approximation is
#'     itself a good fit. A genuinely coupled arm (`cell_coupling !=
#'     "separable"` on that arm, e.g. tulpaObs's `occu_cover`) has no per-obs
#'     likelihood for the separable formula to score; it is scored instead by
#'     the contraction of the cell third-derivative tensor, which differences
#'     the cross-arm Hessian the coupling spec already returns. A fit that can
#'     carry neither reports `NaN`, not a silently wrong 0, and
#'     `$inner_skew_declined` says WHY -- `"coupled_arm"` marks arms the inner
#'     layer could not score at all, distinct from a diagnostic that was simply
#'     switched off -- while `$inner_skew_arms_declined` names the arms with no
#'     oracle, including on a partially scored fit. See
#'     [diagnostics()] for the combined verdict. Wired
#'     for both the single-block backends (icar/bym2/car_proper) and the
#'     multi-block path (a per-group RE, a trend field, or an arm-specific
#'     field block).
#'   * `within_cell` (`"log_quadratic"`) -- the WITHIN-CELL construction the
#'     reported per-axis hyperparameter intervals are read with. 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.
#'     `"log_quadratic"` joins each row's node log densities by overlapping
#'     quadratics, exact for a Gaussian at any spacing, and continues the end
#'     quadratic past the outer node instead of stopping at the box edge; it
#'     runs on the rows that resolve their own conditional (`h / sd <= 2`) and
#'     reads every other row as boxes, declining to `"box_uniform"` where no
#'     row qualifies. `"box_uniform"` puts the cumulative FULL mass at each
#'     cell EDGE and interpolates between edges; `"chord"` puts the cumulative
#'     MID-mass at each cell coordinate and interpolates between coordinates.
#'     The default moved from `"box_uniform"` at gcol33/tulpa#932: a placed
#'     axis sits at 1.25 posterior SDs per cell, where the box read is narrow
#'     by construction, and against dense references of a two-arm ICAR
#'     fixture the log-quadratic read's summed error is 0.038 against the box
#'     read's 0.110 on the same grids. `outer_grid_h_over_sd` is how wide a
#'     cell is on each axis, and `theta_within_cell` is what each axis was
#'     actually read with. Only a `"density"` support admits the box and
#'     log-quadratic reads -- a CCD design, a locally refined grid and a
#'     posterior sample are not cell partitions that tile -- and an axis they
#'     decline on reports the read that ran with a reason rather than erroring.
#'     Nothing else moves: point estimates, moments, draws and weights are
#'     untouched.
#'   * `skew_correct` (`TRUE`) -- consume the inner-Laplace expansion instead of
#'     only grading it: report
#'     Cornish-Fisher marginal quantiles at each coefficient's own `gamma_3`,
#'     about the centre `gamma_1 + gamma_3 / 2`, from `summary()` / `confint()`
#'     wherever the combined inner band (`gamma_3` and the inner importance
#'     k-hat) says the leading-order expansion is in its regime, and the Gaussian
#'     quantiles everywhere else. It is post-processing on the reported
#'     quantiles: draws, modes and weights are untouched, so a fit run with it
#'     off is bit for bit the fit it was before. The band that bounded the
#'     relocation itself (`centre_unreliable`) is off (`Inf`):
#'     every finite cutoff was measured to decline the coefficients the
#'     correction helps most.
#'     `$skew_correction` records the per-coefficient `gamma_3`, `gamma_1` and
#'     the centre they form, the band, the k-hat, the combined
#'     band, the eligibility and the reason behind it; the `skew_applied`
#'     attribute on `summary()` / `confint()` records what was actually used at
#'     the requested level. A FULLY COUPLED fit declines it: the location term's
#'     contraction against a covariance block is not reachable from the cell
#'     third-derivative oracle, and an absent `gamma_1` is never read as zero.
#'     Such a fit therefore reports what it reported before the default moved, to
#'     the bit -- a declined coefficient keeps the grid-mixture read. See
#'     [tulpa_nested_laplace()] for the measurement behind the default.
#'   * `subspace_debias` (`FALSE`) -- correct only the latent directions the
#'     inner-layer diagnostics flagged, by exact Metropolis along the
#'     Gaussian-conditional-mean surface through each cell's mode, and leave the
#'     rest at their Gaussian conditional. Settings and semantics are the ones
#'     documented on [tulpa_nested_laplace()]. On a fully coupled fit `gamma_3`
#'     is `NaN` for every arm, so the selector rests on the derivative-free
#'     inner Pareto-k-hat; where that also bands the
#'     coordinate reliable, `idx` pins the set explicitly. An EMPTY selection
#'     leaves the fit bit-for-bit identical to the plain path. A fit whose
#'     inner solve took the s2z rank-1 or the PSD eigen-clamp path carries no
#'     usable factor to build the surface from and is left uncorrected, the same
#'     two paths `diagnose_skew` declines on.
#'   * `cila` (`FALSE`) -- corrected integrated Laplace, the second inner-layer
#'     debias (after Lai, Margossian and Sheldon,
#'     arXiv:2605.20345). Where `subspace_debias` selects coordinates and runs
#'     exact Metropolis on them, this selects nothing: at every outer cell it
#'     draws `n_points` points from the whole inner Gaussian, weights each by the
#'     exact joint density it came from, and reports the weighted particles. The
#'     cell marginals and the latent posterior both converge to the exact ones as
#'     the effort grows, so `n_points` is the only dial. `TRUE` takes the
#'     defaults; a list overrides `n_points` (`1024L`), `variant`
#'     (`"qmc"`, a Sobol net; `"is"` for iid draws, `"rqmc"` for the net under
#'     `n_shift` random shifts), `n_shift` (`8L`), `n_draws` (the reported draw
#'     count) and `seed` (the auxiliary stream, engine-owned so requesting the
#'     correction leaves every other posterior draw unchanged). Below 512 points
#'     a cell's particle set is too coarse to be a marginal at all and the
#'     request is refused. The correction reports DRAWS, so every
#'     coefficient-facing method reads them instead of the grid's Gaussian
#'     mixture, and the corrected per-cell masses become the fit's own
#'     `weights` / `log_marginal` with `weights_source` reporting `"cila"`; the
#'     pre-correction pair is kept as `$cila$laplace` and `$cila$retained_mass`
#'     is the original share of the cells that produced a usable particle set.
#'     `$cila` also carries the variant actually run and the PSIS grade
#'     (`pareto_k`, `rel_ess`) of the correction's own weights. A cell whose
#'     inner solve factorized SPARSELY draws through the CHOLMOD factor's own
#'     triangular and permutation solves; an LDL' factor carries no square root
#'     to draw with and is declined with `"sparse_factor_not_ll"`.
#'   * `auto_recenter` (`TRUE`) -- re-centre a default outer grid axis on its
#'     posterior mode and refit when the fit rails against a boundary node.
#'     `FALSE` integrates over the grid exactly as given, whatever it is, and
#'     records `outer_grid_recenter_declined = "auto_recenter_disabled"`. The
#'     joint FIELD rescues trigger on the whole grid's collapsed-edge regime
#'     rather than on a per-axis rail, so the per-axis policy names
#'     [tulpa_nested_laplace()] takes (`"rail"`, `"resolve"`, `"always"`) are
#'     refused here with an error rather than accepted and ignored. A per-arm
#'     dispersion axis (`phi_grid`) is the exception and fires on its own
#'     sizing: it is crossed onto the tensor independently of the field's
#'     geometry, so whether the field's grid collapsed says nothing about
#'     whether the dispersion axis brackets its own posterior. The copy scale
#'     `alpha` fires the same way when its grid is the engine's own or a
#'     wrapper's `auto_grid()` default: its continuum is laid at the mode and
#'     its "no coupling" atom at 0 kept, and an axis whose posterior sits on
#'     that atom is left where it is. Every axis that fires is laid from the ONE
#'     outer mode the placement mode-find reached, and the grid is refit once
#'     with all of them (gcol33/tulpa#925). `FALSE` holds these axes too.
#'   * `recenter_pilot` (`FALSE`) -- detect the placement above on a THINNED
#'     grid rather than on the full one. Placement reads two things, the argmax
#'     cell and an FD curvature stencil at it, and reads both off
#'     `log_marginal`; neither is read off the integration the detecting pass
#'     paid for, so every inner Newton solve of that pass is discarded the
#'     moment a placement fires. With the pilot on, the detecting grid keeps
#'     each axis's endpoints and at most `TRUE`'s three nodes between them (an
#'     integer names its own resolution), and the full grid is solved once, at
#'     the placed axes.
#'
#'     Off by default because the trade is a property of the WORKLOAD. With `F`
#'     full cells, `P` pilot cells and a firing rate `p`, the un-piloted cost is
#'     `F + p (S + F)` against the pilot's `P + p S + F`, so the pilot pays
#'     exactly when `P < p F`. Measured on a 3-axis (sigma, alpha, phi) donor +
#'     copy fixture over 48 paired seeds: 0.78x the cells solved on the
#'     configuration that places on 10 of 12 seeds, and 1.29x on the one that
#'     never places.
#'
#'     A coarse detector is not the same detector. The collapse trigger is
#'     `ess_grid < 2`, which falls with the cell count, so a pilot fires
#'     somewhat more readily than the full grid: on that same sweep the two
#'     disagreed on 9 of 48 pairs, every one of them the pilot placing where the
#'     full grid did not, and 7 of the 9 reached a HIGHER maximum inner
#'     log-marginal after placing (+0.35 to +2.14 nats) while 2 lost 0.46 and
#'     0.72. It cannot go the other way: a declined pilot is followed by the
#'     full fit, which then gets its own detection, so the pilot only ever adds
#'     a placement. A fit that declines both is bit-identical to the same fit
#'     with the knob off. What the pilot detected on is recorded in
#'     `outer_grid_pilot`.
#'   * `max_grid_cells` (`NULL`, no ceiling) -- a ceiling on a multi-block
#'     tensor outer grid's cell count, refused before any inner solve with the
#'     axes that produced the count. The count is the number of inner solves,
#'     so it includes an active `phi_grid`: 2500 latent cells crossed with 4
#'     dispersion values are 10000 cells. It bounds the dense tensor, including
#'     the one a declined `"grid_adaptive"` falls back to. Unset, a tensor of
#'     any size is solved, and one past 2048 cells is announced with a warning
#'     first.
#'   * `checkpoint` (`NULL`) -- grid-cell checkpoint/resume. Set
#'     `list(path = "fit.ckpt", resume = TRUE)` to make a killed or interrupted
#'     fit resumable: each completed outer-grid cell is appended to `path`, and
#'     a later call with the same responses + grid + control loads the finished
#'     cells and solves only the rest. EVA-scale joint fits run for hours, so a
#'     wrapper teardown, reboot, or OOM near the end otherwise loses the whole
#'     run. `resume = TRUE` (the default when a `path` is given) continues from
#'     an existing file; `resume = FALSE` removes it first and starts over.
#'     Adaptive-grid refinement cells are checkpointed under their own
#'     coordinate keys, so resume covers them too. A file written for different
#'     data or solver settings is rejected (fingerprint mismatch) rather than
#'     resumed onto a stale result.
#'
#' @return A list of class `c("tulpa_nested_laplace_joint",
#'   "tulpa_nested_laplace", "list")` with:
#'   * `theta_grid`, `theta_names` -- outer-grid hyperparameter values
#'     (includes the `alpha` axis when `copy` is set).
#'   * `log_marginal`, `weights` -- per-grid-point log-marginal and integration
#'      weights (sum to 1).
#'   * `theta_mean`, `theta_sd` -- posterior moments per hyperparameter,
#'      including `alpha` when `copy` is set.
#'   * `theta_median`, `theta_ci_lo`, `theta_ci_hi` -- weighted-quantile
#'      median and 2.5/97.5 empirical CI per hyperparameter axis (same
#'      names as `theta_mean`). Recommended summary for right-skewed
#'      scale-like axes (alpha at small n_pos, sigma/range/phi near
#'      a boundary), where the posterior mean is pulled by the right
#'      tail away from the bulk and `mean +/- 1.96 sd` mis-states the
#'      uncertainty.
#'   * `theta_sd_source` -- per axis, what the reported SD was read off:
#'      `"weighted"` (the spread of the weights, on a resolved axis),
#'      `"stencil"` (a parabola at the modal node) or `"within_cell"` (the
#'      modal node's mass spread over its cell, the read the interval of an
#'      unresolved axis comes from). `theta_sd_stencil_declined` says why a
#'      parabola was not formed.
#'   * `modes` -- `[n_grid x n_x]` matrix of inner modes.
#'   * `fitted_eta` -- `[n_grid x N]` matrix of the linear predictor at each
#'      cell's own mode, offset and latent field included. Present on a
#'      single-arm fit, which is the one that has a single linear predictor to
#'      report; a cell that was pruned or whose block preparation failed reads
#'      `NA`. `fitted_eta_var` is its per-cell within-cell variance, present
#'      under `control$fitted_var` (the default). Together these are the
#'      per-cell Gaussian [posterior_predict()] draws a replicate from.
#'   * `n_iter` -- inner Newton iterations per grid point.
#'   * `arm_layout` -- list with per-arm `beta_start`, `re_start`,
#'      spatial offset(s) and `n_x` for decoding modes.
#'   * `prior`, `responses`, `copy` -- echoed inputs.
#'   * `timing` -- named numeric of wall-clock seconds: `total` plus the
#'      `setup` (validation / encoding / grid construction), `grid` (inner
#'      Laplace solves, including adaptive-refinement and consistency passes),
#'      `postproc` (weight / moment / marginal assembly) and `diagnostics`
#'      (outer Pareto-\eqn{\hat{k}}) phases. The `grid` phase is the one that
#'      scales with grid size and core count. Surfaced one-line in `print`.
#'   * `pareto_k`, `pareto_k_is_ess`, `pareto_k_scope` -- outer
#'      Pareto-\eqn{\hat{k}} accuracy diagnostic and its importance-sampling
#'      ESS (both `NA` when `control$diagnose_k = FALSE` or the fit declines;
#'      see the `diagnose_k` control knob). `pareto_k < 0.7` indicates the
#'      nested integration is reliable; `>= 0.7` that the (skewed / heavy-
#'      tailed) hyperparameter posterior is misfit by the Gaussian grid and
#'      the fit should escalate to an exact debias.
#'   * `pareto_k_declined` -- when `pareto_k` is `NA`, WHY:
#'      `"not_requested"` (`control$diagnose_k = FALSE`; nothing is wrong),
#'      `"unguessable_axis: <axis>"` (a support the engine will not guess, e.g.
#'      car_proper's `rho_car` -- a PERMANENT limitation of that family, so read
#'      the quadrature ESS instead), `"draws_too_few"`, `"grid_too_small"`,
#'      `"no_varying_axis"`, `"degenerate_proposal"` (the outer mode curvature
#'      came back non-finite -- a signal about the fit), `"not_applicable"`, or
#'      `"internal_inconsistency"` (an engine bug worth reporting). `NA` on a
#'      fit whose \eqn{\hat{k}} WAS computed.
#'   * `pareto_k_proposal_source` -- how the outer importance proposal the
#'      \eqn{\hat{k}} scores was built: `"mode_hessian"` from the Laplace
#'      curvature at the hyperparameter mode (the CCD design's, or a
#'      finite-difference Hessian when a sharp posterior collapses the grid),
#'      `"grid_moment"` from the grid-weighted covariance, `"moment_matched"`
#'      when the moment-matching refinement improved on either, `"grid_mixture"`
#'      when the faithful mixture-over-grid-cells proposal won, or
#'      `"skew_normal"` when the skew-normal rescue did. `NA` when the
#'      diagnostic is off or declines. The mode-Hessian source keeps the
#'      \eqn{\hat{k}} meaningful when the grid concentrates on ~1 cell.
#'   * `pareto_k_regime` -- whether the outer grid actually integrated
#'      hyperparameter uncertainty: `"spread"`, or `"collapsed_interior"` /
#'      `"collapsed_edge"` when the quadrature effective sample size falls below
#'      two cells so the integration has degenerated to a point evaluation at
#'      the modal hyperparameter. An interior collapse is benign (the grid
#'      bracketed the mode; the fit is empirical Bayes at that mode and only the
#'      integrated hyperparameter uncertainty is missing); an edge collapse means
#'      the mode sits at an extreme node, so the grid may be too narrow.
#'      Attached from the stored weights, so it is present even with
#'      `control$diagnose_k = FALSE`.
#'   * `pareto_k_grid_edge_axes`, `pareto_k_grid_edge_sides` -- for an edge
#'      collapse, which axes the dominant cell sits against and on which side
#'      (`"lower"` / `"upper"`); widen those axes and refit to confirm the mode
#'      is bracketed. Axes the grid pins to one value are excluded (pinned, not
#'      at a boundary).
#'   * `axis_span` -- per outer axis, the span the fit worked over and the two
#'      terms that separate it from the nodes it was given: `nodes` (the range
#'      of the axis's initial continuum nodes), `declared` (their support, i.e.
#'      those nodes plus the half node step the outermost cells own -- for `k`
#'      equally spaced nodes that is `k / (k - 1)` times the node range on the
#'      axis's integration coordinate, so 2x at two nodes and 1.125x at nine),
#'      `integrated` (the region the final grid's cell measure integrates, the
#'      same interval `axis_support` reports), `refine` (the mode refinement
#'      ran under: `"none"` / `"densify"` / `"extend"`) and `n_nodes` (initial
#'      and final continuum node counts).
#'   * `outer_grid_placement` -- `"fixed"` (the default `sigma_grid` axis was
#'      used as-is) or `"auto_recentered"` when a railed `sigma` axis
#'      triggered the mode-Hessian recenter-and-refit (see the `prior`
#'      argument above). `outer_grid_recenter_attempts` (integer)
#'      and `outer_grid_prior_added` (logical: whether the light default
#'      PC(U=3, alpha=0.01) sigma prior was applied) are set alongside it, plus
#'      `outer_grid_prior_declined = "prior_pinned"` when a second attempt ran
#'      but the caller's pinned `prior_sigma` held that prior back.
#'   * `outer_grid_recenter_declined` -- on a `"fixed"` placement, why the
#'      recenter did not run: `"grid_not_collapsed"` (the grid already brackets
#'      the mode -- the common case, no refit needed), `"axis_pinned"` (the
#'      caller pinned `sigma_grid`; mark it with [auto_grid()] if it is a
#'      default rather than a choice), `"default_axis_pinned"` (a wrapper
#'      package declared the nodes with `auto_grid(place = FALSE)` and asked
#'      for them as written -- nothing you pinned), or `"no_usable_curvature"`
#'      (the outer mode-find the recenter lays the axis from did not reach a
#'      mode with a usable curvature, e.g. a car_proper grid whose `rho_car`
#'      axis has unguessable support).
#'      Absent when the fit WAS recentred. A joint fit's field SD and its
#'      per-arm dispersion are placed by different passes over the same grid,
#'      and this slot reduces over them rather than holding the last to speak:
#'      an axis held because it was declared does not stand as the fit's answer
#'      while a pass with a movable axis has one.
#'   * `outer_grid_axis_declined` -- the same question PER AXIS, as a named
#'      character vector. The slot above holds one reason for the whole fit and
#'      is written only while the fit is unplaced, so on a fit where one axis
#'      moved and another did not it says `auto_recentered` and nothing about
#'      the axis that stayed -- which is the axis `grid_coarsest_axis` then
#'      names. Written by the field-SD passes for `sigma` (`b<k>.sigma` on a
#'      copy block) and by the per-arm dispersion pass for `phi_<arm>`, and
#'      carried across a later pass placing a different axis; absent on a fit
#'      no pass declined an axis of. An axis a pass DID move that still rails
#'      reads `"attempts_exhausted"` (or the reason its last attempt could not be
#'      laid), and is the axis the railed-mode warning names.
#'   * `outer_grid_pilot` -- present only when `control$recenter_pilot` ran: the
#'      pilot's resolution (`n_pilot`), its cell count (`cells`), the axes it
#'      thinned (`axes`) and those it could not (`axes_kept`), and what the
#'      placement was DECIDED from -- the pilot's own `regime`, `ess_grid` and
#'      `edge_axes`. The reported grid is never the grid the trigger was read
#'      on, so a two-pass fit is legible from the object rather than only from
#'      the progress log. `outer_grid_pilot_declined = "pilot_placement_declined"`
#'      records a pilot whose detection placed nothing and which was therefore
#'      followed by the full fit and the full fit's own detection.
#'   * `outer_grid_recenter_sd_clamp`, `outer_grid_recenter_sd_raw`,
#'      `outer_grid_recenter_sd_used` -- on a recentred fit, one entry per moved
#'      axis: which mode-SD bound the placement hit (`"none"` / `"floor"` /
#'      `"ceiling"`), the SD the finite-difference stencil MEASURED, and the SD
#'      the axis was actually laid from. A clamped axis was
#'      laid from a substituted spread rather than a measured one, and the
#'      reported interval is read off that span; the two used to be
#'      indistinguishable on the fit.
#'   * `pareto_k_outer_skew` -- per-axis skewness of the hyperparameter marginal
#'      in the proposal's whitened coordinate, present only when a
#'      \eqn{\hat{k}} above the good band triggered the skew-normal rescue pass.
#'      It is the EXPLANATION for an inflated \eqn{\hat{k}}: a symmetric Gaussian
#'      proposal against a right-skewed variance-component marginal has a heavy
#'      importance-ratio tail whatever the integration's quality. `NULL` when the
#'      Gaussian proposal already fit.
#'   * `pareto_k_se_boot`, `pareto_k_ci_low`, `pareto_k_ci_high`,
#'      `pareto_k_se_formula`, `pareto_k_tail_points`,
#'      `pareto_k_tail_points_requested`, `pareto_k_band_confident` -- the outer
#' Pareto-\eqn{\hat{k}} uncertainty, present whenever the
#'      diagnostic ran. `pareto_k_se_boot` is the bootstrap SE of the k-hat and
#'      `pareto_k_ci_low` / `pareto_k_ci_high` its 2.5\% / 97.5\% bootstrap
#'      quantiles -- the estimator's sampling spread GIVEN the proposal, not a
#'      posterior credible interval; `pareto_k_se_formula` is the closed-form
#'      GPD-shape MLE asymptotic SE cross-check. `pareto_k_tail_points` is the GPD
#'      tail size used and `pareto_k_tail_points_requested` the `k_tail_points`
#'      request (`NA` when automatic). `pareto_k_band_confident` is TRUE iff the
#'      bootstrap CI lies within one reliability band, `NA` when the bootstrap was
#'      off or could not fit. `diagnose_draws` and `diagnose_cost_ratio` (the
#'      diagnostic's draw budget and its wall-clock cost relative to the fit) are
#'      attached at the top level.
#'   * `pareto_k_by_arm`, `pareto_k_by_arm_is_ess`, `pareto_k_by_arm_scope` --
#' present only with `control$diagnose_k = "by_arm"`.
#'      Named (by arm) outer Pareto-\eqn{\hat{k}} restricted to each arm's
#'      hyperparameter axes, the other arms held at their posterior mean, so a
#'      tail-heavy joint k can be localised to one arm. A per-arm entry is `NA`
#'      when that arm carries no varying axis. Each arm carries its own bootstrap
#'      uncertainty: `pareto_k_by_arm_se_boot`, `pareto_k_by_arm_ci_low` /
#'      `pareto_k_by_arm_ci_high`, `pareto_k_by_arm_se_formula`,
#'      `pareto_k_by_arm_tail_points` and `pareto_k_by_arm_band_confident`.
#'   * `k_quality_requested`, `k_quality_reached`, `k_quality_best`,
#'      `k_quality_miss`, `k_quality_reason`, `k_quality_rounds`,
#'      `k_quality_k_trace` -- the reliability verdict for the
#' `control$k_quality` intent. `k_quality_requested`
#'      echoes the intent; `k_quality_best` is the band THIS fit reached
#'      (`"good"` / `"ok"` / `"unreliable"`, or `"uncertain"` when the bootstrap CI
#'      crosses a boundary) -- the band of the returned fit, not the best band
#'      seen while escalating, which `k_quality_k_trace` carries;
#'      `k_quality_reached` is `TRUE`/`FALSE` for an `"ok"` /
#'      `"good"` target (`NA` for `"report"` / `"none"`), never silently
#'      downgraded; `k_quality_miss` classifies a miss as `"resolution"` (the
#'      \eqn{\hat{k}} confidently outside the band -- a grid deficiency) or
#'      `"precision"` (the bootstrap CI straddling a boundary -- the estimator's
#'      variance), `NA` when there is no miss to classify, and is what the
#'      escalation reads once refinement is exhausted; `k_quality_reason`
#'      records why it stopped;
#'      `k_quality_rounds` is the number of escalation re-fits performed (`0` when
#'      the first fit sufficed or escalation was off); and `k_quality_k_trace` is
#'      the outer \eqn{\hat{k}} round by round, starting at the first fit's, so a
#'      chase that did not descend monotonically is visible. The refinement rungs
#'      rebuild the proposal from the grid they just changed, so a rising trace is
#'      a re-founded proposal on a strictly larger grid rather than a worse fit,
#'      which is why the LAST round is returned rather than the lowest-\eqn{\hat{k}}
#'      one.
#'   * `adaptive_grid_info` -- when `adaptive_grid = TRUE`, a list with
#'      `triggered_axes` (character) and `n_points_added` (integer)
#'      describing the refinement passes. NULL otherwise.
#'   * `local_ccd_info` -- when `local_ccd` engaged (or `k_refine = "ccd"`), a
#'      list with `n_cells_refined`, `n_nodes_added`, the refined `cells`, the
#'      `n_cells_before` / `n_cells_after` grid sizes, `n_design_nodes` (the cells
#'      that carry an in-cell design weight rather than their own mass),
#'      `design_mass` (the share of the integration weight sitting on them) and
#'      `cell_share` (the share each refined cell held on the BASE grid, before
#'      any node was placed). The two shares differ: the replacement nodes sit
#'      nearer the peak than the cell's own coordinate did, so refining raises
#'      the share the refined region holds, and reading both separates how
#'      concentrated the base grid already was from how much the refinement
#'      concentrated it. Three per-cell readings come with it, on three
#'      orthogonal axes, none of them gating anything. Shape: `misfit` (the
#'      standardized cubic magnitude of the part of the cell the design cannot
#'      represent, the one reading `skew_max` gates on, with `misfit_declined` /
#'      `cells_declined` for the cells put back as their own mass atoms).
#'      Centring: `offset` (the norm of the whitened gradient, i.e. how far the
#'      cell's own peak sat from the cell's coordinate in units of its marginal
#'      spread) and `mode_gain` (the same displacement in the cell's own
#'      curvature units, `0.5 * g' (-H)^-1 g` in nats, read off the same fit's
#'      quadratic coefficients and comparable across cells whose curvature
#'      differs; `NA` where `-H` is not positive definite, a cell with no
#'      interior peak to be displaced from). Mass: `log_mass_ratio` (the cell's
#'      refined mass over the coarse atom it replaced, `log_mass_refined -
#'      log_mass_coarse`, both carried on the log scale with the cell's outer
#'      design weight in them) and `max_node_weight` (the share the single
#'      largest node takes of its own cell's refined mass). Each carries a
#'      `_declined` twin for the cells the gate put back, whose nodes were
#'      evaluated before the score was read. `misfit` certifies that the design
#'      can represent the cell's shape and nothing more: a cell can be exactly
#'      quadratic and still be read at a point that is not representative of it
#'      (`offset` / `mode_gain`), and refining it can still move how much mass it
#'      competes for against its unrefined neighbours (`log_mass_ratio`, the
#'      embedded-rule local error indicator of adaptive cubature). NULL when
#'      local CCD was off or declined (single-block, `< 4` axes, an active
#'      `phi_grid`, or no peaked interior cell).
#'   * `integration_requested`, `integration_declined` -- what `integration`
#'      asked for, and why the requested CCD or adaptive lattice did not run. `$integration` names the
#'      integrator that RAN, and `.nl_node_support()` keys the interval
#'      construction off it, so a caller who asked for a moment rule and
#'      received a density grid reads the reason here rather than inferring it.
#'      `integration_declined` is `NA_character_` when nothing was declined, and
#'      otherwise one of `"axis_count"` (fewer transformable latent axes than the
#'      requested mode's threshold), `"unguessable_axis"` (a CAR_proper
#'      `rho_car` or a non-BYM2 `rho`), `"degenerate_axis"` (a single-valued
#'      axis), `"modefind_ridge"` / `"modefind_boundary"` /
#'      `"modefind_degenerate"` / `"modefind_failed"`, `"hessian_singular"`,
#'      `"hessian_not_pd"`, `"copy_atom_components"` (more copy atoms than the
#'      design splits into, see `copy_atom_mass`), `"copy_atom_mass"` (a
#'      declared atom mass outside `[0, 1)`) or `"placement_budget"` (placing
#'      the design would cost more inner solves than integrating the tensor grid
#'      it replaces; see the `ccd_budget` control knob). Under
#'      `integration = "grid_adaptive"` it is one of
#'      `"adaptive_small_grid"` (fewer outer cells than
#'      `adaptive_grid_min_cells`), `"adaptive_degenerate_lattice"`,
#'      `"adaptive_seed_failed"` (no finite marginal in the coarse seed),
#'      `"adaptive_diffuse_seed"` / `"adaptive_flood_cap"` /
#'      `"adaptive_kept_cap"` (the posterior mass spreads over more than
#'      `adaptive_grid_max_frac` of the dense grid, read at the seed, during the
#'      flood, or on the closed kept region), or `"adaptive_error"` (the flood
#'      raised an error, reported as a warning carrying its message); the fit
#'      then integrates the dense tensor. Multi-block fits only.
#'   * `ccd_modefind_evals`, `ccd_modefind_budget`, `ccd_modefind_rounds`,
#'      `ccd_modefind_seconds` -- what PLACING the CCD cost: inner Newton solves
#'      spent by the mode-find (seed, step calibration and every round of every
#'      mixture component), the ceiling they were spent against, the rounds
#'      taken and the wall-clock seconds. Present on any multi-block fit that
#'      attempted a placement, including one that declined -- a decline reports
#'      what it had already paid -- and absent when no placement was attempted,
#'      which is not the same as a placement that cost nothing. Read alongside
#'      `$integration` to judge whether the CCD earned its design on this model
#'      rather than arguing it from the round caps.
#'   * `weight_kind` -- one entry per outer-grid cell, `"mass"` or `"design"`.
#'      A tensor cell holds the mass of its own cell and a CCD node holds a
#'      design weight, so a fit integrated by one rule reports one value
#'      throughout; a locally refined grid carries both, and this says which per
#'      cell instead of leaving it to be read off `integration`.
#'   * `theta_interval_read`, `theta_interval_design_mass` -- what the reported
#'      per-axis `theta_median` / `theta_ci_lo` / `theta_ci_hi` were read off, and
#'      how much of the integration weight sits on nodes whose cumulative sum is
#'      not a CDF. `"density"` is a weighted quantile over cells that discretize
#'      the posterior (`design_mass` 0), `"moment_rule"` an interval from the
#'      moments a central-composite design delivers (`design_mass` 1), and
#'      `"mixed"` the locally CCD-refined grid, which carries both kinds at once.
#'      A mixed support still reports the weighted quantile, which measured best
#'      against a converged reference in both `design_mass` regimes; the share is
#'      the regime variable, since on that part of the support the quantile is
#'      bounded by the refined cells' own grid neighbourhoods rather than by the
#'      posterior. Every nested path stamps the pair, not only
#'      the multi-block driver.
#'   * `within_cell_requested`, `theta_within_cell`,
#'      `theta_within_cell_declined` -- the WITHIN-CELL construction the same
#'      intervals were read with. The kind above says what the
#'      integrator left; this says how each cell's mass was spread inside its own
#'      box when the grid was read back. `"log_quadratic"` is the default and
#'      interpolates each row's node log densities; `"box_uniform"`
#'      (`control$within_cell`) puts the cumulative full mass at each cell edge;
#'      `"chord"` puts the cumulative mid-mass at each cell coordinate, the same
#'      masses over the same boxes with the knots moved half a cell. The
#'      construction is recorded per axis, and an axis whose cell partition could
#'      not be built falls back to `"chord"` on its own with the reason in
#'      `theta_within_cell_declined` (`"support_<kind>"`, `"single_node"`,
#'      `"boxes_do_not_tile"`, `"no_usable_node"`). A `"log_quadratic"` request
#'      falls back to `"box_uniform"` first, with `"unresolved"` (no row
#'      resolves its own conditional), `"not_a_grid_axis"` or
#'      `"no_node_density"` as the reason.
#'   * `theta_cell_edge_coord`, `theta_cell_edge_declined` -- per axis, the
#'      COORDINATE that axis's outer cell edges were mirrored in, and why the
#'      support the axis declares did not produce them. A
#'      declared support is authoritative: its mirrored edge is used when finite
#'      and inside it, and otherwise the axis reports its extreme grid
#'      coordinate -- which the containment test has already placed inside the
#'      support -- rather than falling through to a coordinate guessed from the
#'      node values, which is what put a lower bound of 0 on a `positive` axis
#'      and a bound above 1 on a `unit` one. An axis that declares no support at
#'      all keeps the guess by design, but its mirrored edge is checked for
#'      being a representable double that brackets the coordinates, and falls
#'      back to the extreme grid coordinate with
#'      `"mirrored_edge_not_representable"` when it is not.
#'      `theta_cell_edge_declined` is `NA` where nothing was declined.
#'   * `outer_grid_cell_width`, `outer_grid_axis_sd`, `outer_grid_h_over_sd` --
#'      per axis, the cell width and the posterior SD in that axis's own
#'      coordinate, and their ratio. A within-cell reconstruction resolves an
#'      interval endpoint to within one cell, so how much of the reported width
#'      and of its realized coverage is a property of where the grid fell rather
#'      than of the posterior is governed by this ratio; below
#'      `.nl_diag("grid_resolved")` the cells are narrower than the posterior
#'      they discretize and the two constructions converge to each other. NA on
#'      an axis whose own marginal is maximal at an endpoint, which is the
#'      placement question `outer_grid_railed_axes` reports.
#'   * `outer_grid_resolution_declined` -- per axis, why that axis carries no
#'      ratio, from the closed vocabulary `too_few_nodes`, `mode_at_edge`,
#'      `coord_not_finite`, `spacing_not_finite`, `stencil_degenerate`,
#'      `curvature_not_negative`; `NA` on an axis that scored. The reasons are
#'      not interchangeable: `mode_at_edge` says the grid does not contain that
#'      axis's own posterior mode, which is a stronger statement about the fit
#'      than any ratio, and the whole-grid `resolved` verdict is withheld while
#'      one is present rather than being read off the axes that did score.
#'   * `prune_cheap_log_marginal`, `prune_mask`, `prune_n_pruned`,
#'      `prune_tol`, `prune_screen_iters` -- present only when `prune = TRUE`
#'      and the safety gate did not fall back. Cheap-pass log-marginals at every
#'      cell -- the truncated solve's value plus half its Newton decrement, the
#'      second-order estimate of the converged value the screen ranks by, so a
#'      cell whose inner mode had far to move from its warm start is not ranked
#'      on where the truncation left it (gcol33/tulpa#919); that decrement is
#'      `prune_screen_decrement` --
#'      a logical mask of pruned cells, the pruned-cell count, the
#'      threshold actually applied, and the screening depth the driver actually
#'      ran at. Pruned cells have `log_marginal = -Inf` so they get zero
#'      weight under `.nl_normalise_weights_safe`.
#'   * `prune_log_gap_cut`, `prune_cheap_lm_spread`, `prune_min_keep`,
#'      `prune_n_floor_restored` -- the screen read on the scale it operates on:
#'      how many nats below the best screened cell the tolerance cut at, how
#'      many nats the screened surface spans, the kept-cell floor applied, and
#'      how many cells that floor put back. A cut that is a sliver of the spread
#'      is a tolerance with no resolution on this grid.
#'   * `prune_dropped_mass_bound`, `prune_n_repaired`, `prune_gate_mass` -- the
#'      posterior mass the dropped cells may carry, with each lifted by the
#'      worst screening error measured on a kept cell; how many dropped cells
#'      the driver solved in full to bring that under the gate; and the gate.
#'   * `outer_mode_carried` -- `TRUE` on a placement refit that read the outer
#'      mode it was laid from rather than finding it again, which it does unless
#'      the placed grid rails.
#'   * `prune_fallback_triggered`, `prune_fallback_reason` -- present only when
#'      a fallback to the full grid fired. The returned fit is the full-grid
#'      (unpruned) result; the reason string records whether the cheap-pass
#'      safety gate tripped or the screened grid left the placement pass no
#'      curvature.
#'   * `n_threads_outer_requested`, `n_threads_outer_realised` -- the outer-grid
#'      width asked for and the width the grid ran at. The driver clamps the
#'      request to the team the OpenMP environment hands out, so a fit launched
#'      under `OMP_NUM_THREADS=1` with `n_threads_outer = 10` runs serial; the
#'      pair is what a timing has to be recorded against.
#'   * `converged`, `n_iter`, `score_max` -- per cell, whether the inner Newton
#'      reached a mode within `control$max_iter`, the iterations it took, and
#'      the largest absolute component of the joint penalized score at the
#'      mode it reported (the residual a cell that stopped short stopped at). A cell that
#'      stopped without a mode is solved again from the mode of its nearest
#'      converged cell (`nonconverged_resolve` counts the cells that first
#'      stalled and those the re-solve recovered); one that still has no mode
#'      stays in the integration at the value it stopped at, is listed in
#'      `nonconverged_cells` with the posterior mass such cells carry in
#'      `nonconverged_mass`, and is named by [diagnostics()]. Past the screen's
#'      mass gate the fit warns.
#'
#' @seealso [tulpa_nested_laplace()] for the single-arm engine.
#' @references
#' Rue, Martino & Chopin (2009). Approximate Bayesian inference for latent
#' Gaussian models by using integrated nested Laplace approximations.
#' \emph{JRSS-B} 71(2):319-392.
#' @examples
#' \donttest{
#' set.seed(1)
#' S <- 25L                                   # shared spatial units (chain graph)
#' nb <- lapply(seq_len(S), function(s) setdiff(c(s - 1L, s + 1L), c(0L, S + 1L)))
#' nn <- lengths(nb)
#' field <- as.numeric(scale(cumsum(rnorm(S, 0, 0.4))))
#' mk_arm <- function(m, fam) {
#'   si <- sample(S, m, replace = TRUE); x <- rnorm(m)
#'   lin <- 0.2 + 0.5 * x + 0.8 * field[si]
#'   y <- if (fam == "binomial") rbinom(m, 1L, plogis(lin)) else lin + rnorm(m, 0, 0.5)
#'   list(y = as.numeric(y), n_trials = rep(1L, m), X = cbind(1, x),
#'        spatial_idx = si, family = fam, phi = if (fam == "gaussian") 0.5 else 1)
#' }
#' prior <- list(type = "icar", n_spatial_units = S,
#'               adj_row_ptr = c(0L, cumsum(nn)), adj_col_idx = unlist(nb) - 1L,
#'               n_neighbors = nn, sigma_grid = c(0.3, 0.7, 1.4))
#' fit <- tulpa_nested_laplace_joint(
#'   responses = list(occ = mk_arm(200L, "binomial"), pos = mk_arm(200L, "gaussian")),
#'   prior = prior)
#' fit$theta_mean        # shared field amplitude, integrated across both arms
#' }
#' @export
tulpa_nested_laplace_joint <- function(responses,
                                       prior,
                                       copy = NULL,
                                       phi_grid = NULL,
                                       prior_sigma = NULL,
                                       prior_alpha = NULL,
                                       prior_phi = NULL,
                                       hyperprior = c("proper", "flat"),
                                       cell_coupling = "separable",
                                       control = list()) {
    # Renamed knob first (a targeted message beats the generic whitelist one),
    # then the whitelist, both before any response validation so a bad knob is
    # reported regardless of the rest of the call.
    if (!is.null(control$diagnose_draws)) {
        stop("`control$diagnose_draws` was renamed: use `control$k_samples` ",
             "(the same knob on every nested-Laplace fitter).", call. = FALSE)
    }
    tulpa_check_control(control, .CONTROL_KEYS$nested_laplace_joint,
                   "tulpa_nested_laplace_joint")
    hyperprior <- .hp_choice(match.arg(hyperprior))
    # k_quality reliability front door + escalation. The
    # single fit lives in .tulpa_nl_joint_once(); this wrapper resolves the
    # reliability intent, attaches the honest verdict (for BOTH the single- and
    # multi-block paths), and -- for "ok" / "good" -- chases the band by REFINING
    # THE INTEGRATION GRID, driven by the bad outer k (with k_refine = "grid"):
    # each round widens / densifies the grid where the posterior mass escapes its
    # bounds, re-fits + re-diagnoses, until the requested band is confidently
    # reached or the round budget is hit. k_samples is the separate
    # estimate-precision knob, not an escalation lever. Never silently downgrades.
    k_quality <- control$k_quality %||% "report"
    if (!is.character(k_quality) || length(k_quality) != 1L ||
        !k_quality %in% c("none", "report", "ok", "good")) {
        stop("`control$k_quality` must be one of \"none\", \"report\", \"ok\", \"good\".",
             call. = FALSE)
    }
    k_refine <- control$k_refine %||% "grid"
    if (!is.character(k_refine) || length(k_refine) != 1L ||
        !k_refine %in% c("none", "grid", "ccd")) {
        stop("`control$k_refine` must be \"none\", \"grid\", or \"ccd\".", call. = FALSE)
    }
    k_max_rounds <- control$k_max_rounds %||% 2L
    if (length(k_max_rounds) != 1L || is.na(k_max_rounds) ||
        k_max_rounds != round(k_max_rounds) || k_max_rounds < 0L) {
        stop("`control$k_max_rounds` must be a single integer >= 0.", call. = FALSE)
    }
    k_max_rounds <- as.integer(k_max_rounds)
    k_conf_bands <- control$k_conf_bands
    diag_raw     <- control$diagnose_k %||% TRUE
    diagnose_k   <- (identical(diag_raw, "by_arm") || isTRUE(diag_raw)) &&
                    !identical(k_quality, "none")

    attach_q <- function(res)
        .joint_attach_k_quality(res, k_quality, diagnose_k,
                                res$diagnose_draws %||% NA_integer_, k_conf_bands)

    # Outer-axis provenance, taken ONCE before the first fit: which grid axes
    # the caller declared as defaults with `auto_grid()` (as opposed to pinned
    # deliberately), and the same prior with the markers stripped so nothing
    # downstream sees an attributed numeric. The auto-recenter passes below read
    # the record; everything else reads plain grids.
    prov  <- .nl_grid_provenance(prior)
    prior <- prov$prior
    # The same question for the per-arm dispersion axes, which live on
    # `phi_grid` rather than on a prior block (gcol33/tulpa#663). Taken BEFORE
    # `.normalise_phi_grid()`, which coerces each entry with `as.numeric()` and
    # would drop the marker; normalising here rather than leaving it to
    # `.tulpa_nl_joint_once()` (which does it again, idempotently) is what gives
    # the rescue below a stable `phi_grid[[arm]]` write target whichever spelling
    # the caller used.
    arm_names <- names(responses) %||% paste0("arm", seq_along(responses))
    phi_prov  <- .nl_phi_provenance(phi_grid, arm_names)
    phi_grid  <- .normalise_phi_grid(phi_prov$phi_grid, arm_names)
    # Refuse an axis the resolved per-block path cannot read, and record a
    # dropped default.
    .op_axis <- .nl_publish_axis_dropped(
        .nl_check_axis_fields(prior, "joint", auto = prov$auto,
                              copy = copy, responses = responses))
    on.exit(options(.op_axis), add = TRUE)
    # `control$auto_recenter = FALSE` holds every grid exactly as given -- the
    # opt-out for a caller who wants the default axis integrated as-is. The
    # joint rescues trigger on the whole grid's `collapsed_edge` regime rather
    # than on a per-axis rail, so the per-axis policy names (`"rail"`,
    # `"resolve"`, `"always"`) have no measured meaning here and are refused
    # rather than accepted and ignored.
    if (is.character(control$auto_recenter)) {
        stop("control$auto_recenter = \"",
             paste(control$auto_recenter, collapse = "\", \""),
             "\" names a per-axis placement policy implemented on the ",
             "standalone tulpa_nested_laplace() registry path only; the joint ",
             "path recentres on the grid's collapsed-edge regime, not a ",
             "per-axis rail. Use TRUE or FALSE here.", call. = FALSE)
    }
    auto_recenter <- !isFALSE(control$auto_recenter)

    # Outer-grid width, resolved ONCE for the fit: the driver clamps the request
    # to the team the OpenMP environment will hand out, and every kernel call
    # below goes through that clamp, so reporting it here says it once rather
    # than once per grid solve.
    outer_threads <- .nl_outer_threads(control$n_threads_outer %||% 1L,
                                       fn = "tulpa_nested_laplace_joint()")

    # The call state the placement pass reads and rewrites
    # (R/joint_placement_pass.R): every argument a placed axis lives on, and the
    # resolved copy spec a copy block's axes are found through.
    state <- list(prior = prior, prior_sigma = prior_sigma, phi_grid = phi_grid,
                  responses = responses, copy = copy,
                  cp = tryCatch(.resolve_copy_multi(copy, responses, prior),
                                error = function(e) NULL))
    families <- .joint_placement_families(auto = prov$auto,
                                          phi_auto = phi_prov$auto)

    # The axes beyond a field SD the placement pass may move -- the per-arm
    # dispersions and the copy scale -- named as the grid names them. Handed to
    # every fit so the mode-find runs on a grid that concentrated without
    # railing; empty when every one is held, which is the gate that keeps such a
    # fit paying nothing for a pass that could not run.
    placement_axes <- if (!auto_recenter) character(0) else
        .joint_movable_extra_axes(state, phi_prov$auto, prov$auto)

    # Grid-cell checkpoint. `resume = FALSE` starts THIS FIT over, so any prior
    # file is removed once, here, before the first solve. An outer-grid
    # placement refits, so the per-solve `.tulpa_nl_joint_once()` runs several
    # times within one fit; taking the removal there deleted the cells an
    # earlier solve of the same fit had already written, and a later resume
    # then re-appended the whole pre-placement grid on top of the survivors.
    .ckpt_fit <- .nl_checkpoint_args(control)
    if (nzchar(.ckpt_fit$path) && !isTRUE(.ckpt_fit$resume) &&
        file.exists(.ckpt_fit$path)) {
        file.remove(.ckpt_fit$path)
    }

    ctrl <- control
    fit_once <- function(prior_i, prior_sigma_i, ctrl_i = ctrl,
                         phi_grid_i = phi_grid, copy_i = copy,
                         responses_i = responses)
        attach_q(.tulpa_nl_joint_once(responses_i, prior_i, copy_i, phi_grid_i,
                                      prior_sigma_i, prior_alpha, prior_phi,
                                      cell_coupling, ctrl_i,
                                      placement_axes = placement_axes,
                                      hyperprior = hyperprior,
                                      phi_auto = phi_prov$auto,
                                      grid_auto = prov$auto))
    fit_state <- function(st, ctrl_i = ctrl)
        fit_once(st$prior, st$prior_sigma, ctrl_i = ctrl_i,
                 phi_grid_i = st$phi_grid, copy_i = st$copy,
                 responses_i = st$responses)

    # Placement pilot (gcol33/tulpa#636). Placement reads an argmax cell and an
    # FD curvature stencil, and reads them off `log_marginal` -- not off the
    # integration the detecting pass paid for, every cell of which is discarded
    # the moment a placement fires. With the pilot on, the DETECTING grid is a
    # thinned read of the same spans and the full grid is solved once, at the
    # placed axes. The placement pass is untouched by this: it detects on the
    # fit it is handed and writes onto the STATE it is handed, so passing the
    # pilot fit alongside the full state places the full grid. What the pilot
    # obliges this caller to add is the other half -- a pilot grid is too coarse
    # to integrate, so a placement that declines is followed by the full fit the
    # pilot stood in for (`pilot_fallback` below).
    pilot <- .nl_recenter_pilot(
        prior, phi_grid, copy, responses,
        copy_blocks = as.integer(state$cp$copy_blocks_zero %||% integer(0)) + 1L,
        n = .nl_pilot_n(control), enabled = auto_recenter)

    res <- if (isTRUE(pilot$active)) {
        fit_once(pilot$prior, prior_sigma, .nl_pilot_control(ctrl),
                 pilot$phi_grid, pilot$copy, pilot$responses)
    } else {
        fit_once(prior, prior_sigma)
    }
    # What the pilot DETECTED on, kept before its fit is replaced: the trigger
    # is read off a grid nothing downstream will hold, so without this a fit
    # placed from a pilot cannot say what the placement was decided from.
    pilot_detect <- if (!isTRUE(pilot$active)) NULL else list(
        cells      = nrow(res$theta_grid),
        regime     = res$pareto_k_regime,
        ess_grid   = .nl_grid_ess(res$weights),
        edge_axes  = res$pareto_k_grid_edge_axes,
        screened   = isTRUE(res$prune_screen_only))
    res$k_quality_rounds <- 0L
    res$outer_grid_placement <- res$outer_grid_placement %||% "fixed"

    # Place the grid before the k_quality escalation below: that loop assumes
    # the grid's BOUNDS are right and only needs densifying / refining, so a
    # mis-centred axis (mode beyond the fixed ceiling) must be fixed first, not
    # densified. A no-op (zero extra fit) unless an axis fires; the call state
    # is reassigned to the placed axes so escalation, if it still runs,
    # continues from there.
    #
    # A refit is handed the outer mode the placement was read off, so the
    # var-of-means consistency pass inside it lays an axis still collapsed onto
    # one node at that mode instead of bisecting towards it
    # (`.hyper_propose_at_mode()`).
    place <- function(res_i, st_i, ctrl_i = ctrl)
        .joint_place_axes(
            res_i, st_i,
            refit = function(st, from) {
                op <- options(tulpa.nl_outer_mode = .nl_fit_outer_mode(from))
                on.exit(options(op), add = TRUE)
                fit_state(st, ctrl_i)
            },
            families = families, enabled = auto_recenter)
    adopt <- function(placed) {
        state       <<- placed$st
        prior       <<- state$prior
        prior_sigma <<- state$prior_sigma
        phi_grid    <<- state$phi_grid
        responses   <<- state$responses
        copy        <<- state$copy
        placed$res
    }

    res <- adopt(place(res, state))

    # The pilot's other half. A placement that fired has already produced a full
    # fit at the placed axes; one that DECLINED has left `res` on the pilot's own
    # coarse grid, which nothing may integrate -- the fit the caller asked for is
    # the full grid the pilot stood in for. That full fit then gets its OWN
    # detection, which is what keeps the pilot a pure PRE-screen: a coarse grid
    # can only ADD a placement, never cost one, so a fit is never placed less
    # well than the same fit with the knob off. When the full grid declines too
    # (the common case), the second detection reads stored weights and costs
    # nothing, and the fit is bit-identical to the un-piloted one.
    if (isTRUE(pilot$active)) {
        if (!identical(res$outer_grid_placement, "auto_recentered")) {
            res <- fit_state(state)
            res$k_quality_rounds     <- 0L
            res$outer_grid_placement <- res$outer_grid_placement %||% "fixed"
            res <- adopt(place(res, state))
            pilot_placed <- FALSE
        } else {
            pilot_placed <- TRUE
        }
        res <- .nl_pilot_attach(res, pilot, pilot_detect,
                                declined = if (pilot_placed) NULL else
                                    "pilot_placement_declined")
    }

    # The cheap-pass screen's other failure, and the one its own safety gate
    # cannot see: the kept set carries the posterior but no curvature, so the
    # placement pass above declined and the outer axis stayed where it was laid.
    # A hyperparameter read off such a fit is an axis endpoint, not an estimate,
    # so the screened attempt is replaced by the full grid exactly as a
    # distrusted ranking is. The kept-set floor makes this rare; it is the net
    # under it. Everything after this point runs with screening off, so the
    # k_quality escalation below does not re-screen the grid it is refining.
    if (.nl_prune_placement_lost(res)) {
        ctrl_full <- utils::modifyList(ctrl, list(prune = FALSE))
        res <- .nl_prune_placement_fallback(res, function() {
            r <- fit_state(state, ctrl_full)
            r$k_quality_rounds     <- 0L
            r$outer_grid_placement <- r$outer_grid_placement %||% "fixed"
            adopt(place(r, state, ctrl_full))
        })
        ctrl <- ctrl_full
    }

    res$k_quality_rounds     <- res$k_quality_rounds %||% 0L
    res$outer_grid_placement <- res$outer_grid_placement %||% "fixed"
    res$k_quality_k_trace    <- res$pareto_k

    if (k_quality %in% c("ok", "good") && isTRUE(diagnose_k) &&
        k_max_rounds > 0L && k_refine %in% c("grid", "ccd")) {
        # A miss has two causes, and the ladder is ordered by which lever
        # reduces which quantity -- not by which cause the verdict named.
        #
        # Refining the grid lowers the k-hat ITSELF: the proposal is rebuilt
        # from a grid that now covers the posterior better. Buying importance
        # draws lowers only the k-hat ESTIMATOR's variance -- it narrows the
        # bootstrap CI around whatever the k-hat already is. Refinement
        # therefore goes FIRST on either miss: on a "resolution" miss it is the
        # only lever, and on a "precision" miss (the CI straddling a band
        # boundary) it can move the underlying k INTO the band rather than
        # merely resolving whether it was already there. Measured on the
        # fixture in `test-nested-laplace-joint-pareto-k.R`, which reaches the
        # good band from an "uncertain" first fit by refining and does not by
        # sharpening.
        #
        # The two refinement rungs, re-diagnosed each round:
        #   * "grid" -- the grid's bounds let posterior mass escape: each round
        #     runs the adaptive boundary-extension / interior-densification pass
        #     (R/hyper_grid_refine.R) where the integrand mass piles at an edge,
        #     one more pass per round (recursive in the bad-k direction).
        #   * "ccd"  -- the grid is too coarse to resolve a sharply-peaked
        #     direction: each round refines more high-weight, mutually
        #     non-adjacent cells with local curvature-aware CCD node clouds
        #     (R/nested_laplace_joint_ccd_local.R). A tensor base is forced (the
        #     finite-difference curvature stencil needs axis neighbours).
        # When refinement stops producing -- no boundary mass left to chase, no
        # peaked interior cell to split -- the rung is EXHAUSTED, and that is
        # where the two causes finally separate. On a "resolution" miss nothing
        # remains: the k-hat is confidently outside the band and the grid cannot
        # be improved, so the loop stops and says so. On a "precision" miss the
        # estimate is still ambiguous and there IS a lever left, the one
        # refinement never touched: fit the GPD tail on more actual tail ratios.
        # Those rounds buy draws (`k_precision_growth` x the realised budget)
        # and hold the grid where the last refinement left it. Stopping there
        # instead -- which is what a loop reading only "reached / not reached"
        # does -- abandons a fit whose integration may already meet the request,
        # for want of the cheaper half of the diagnostic.
        #
        # The refinement level advances on refinement rounds only, counted by
        # `passes` rather than by `round`, so a draw round does not silently
        # skip a pass. `k_quality_k_trace` records the k-hat the chase actually
        # walked, round by round: the refinement rungs rebuild the proposal from
        # the grid they just changed, so a rising trace is a re-founded proposal
        # on a strictly larger grid rather than a worse fit, and the trace is
        # what makes that visible instead of leaving only the last round's
        # number.
        is_ccd_refine <- identical(k_refine, "ccd")
        refined_field <- if (is_ccd_refine) "local_ccd_info" else "adaptive_grid_info"
        lc_base      <- if (is.list(control$local_ccd)) control$local_ccd else list()
        lc_max_cells <- as.integer(lc_base$max_cells %||% 6L)
        round     <- 0L
        passes    <- 0L                            # refinement level so far
        n_draws_r <- 0L                            # rounds spent on tail ratios
        exhausted <- FALSE                         # refinement has nothing left
        k_trace   <- res$pareto_k
        while (!isTRUE(res$k_quality_reached) && round < k_max_rounds) {
            round <- round + 1L
            draws_round <- exhausted
            if (draws_round) {
                # Buy tail ratios and leave `ctrl`'s integration settings
                # exactly where the last refinement round left them, so the
                # round holds the grid rather than advancing it by the default
                # one pass. The refinement switches are turned on by the branch
                # below, never here, and a caller who asked for `adaptive_grid`
                # themselves keeps it.
                #
                # The GPD tail size is PINNED to the fraction the fit's own
                # first pass used. The automatic PSIS rule is
                # `min(S/5, 3 sqrt(S))`, so the fitted tail FRACTION shrinks as
                # `3 / sqrt(S)`: left to itself, buying draws would fit the
                # shape at a deeper quantile of the weight distribution rather
                # than fit the same one more precisely -- it would move the
                # estimand, which is the one thing this rung must not do
                # (gcol33/tulpa#631). Measured on synthetic outer targets, the
                # automatic rule takes a k-hat from 0.57 to 7.94 over a 100x
                # budget while a held fraction holds it at 0.57 to 0.69 with the
                # seed spread narrowing. The fraction is at most 1/5 by
                # construction, so this never trips the 20% cap warning.
                n_draws_r <- n_draws_r + 1L
                k_prev <- as.integer(res$diagnose_draws %||% .nl_diag("k_samples"))
                k_next <- as.integer(ceiling(
                    .nl_diag("k_precision_growth") * k_prev))
                tp_prev <- res$pareto_k_tail_points
                if (is.null(ctrl$k_tail_points) && length(tp_prev) == 1L &&
                    is.finite(tp_prev) && k_prev > 0L) {
                    ctrl$k_tail_points <- max(
                        .PSIS_MIN_EVAL,
                        as.integer(floor(k_next * (tp_prev / k_prev))))
                }
                ctrl$k_samples <- k_next
            } else {
                passes <- passes + 1L
                if (is_ccd_refine) {
                    ctrl$integration <- "grid"
                    ctrl$local_ccd   <- utils::modifyList(
                        lc_base, list(max_cells = passes * lc_max_cells))
                } else {
                    ctrl$adaptive_grid            <- TRUE
                    ctrl$adaptive_grid_max_passes <- passes
                }
            }
            res <- attach_q(.tulpa_nl_joint_once(responses, prior, copy, phi_grid,
                                                 prior_sigma, prior_alpha, prior_phi,
                                                 cell_coupling, ctrl,
                                                 hyperprior = hyperprior))
            res$k_quality_rounds <- round
            k_trace <- c(k_trace, res$pareto_k)
            if (!isTRUE(res$k_quality_reached) && !draws_round &&
                is.null(res[[refined_field]])) {
                # The refinement produced no info object, so no further round
                # can add a node. For the CCD rung this means no peaked interior
                # cell was found. For the grid rung it means the grid was not
                # extended this round -- either no boundary mass beyond the
                # current width, or (on a multi-block fit) the outer integrator
                # is a CCD design or the adaptive lattice, whose nodes the
                # boundary pass does not extend.
                exhausted <- TRUE
                if (identical(res$k_quality_miss, "precision")) {
                    # What is left is the interval's width, not the grid, and
                    # that lever is untouched. Keep going on tail ratios.
                    next
                }
                # A confidently-bad k-hat on a grid that cannot be refined. The
                # message states what happened without asserting the bad k is
                # definitely not a grid deficiency.
                res$k_quality_reason <- if (is_ccd_refine) paste0(
                    "local CCD found no peaked interior cell to refine (needs a ",
                    ">= 4-axis multi-block tensor grid); the bad k is not a ",
                    "coarse-grid resolution deficiency") else paste0(
                    "grid refinement did not extend the outer grid this round ",
                    "(no boundary mass beyond the current width, or the outer ",
                    "integrator is a CCD design or an adaptive lattice, whose ",
                    "nodes are not extended)")
                break
            }
        }
        if (!isTRUE(res$k_quality_reached) && round >= k_max_rounds) {
            # Name what the budget was actually spent on, so a caller raising
            # k_max_rounds knows which lever the extra rounds will pull.
            rung <- if (is_ccd_refine) "local-CCD-refinement" else "grid-refinement"
            res$k_quality_reason <- sprintf(
                "requested band not confirmed within the k_max_rounds budget (%s)",
                if (n_draws_r == 0L) sprintf("%d %s", passes, rung)
                else if (passes == 0L) sprintf("%d draw-budget", n_draws_r)
                else sprintf("%d %s + %d draw-budget", passes, rung, n_draws_r))
        }
        res$k_quality_k_trace <- k_trace
    }
    res$outer_grid_single_node_axes <- .nl_single_node_axes(res)
    .nl_warn_unplaced_rail(res, "tulpa_nested_laplace_joint()")
    .nl_warn_single_node_axes(res, "tulpa_nested_laplace_joint()")
    .nl_warn_unconverged(res, "tulpa_nested_laplace_joint()")
    .nl_warn_indefinite(res, "tulpa_nested_laplace_joint()")
    .nl_attach_outer_threads(res, outer_threads)
}

.tulpa_nl_joint_once <- function(responses, prior, copy = NULL, phi_grid = NULL,
                                 prior_sigma = NULL, prior_alpha = NULL,
                                 prior_phi = NULL,
                                 cell_coupling = "separable", control = list(),
                                 placement_axes = character(0),
                                 hyperprior = "proper",
                                 phi_auto = NULL,
                                 grid_auto = logical(0)) {
    tm <- .tulpa_timer()
    # Resolve and validate the cell-coupling spec name against the C++
    # registry (separable default is auto-registered on first touch). The
    # resolved spec is held on the result but not yet routed through the
    # inner-Newton scatter -- the per-cell branch in
    # scatter_arm_obs_joint_multi[_sparse] lands with Layer B of Change 2.
    # In Layer A the default sentinel ("separable", arm_ids() empty) keeps
    # every existing joint fit on the per-obs path with no behavioural change.
    if (!is.character(cell_coupling) || length(cell_coupling) != 1L ||
        is.na(cell_coupling) || !nzchar(cell_coupling)) {
        stop("`cell_coupling` must be a single non-empty character string ",
             "naming a registered spec (default: \"separable\").",
             call. = FALSE)
    }
    if (!cpp_cell_coupling_registry_has(cell_coupling)) {
        stop("`cell_coupling = \"", cell_coupling, "\"` is not registered. ",
             "Consumer packages register specs from R_init_<pkg> via the ",
             "tulpa_register_cell_coupling C callable; the built-in ",
             "\"separable\" default is always available.",
             call. = FALSE)
    }

    # Perf/numerical knobs live in `control = list()` (matching tulpa()); the
    # top-level signature carries only statistical arguments. Names were
    # validated at the front of tulpa_nested_laplace_joint().
    max_iter                  <- control$max_iter %||% 50L
    tol                       <- control$tol %||% 1e-6
    n_threads_outer           <- control$n_threads_outer %||% 1L
    # The kernel's inner thread count is the fit's whole grant. A solve inside
    # the outer team takes its share of the pool whatever this says; a solve
    # that runs alone -- a grid's pilot, a refinement or consistency round, the
    # skew probe -- takes all of it (`nl_cell_solve_threads()`).
    n_threads                 <- .tulpa_inner_threads(
        .tulpa_thread_grant(n_threads_outer, control$n_threads %||% 1L),
        control$n_threads_scatter)
    tile_warm                 <- control$tile_warm %||% TRUE
    # Local CCD refinement reads a finite-difference stencil off each refined
    # cell's axis neighbours on the solved grid, and a dropped neighbour carries
    # no log-marginal to read, so a fit that asks for it solves every cell
    # unless its caller asked for the screen outright.
    prune                     <- control$prune %||%
        (.nl_screen("prune")[["joint"]] && is.null(control$local_ccd))
    # `prune_log_gap` states the same screening cut in nats and replaces the
    # tolerance when set; the two together are refused rather than ranked.
    prune_tol                 <- .nl_prune_tol_from_control(
        control, control$prune_tol %||% .nl_screen("prune_tol"))
    screen_iters              <- control$screen_iters %||% .nl_screen("iters")
    x_init                    <- control$x_init
    verbose                   <- control$verbose %||% FALSE
    store_Q                   <- control$store_Q %||% FALSE
    keep_grid_hessians        <- isTRUE(control$keep_grid_hessians %||% TRUE)
    # The subspace-debias coupling closure reads the modal cell's joint
    # precision, which only `store_Q` produces. Requesting the closure
    # therefore requests the precision, the same way requesting the debias
    # already turns `keep_grid_hessians` on (gcol33/tulpa#862).
    keep_joint_prec <- {
      sd_cfg <- .subspace_debias_config(control$subspace_debias)
      !is.null(sd_cfg) && !identical(sd_cfg$closure, FALSE)
    }
    if (keep_joint_prec) store_Q <- TRUE
    adaptive_grid             <- control$adaptive_grid %||% FALSE
    adaptive_grid_edge_thresh <- control$adaptive_grid_edge_thresh %||% 0.02
    adaptive_grid_max_passes  <- control$adaptive_grid_max_passes %||% 1L
    axis_refine               <- control$axis_refine
    var_of_means_consistency  <- control$var_of_means_consistency %||% TRUE
    copy_atom_mass            <- control$copy_atom_mass %||% .TULPA_COPY_ATOM_MASS
    copy_slab                 <- .hyper_check_copy_slab(control$copy_slab)
    force_sparse              <- control$force_sparse %||% FALSE
    inner_sparse_override     <- .resolve_inner_factorization(
        control$inner_factorization)
    # Local CCD refinement of the multi-block outer grid. A
    # coarse tensor base grid is refined by a small curvature-aware node cloud on
    # a few high-weight, mutually non-adjacent cells, so the base grid can stay
    # coarse at moderate-to-high latent dimension where a uniformly fine tensor is
    # k^d-expensive. `control$local_ccd` is TRUE (defaults) or a list(max_cells=,
    # f0=, skew_max=, rank=); NULL (the default) keeps the unrefined grid. Engages only on
    # the multi-block tensor path at >= 4 transformable latent axes with no active
    # phi grid; the single-block path always integrates on the tensor grid
    # unrefined.
    local_ccd                 <- control$local_ccd
    if (!is.null(local_ccd) && !isTRUE(local_ccd) && !is.list(local_ccd)) {
        stop("`control$local_ccd` must be NULL, TRUE, or a ",
             "list(max_cells=, f0=, skew_max=, rank=).", call. = FALSE)
    }
    # Adaptive-lattice integrator tuning (integration = "grid_adaptive", the
    # low-dimensional multi-block companion to the CCD; see
    # nested_laplace_joint_adaptive.R). `adaptive_grid_cutoff` is the log-density
    # keep/expand radius from the peak (larger = closer to the dense tensor, fewer
    # cells skipped); `adaptive_grid_stride` the coarse-seed subsample stride;
    # `adaptive_grid_max_frac` the fraction of the dense grid past which the
    # builder declines back to the tensor.
    adaptive_cutoff           <- control$adaptive_grid_cutoff   %||% 10
    adaptive_stride           <- as.integer(control$adaptive_grid_stride %||% 2L)
    adaptive_max_frac         <- control$adaptive_grid_max_frac %||% 0.75
    adaptive_min_cells        <- control$adaptive_grid_min_cells %||% 48
    # Outer Pareto-k-hat accuracy diagnostic. `diagnose_k` (default TRUE)
    # importance-samples the joint hyperparameter posterior against the proposal
    # the integrator fits; `k_samples` (default 500) is the number of draws,
    # each one extra inner joint solve. RNG-restored, so the fit's draws are
    # unchanged whether or not it runs. `diagnose_k = "by_arm"`
    # additionally computes a k-hat restricted to each arm's hyperparameter axes
    # (other arms held at their posterior mean), to localise which arm drives a
    # tail-heavy joint k; the joint k is unchanged and stays
    # the default.
    # k_quality within ONE fit: "none" disables the diagnostic,
    # and "ok" / "good" raise the default draw budget so the bootstrap CI can
    # resolve the requested band. The verdict quartet, the band classification, and
    # the multi-round escalation live in the front-door wrapper
    # tulpa_nested_laplace_joint(), which drives this single-fit engine.
    k_quality                 <- control$k_quality %||% "report"
    if (!is.character(k_quality) || length(k_quality) != 1L ||
        !k_quality %in% c("none", "report", "ok", "good")) {
        stop("`control$k_quality` must be one of \"none\", \"report\", \"ok\", \"good\".",
             call. = FALSE)
    }
    diagnose_k_raw            <- control$diagnose_k %||% TRUE
    pareto_k_by_arm           <- identical(diagnose_k_raw, "by_arm")
    diagnose_k                <- (pareto_k_by_arm || isTRUE(diagnose_k_raw)) &&
                                 !identical(k_quality, "none")
    # Diagnostic importance-draw budget. `k_samples` is the
    # single precision knob -- the SAME name every nested-Laplace fitter uses:
    # the outer Pareto-k is scored ONCE over this many importance draws, and a
    # tighter k needs MORE actual tail ratios, i.e. a larger `k_samples`.
    diagnose_draws_user       <- control$k_samples
    diagnose_draws            <- diagnose_draws_user %||% .nl_diag("k_samples")
    if (length(diagnose_draws) != 1L || is.na(diagnose_draws) ||
        diagnose_draws != round(diagnose_draws) || diagnose_draws < 1L) {
        stop("`control$k_samples` must be a single integer >= 1.", call. = FALSE)
    }
    diagnose_draws            <- as.integer(diagnose_draws)
    # k_quality "ok" / "good" raise the default draw budget so the bootstrap CI can
    # resolve the requested band (a tighter CI needs more actual tail ratios); an
    # explicit `k_samples` is always respected.
    if (is.null(diagnose_draws_user) && k_quality %in% c("ok", "good"))
        diagnose_draws <- as.integer(.nl_diag(
            if (identical(k_quality, "good")) "k_samples_good" else "k_samples_ok"))
    # Outer-thread width for the diagnostic's importance batch.
    # The `diagnose_draws` re-solves are independent and run after the grid (all
    # cores free), each solved single-threaded, so widening this pool is a
    # bit-identical wall-clock speedup. `NULL` (default) follows the fit's thread
    # grant; "auto" uses the performance cores; an integer pins the width.
    k_threads                 <- control$k_threads
    pareto_k_threads          <- .tulpa_pareto_k_threads(n_threads_outer,
                                                         n_threads, diagnose_draws,
                                                         k_threads)
    # Bootstrap outer Pareto-k uncertainty. The chosen proposal
    # is scored once; the k-hat's sampling uncertainty is then estimated by
    # bootstrapping its raw importance log-ratios -- `k_bootstrap` replicates, each
    # re-fitting the GPD tail, NO new inner solves. The bootstrap reports how
    # UNSTABLE the current estimate is, not a tighter k (raise `k_samples` for
    # that, NOT `k_bootstrap`). `k_tail_points` (NULL = the automatic PSIS rule) is
    # an expert tail-threshold control, capped at the 20%-of-draws ceiling.
    # `k_conf_bands` are the reliability-band boundaries the bootstrap CI is tested
    # against for `pareto_k_band_confident`; NULL (default) uses the
    # sample-size-dependent boundaries c(0.5, min(1 - 1/log10(S), 0.7)) at the
    # realised draw count S.
    k_bootstrap               <- control$k_bootstrap %||% .nl_diag("k_bootstrap")
    if (length(k_bootstrap) != 1L || is.na(k_bootstrap) ||
        k_bootstrap != round(k_bootstrap) || k_bootstrap < 0L) {
        stop("`control$k_bootstrap` must be a single integer >= 0.", call. = FALSE)
    }
    k_bootstrap               <- as.integer(k_bootstrap)
    k_tail_points             <- control$k_tail_points
    if (!is.null(k_tail_points)) {
        k_tail_points <- suppressWarnings(as.integer(k_tail_points))
        if (length(k_tail_points) != 1L || is.na(k_tail_points) || k_tail_points < 1L) {
            stop("`control$k_tail_points` must be a single positive integer, or NULL.",
                 call. = FALSE)
        }
    }
    k_conf_bands              <- control$k_conf_bands
    if (!is.null(k_conf_bands) &&
        (!is.numeric(k_conf_bands) || anyNA(k_conf_bands) || length(k_conf_bands) < 1L ||
         is.unsorted(k_conf_bands, strictly = TRUE))) {
        stop("`control$k_conf_bands` must be NULL (the sample-size-dependent ",
             "default) or a strictly increasing numeric vector.", call. = FALSE)
    }
    # Inner-Laplace skewness diagnostic, the complementary
    # layer to diagnose_k above: whether the Gaussian inner Laplace on the
    # latent field is itself a good fit, scored at the fitted MAP grid cell
    # (one extra Newton solve via kernel_fn, not an importance batch). Default
    # probe scope is every arm's fixed-effects coefficients (see
    # .nlj_inner_skew_at_theta()); pass `control$skew_idx` (1-based, in the
    # joint [arm1_beta | arm1_re | arm2_beta | ... | blocks] latent layout,
    # see res$arm_layout) to probe additional indices. Declines to NA for a
    # genuinely coupled arm (a non-"separable" cell_coupling): the coupled
    # arm's per-obs likelihood is not what is being scored (see
    # build_joint_curvature3_fns), so this is honest, not a gap.
    diagnose_skew              <- isTRUE(control$diagnose_skew %||% TRUE)
    skew_idx                   <- control$skew_idx
    skew_correct               <- isTRUE(control$skew_correct %||%
                                          .nl_diag("skew_correct"))
    # Outer-grid node layout for the multi-block path. A CCD
    # places a central-composite design around the joint hyperparameter mode --
    # far fewer inner solves than the full tensor product (1 + 2d + 2^d vs k^d).
    # "auto" (default) uses the CCD only where the tensor blow-up bites hardest,
    # >= 4 transformable axes, and the tensor grid at <= 3 axes (the d = 3 grid
    # is cheap and a ridged 3-axis posterior integrates more reliably on it).
    # "ccd" lowers the threshold to >= 3 axes; "grid" always forces the tensor
    # product. The CCD auto-falls back to the tensor grid for an unguessable axis
    # (CAR_proper rho_car / non-BYM2 rho), a flat / ridged outer posterior, or a
    # degenerate mode-find; an active phi grid rides as a
    # tensor axis crossed on top of the CCD. Single-block joint
    # backends always use the tensor grid (CCD applies to the multi-block path).
    integration               <- match.arg(control$integration %||% "auto",
                                            c("auto", "ccd", "grid",
                                              "grid_adaptive"))
    # WITHIN-CELL construction for the reported per-axis hyperparameter
    # intervals. `integration` names how the grid was BUILT;
    # this names how each cell's mass is spread inside its own box when the
    # grid is read back. Both are numerical reconstructions of the same design,
    # so both are control knobs; the default reproduces the shipped read
    # byte-for-byte.
    within_cell               <- .nl_within_cell_mode(control$within_cell)
    # Inner-Newton curvature + PD enforcement for the (possibly indefinite)
    # joint mixture Hessian. "auto" (default) steps on the observed Hessian
    # wherever it factors as it stands and on the coupled cells' complete-data
    # expected information where it does not; "lm" escalates a diagonal ridge
    # until CHOLMOD factorizes the observed Hessian; "psd" eigen-clamps the
    # dense observed Hessian; "fisher" steps on the expected information
    # throughout. Every mode keeps the observed Hessian for the final
    # mode-pass log_det and SEs, and without a coupled cell "auto" and
    # "fisher" are "lm". Two C++ axes: hessian_pd_mode (0 = LM, 1 = PSD) and
    # step_curvature_mode (0 = observed, 1 = expected, 2 = auto).
    hessian                   <- match.arg(control$hessian %||% "auto",
                                           c("auto", "lm", "psd", "fisher"))
    hessian_pd_mode           <- if (identical(hessian, "psd")) 1L else 0L
    step_curvature_mode       <- switch(hessian, auto = 2L, fisher = 1L, 0L)
    # Inner-Newton Cholesky factor reuse (Shamanskii / chord method). For a
    # non-quadratic positive arm (e.g. beta cover) the latent Hessian changes
    # every inner iteration, so the plain Newton loop re-factorizes the sparse
    # Cholesky on each step -- the dominant per-grid-cell cost.
    # `inner_refresh = m` re-factorizes only every m-th inner step and reuses
    # the cached factor in between; the gradient stays exact and each step is
    # line-search safeguarded, so the converged mode is unchanged and the final
    # mode-pass Hessian (log_det / SEs) is always fresh. Applies to the sparse
    # joint path (large fields, LM curvature); the dense small-n_x path
    # re-factorizes cheaply and ignores it. Default 1 = re-factorize every step.
    inner_refresh             <- as.integer(control$inner_refresh %||% 1L)
    if (length(inner_refresh) != 1L || is.na(inner_refresh) ||
        inner_refresh < 1L) {
        stop("`control$inner_refresh` must be a single integer >= 1.",
             call. = FALSE)
    }

    # Outer-grid progress. The cpp entry is reached through
    # the polymorphic joint backends (`backend$call_kernel`) and the adaptive-
    # refinement closures, so the four progress knobs travel via a scoped
    # option rather than every backend signature; `.cpp_joint_multi` reads it at
    # the cpp boundary. Restored on exit so the option never leaks past the fit.
    .op_progress <- options(tulpa.nl_progress = .nl_progress_args(control))
    on.exit(options(.op_progress), add = TRUE)

    # CCD placement knobs (the ceiling on what placing a design may spend in
    # inner solves, and the curvature-reuse switch). `.joint_ccd_grid()` is
    # reached through the multi-block dispatch rather than from this frame, so
    # the resolved record travels on a scoped option like progress, the
    # checkpoint and the screening depth. Restored on exit so the option never
    # leaks past the fit.
    .op_ccd <- options(tulpa.ccd_placement = .ccd_placement_args(control))
    on.exit(options(.op_ccd), add = TRUE)

    # Grid-cell checkpoint/resume. Threaded to the cpp
    # boundary via a scoped option, like progress, so every backend / adaptive-
    # refinement kernel call within this fit shares one checkpoint file. A
    # `resume = FALSE` run starts over from an empty file; that removal belongs
    # to the fit, not to one solve, so it is taken by the caller
    # (`tulpa_nested_laplace_joint()`) before the first solve. This function
    # runs once per SOLVE -- a placement rescue refits -- and only publishes the
    # path. A resume (the default) keeps the file and the C++ layer loads its
    # completed cells.
    .ckpt <- .nl_checkpoint_args(control)
    .op_checkpoint <- options(tulpa.nl_checkpoint = .ckpt)
    on.exit(options(.op_checkpoint), add = TRUE)

    # The caller's outer-grid ceiling, on the same scoped-option transport: the
    # tensor grid is built inside .joint_dispatch_multi() and again on each
    # refinement pass, so the value is read where the grid is, not carried
    # through every backend signature.
    .op_grid_cap <- options(
        tulpa.nl_max_grid_cells = .nl_max_grid_cells(control))
    on.exit(options(.op_grid_cap), add = TRUE)

    if (!is.list(responses) || length(responses) < 1L) {
        stop("`responses` must be a non-empty list of arm specs.", call. = FALSE)
    }
    if (!is.list(prior)) {
        stop("`prior` must be a list with a `type` field, or a list of block specs (empty for none).",
             call. = FALSE)
    }
    # Parse user-specified regularizing hyperpriors on (sigma, alpha) once,
    # at the entry point. Validation errors raised here surface to the user
    # without going through the multi-block / single-block fork. `prior_sigma`
    # and `prior_alpha` also take the per-block shape on a list-of-blocks
    # `prior`; `prior_phi` names per-arm dispersion axes, which carry no block.
    multi_block <- .is_multi_block_prior(prior)
    fn_sigma <- .joint_parse_hyperprior(prior_sigma, "prior_sigma", multi_block)
    fn_alpha <- .joint_parse_hyperprior(prior_alpha, "prior_alpha", multi_block)
    fn_phi   <- .joint_parse_sigma_prior(prior_phi,   "prior_phi")

    # Detect multi-block prior (list-of-blocks). Routed to a separate
    # dispatch path that builds vector<LatentBlock> on the joint side.
    # Pruning toggle: the kernel takes prune_tol > 0 to mean "screen and
    # prune"; gate it here with the user-facing `prune` boolean so the toggle
    # composes cleanly (prune=FALSE always disables, prune=TRUE uses
    # prune_tol). prune_tol must be in [0, 1) -- clamp out-of-range to 0
    # rather than erroring, because a zero tolerance is the safe no-op.
    if (!isTRUE(prune)) {
        prune_tol_eff <- 0.0
    } else {
        prune_tol_eff <- as.numeric(prune_tol)
        if (!is.finite(prune_tol_eff) || prune_tol_eff < 0 ||
            prune_tol_eff >= 1) {
            stop("`prune_tol` must be a finite numeric in [0, 1).",
                 call. = FALSE)
        }
    }

    # Cheap-screen depth. Validated whatever `prune` says, so a misspecified
    # value is reported where the user set it rather than at the next fit that
    # happens to switch pruning on.
    if (length(screen_iters) != 1L || !is.numeric(screen_iters) ||
        !is.finite(screen_iters) || screen_iters < 1 ||
        screen_iters > .Machine$integer.max ||
        screen_iters != round(screen_iters)) {
        stop("`screen_iters` must be a single integer >= 1.", call. = FALSE)
    }
    screen_iters <- as.integer(screen_iters)

    # Cheap-screen depth reaches the cpp entry on the same scoped-option
    # transport as progress and the checkpoint: the entry is called through the
    # polymorphic backends (`backend$call_kernel`), the multi-block dispatch and
    # the refinement closures, so one option reaches every kernel call of this
    # fit instead of one scalar threaded through each of those signatures.
    # `.cpp_joint_multi` reads it at the cpp boundary. Restored on exit so the
    # option never leaks past the fit.
    .op_screen <- options(tulpa.nl_screen_iters = screen_iters)
    on.exit(options(.op_screen), add = TRUE)

    # A placement pilot's detecting grid (`.nl_pilot_control()`), read by the
    # multi-block main grid solve alone: that grid is solved only to find where
    # the posterior sits, so its screened surface is returned without the full
    # pass. The placement probes and every refinement solve of the same fit go
    # through other call sites and are full solves.
    .op_pilot_screen <- options(
        tulpa.nl_pilot_screen = isTRUE(control$pilot_screen))
    on.exit(options(.op_pilot_screen), add = TRUE)

    # The per-row predictive variance of the linear predictor rides the same
    # transport. It is what a grid-mixture replicate draws its WITHIN-cell
    # spread from (`posterior_predict()`), and a real per-cell solve sweep, so
    # a fit that will not be predicted from can decline it.
    .op_fitted_var <- options(
        tulpa.nl_fitted_var = isTRUE(control$fitted_var %||% TRUE))
    on.exit(options(.op_fitted_var), add = TRUE)

    # The override reaches the dense inner Newton through the multi-block
    # driver's call factory only. Refusing it elsewhere is what keeps it from
    # reading as a setting that was applied.
    if (inner_sparse_override != 0L && !.is_multi_block_prior(prior)) {
        stop("`control$inner_factorization` applies to the multi-block joint ",
             "path; this fit takes the single-block one.", call. = FALSE)
    }

    if (.is_multi_block_prior(prior)) {
        return(.joint_dispatch_multi(
            responses = responses, prior_list = prior, copy = copy,
            phi_grid = phi_grid,
            fn_sigma = fn_sigma, fn_alpha = fn_alpha, fn_phi = fn_phi,
            copy_atom_mass = copy_atom_mass, copy_slab = copy_slab,
            max_iter = max_iter, tol = tol, n_threads = n_threads,
            n_threads_outer = n_threads_outer,
            tile_warm = tile_warm,
            prune_tol = prune_tol_eff,
            x_init = x_init, verbose = verbose, store_Q = store_Q,
            keep_grid_hessians = keep_grid_hessians,
            force_sparse = force_sparse,
            inner_sparse_override = inner_sparse_override,
            cell_coupling = cell_coupling,
            hessian_pd_mode = hessian_pd_mode,
            step_curvature_mode = step_curvature_mode,
            diagnose_k = diagnose_k, diagnose_draws = diagnose_draws,
            pareto_k_by_arm = pareto_k_by_arm,
            pareto_k_threads = pareto_k_threads,
            k_bootstrap = k_bootstrap,
            k_tail_points = k_tail_points,
            k_conf_bands = k_conf_bands,
            diagnose_skew = diagnose_skew,
            skew_idx = skew_idx,
            skew_correct = skew_correct,
            subspace_debias = .subspace_debias_config(control$subspace_debias),
            cila = .cila_config(control$cila),
            inner_refresh = inner_refresh,
            integration = integration,
            within_cell = within_cell,
            local_ccd = local_ccd,
            adaptive_cutoff = adaptive_cutoff,
            adaptive_stride = adaptive_stride,
            adaptive_max_frac = adaptive_max_frac,
            adaptive_min_cells = adaptive_min_cells,
            placement_axes = placement_axes,
            hyperprior = hyperprior,
            adaptive_grid = adaptive_grid,
            adaptive_grid_edge_thresh = adaptive_grid_edge_thresh,
            adaptive_grid_max_passes = adaptive_grid_max_passes,
            var_of_means_consistency = var_of_means_consistency,
            axis_refine = axis_refine,
            grid_auto = grid_auto, phi_auto = phi_auto,
            timer = tm
        ))
    }
    if (!is.null(copy)) {
        stop("`copy` is not used on the single-block ",
             "tulpa_nested_laplace_joint() path. Declare the copy coefficient on ",
             "the arm itself: field_coef = list(name = \"alpha\", grid = <grid>). ",
             "(`copy` remains a multi-block argument -- a list of block specs, ",
             "each naming a distinct copy block.)", call. = FALSE)
    }
    if (is.null(prior$type)) {
        stop("`prior` must be a list with a `type` field, or a list of block specs (empty for none).",
             call. = FALSE)
    }
    .nl_refuse_projector(prior, "The single-block tulpa_nested_laplace_joint() path")
    type <- tolower(prior$type)
    backend <- .joint_backends[[type]]
    if (is.null(backend)) {
        stop("tulpa_nested_laplace_joint(): unsupported prior$type '", type,
             "'. Supported: ",
             paste(shQuote(names(.joint_backends)), collapse = ", "),
             ".", call. = FALSE)
    }

    arms <- lapply(seq_along(responses), function(k) {
        a <- responses[[k]]
        .normalise_joint_arm(a, k)
    })

    # `.resolve_copy` reads the per-arm `field_coef_axis` / `field_coef_const`
    # populated by `.normalise_joint_arm`; pass the normalised `arms` rather
    # than the raw `responses` list.
    cp <- .resolve_copy(arms, prior, type)
    arm_names <- names(responses) %||% paste0("arm", seq_along(responses))
    phi_axes <- .normalise_phi_grid(phi_grid, arm_names)
    grids <- backend$build_grids(prior, cp$has_copy, cp$alpha_grid, phi_axes)
    # How far each axis may be refined follows from where its nodes came from:
    # an axis the caller stated is integrated more finely over exactly that
    # range, one the engine placed may be followed out past its ends.
    # `control$axis_refine` overrides either, per axis. Resolved here, before the
    # first kernel call, so a mis-named axis is refused rather than reported
    # after the grid has been solved.
    #
    # A dispersion axis's provenance travels as the front door's `phi_auto`
    # record, NOT on the grid: the front door strips the markers before the
    # first fit (`.nl_phi_provenance()`, so nothing downstream sees an
    # attributed numeric), which left every `phi_<arm>` axis reading as stated
    # here whoever wrote it.
    axis_refine_modes <- .joint_axis_refine_modes(
        grids, cp, arms, phi_grid = phi_grid, arm_names = arm_names,
        phi_auto = phi_auto, prior = prior, grid_auto = grid_auto,
        user = .joint_check_axis_refine(axis_refine,
                                        .joint_spec_axis_names(grids, cp)))
    force_sparse <- .resolve_force_sparse(force_sparse, function() {
        lay <- if (is.function(backend$layout)) backend$layout(arms, prior)
        if (is.list(lay)) lay$n_x else NULL
    })
    # Per-cell fixed-effect covariance retention. The block
    # size and the field sum-to-zero groups follow from the latent layout, which
    # is fixed before the first solve, so the request goes in with the kernel
    # call and each cell's block is extracted inside that cell's own solve
    # rather than off a grid-wide store of every precision.
    fb_layout <- if (is.function(backend$layout)) backend$layout(arms, prior)
                 else NULL
    fixed_block_p <- if (isTRUE(keep_grid_hessians))
        as.integer(.joint_fixed_layout(responses)$n_fixed) else 0L
    fixed_block_constraints <- if (is.list(fb_layout))
        .joint_constraint_cols(fb_layout, fb_layout$n_x) else list()
    tm$mark("setup")

    # The hyperprior folded into log_marginal at the kernel-call boundary: the
    # caller's density on a role it names, the engine's default on every other
    # axis (`R/hyperprior_default.R`). Refinement decisions (edge scores, modal
    # cell selection, var-of-means thresholds) then all read the posterior, and
    # cells appended in refinement passes get the same prior through `hp_fn`.
    # Known from the grid alone, so a screened solve ranks its cells with it and
    # with their measure.
    hp_families <- stats::setNames(
        vapply(arms, function(a) as.character(a$family %||% ""), character(1)),
        arm_names)
    theta_grid_init <- backend$theta_grid(grids, cp$has_copy)
    hp_declared <- .hp_declare(theta_grid_init)
    hp_record <- function(cells) {
        .joint_hyperprior(cells, list(c(prior, list(type = type))), hp_families,
                          user = list(sigma = fn_sigma, alpha = fn_alpha,
                                      phi = fn_phi),
                          copy_atom_mass = copy_atom_mass,
                          hyperprior = hyperprior, declared = hp_declared)
    }
    hp_fn <- function(new_cells) hp_record(new_cells)$lp
    hp_init <- hp_record(theta_grid_init)
    specs <- .joint_axis_specs(grids, cp,
        user_priors = list(sigma = fn_sigma, alpha = fn_alpha, phi = fn_phi),
        copy_atom_mass = copy_atom_mass, copy_slab = copy_slab,
        axis_refine = axis_refine_modes,
        folded_axes = hp_init$axes)
    screen_offset <- if (prune_tol_eff > 0)
        .nl_screen_log_offset(theta_grid_init, list(hp_init), specs = specs)

    # A placement pilot's detecting grid is read off the screen alone
    # (`.nl_pilot_control()`); the full-grid fallback below never is.
    pilot_screen <- isTRUE(getOption("tulpa.nl_pilot_screen", FALSE))
    call_kernel_with_tol <- function(tol_prune) {
        backend$call_kernel(arms, prior, cp, grids, max_iter, tol,
                            n_threads, x_init, isTRUE(store_Q),
                            fixed_block_p = fixed_block_p,
                            fixed_block_constraints = fixed_block_constraints,
                            arm_names = arm_names,
                            n_threads_outer = n_threads_outer,
                            tile_warm = tile_warm,
                            prune_tol = tol_prune,
                            force_sparse = force_sparse,
                            cell_coupling = cell_coupling,
                            hessian_pd_mode = hessian_pd_mode,
                            step_curvature_mode = step_curvature_mode,
                            inner_refresh = inner_refresh,
                            screen_log_offset = if (tol_prune > 0) screen_offset,
                            screen_only = pilot_screen && tol_prune > 0)
    }
    res <- .joint_main_grid_solve(function() call_kernel_with_tol(prune_tol_eff))
    # Safety gate: if the cheap-pass ranking is unreliable (the screen's
    # argmax disagrees with the full-solve argmax, or the kept posterior
    # collapses onto a cell the screen badly mis-estimated), warn and fall
    # back to the full grid rather than silently returning a pruned answer
    #. A silently-wrong posterior must be impossible.
    if (prune_tol_eff > 0) {
        res <- .joint_prune_safety_gate(
            res, resolve_full = function() call_kernel_with_tol(0.0))
    }
    tm$mark("grid")

    res <- .nl_fold_hyperprior(res, list(hp_init))

    # --- generic-refinement glue (Step 3) -------------------
    # Adaptive grid + var-of-means consistency are now driven by the
    # axis-spec module in `R/hyper_grid_refine.R`. The joint driver builds:
    #   * `specs`     -- one hyper_axis_spec per outer-grid axis (log_scale /
    #                    bounds / refinable / refine_priority encoded once,
    #                    from `.joint_axis_specs`).
    #   * `kernel_fn` -- thin closure around `backend$call_kernel` that
    #                    converts matrix new_cells <-> paired-vector grids and
    #                    packs per-cell modes / n_iter / Q_csc_* into extras.
    #   * `extras_list` -- per-cell side data the kernel returned; refinement
    #                    carries it through, the warm-start chain reads
    #                    `extras[[anchor]]$mode`.
    # After refinement `.joint_glue_extras_to_res` puts the merged extras
    # back into `res$modes` / `res$n_iter` / `res$Q_csc_*_per_grid` so
    # downstream code (`.nl_posterior_moments`, `.nl_attach_axis_sd`,
    # the modes / Q consumers) reads the refined values unchanged.
    kernel_fn <- .joint_make_kernel_fn(arms, prior, cp, backend, max_iter,
                                        tol, n_threads, x_init, store_Q,
                                        arm_names,
                                        fixed_block_p = fixed_block_p,
                                        fixed_block_constraints =
                                            fixed_block_constraints,
                                        cell_coupling = cell_coupling,
                                        hessian_pd_mode = hessian_pd_mode,
                                        step_curvature_mode = step_curvature_mode,
                                        inner_refresh = inner_refresh,
                                        screen_tol = prune_tol_eff,
                                        screen_hp = hp_record)
    # `theta_grid_init` is the grid as it was declared, kept so the fit can
    # report the span it was given apart from the span refinement left it with.
    theta_grid_M  <- theta_grid_init
    log_marginal  <- res$log_marginal
    extras_list   <- .joint_init_extras_from_res(res)

    # Boundary / interior refinement (opt in), then the var-of-means consistency
    # pass. Sharply peaked axes (gaussian noise SD, beta phi at high n_pos, a
    # field SD narrower than its placed cell) collapse joint weight onto a single
    # grid cell, so `sum(w*x^2) - mean^2` on that axis is a floor at zero rather
    # than a spread. The consistency pass bisects the gaps the axis's mass sits
    # across, as levels laid in every row of the other axes that holds the
    # posterior, until the axis marginal's
    # ESS reaches the floor, so the merged grid carries the support the spread is
    # read off. The outer mode a placement refit was laid from is carried in by
    # the front door (`.joint_place_axes()`), which the pass lays a collapsed
    # axis's points at.
    # A round's cells are levels laid across the rows that hold the posterior,
    # independent of each other, so a round runs on the fit's outer width: the
    # kernel puts min(cells, width) of them in the team and gives each the
    # rest of the grant as inner threads, a three-cell round three cells of
    # ten threads on a 32-thread grant and a 78-cell round 32 cells of one.
    refine_kernel <- function(new_cells, warm_start = NULL, store_extras = FALSE,
                              ...)
        kernel_fn(new_cells, warm_start = warm_start, store_extras = store_extras,
                  n_threads_outer = n_threads_outer, ...)
    ref <- .joint_refine_outer_grid(
        theta_grid_M, log_marginal, extras_list, specs, refine_kernel, hp_fn,
        adaptive_grid = adaptive_grid,
        edge_thresh   = adaptive_grid_edge_thresh,
        max_passes    = adaptive_grid_max_passes,
        consistency   = var_of_means_consistency,
        axis_modes    = .nl_outer_mode_axes(getOption("tulpa.nl_outer_mode")))
    theta_grid_M  <- ref$theta_grid
    refining_axis <- ref$refining_axis
    refine_info   <- ref$adaptive_info
    # A cell whose inner Newton stopped short of a mode, on the base grid or on
    # a level a pass added, is solved again from its nearest converged
    # neighbour's mode before the grid is read.
    rs <- .joint_resolve_unconverged(theta_grid_M, ref$log_marginal, ref$extras,
                                     refine_kernel, hp_fn)
    log_marginal  <- rs$log_marginal
    extras_list   <- rs$extras
    res <- .joint_glue_extras_to_res(res, theta_grid_M, log_marginal,
                                      extras_list, refining_axis)
    res$nonconverged_resolve <- rs$info
    tm$mark("grid")                          # refinement inner solves

    res$theta_grid  <- theta_grid_M
    res$theta_names <- colnames(res$theta_grid)
    res$log_quad     <- .hyper_log_quad_weights(res$theta_grid, specs,
                                                refining = refining_axis)
    res$axis_support <- .hyper_grid_supports(res$theta_grid, specs,
                                             refining = refining_axis)
    res$axis_span    <- .joint_axis_span(theta_grid_init, res$theta_grid, specs,
                                         refining = refining_axis)
    res <- .nl_record_unconverged(res)
    res$weights     <- .nl_normalise_weights_safe(res$log_marginal, "outer grid",
                                                  log_quad = res$log_quad)
    res$log_hyperprior <- hp_fn(res$theta_grid)
    res             <- .nl_prune_merged_record(res)
    res             <- .nl_attach_evidence(res, res$theta_grid, specs)
    res             <- .nl_posterior_moments(res, paste0("joint_", type),
                                             within = within_cell)
    if (isTRUE(var_of_means_consistency)) {
        res$var_of_means_consistency_info <- ref$consistency_info
    }
    tm$mark("postproc")
    # The copy scale's point mass, reported beside the posterior mass left on
    # it. A point mass competes against a continuum, so how much of a "no
    # coupling" read is prior and how much is data is not visible from the
    # posterior alone; both numbers are needed to say.
    if ("alpha" %in% colnames(res$theta_grid)) {
        av <- as.numeric(res$theta_grid[, "alpha"])
        wv <- res$weights
        wv[is.na(wv)] <- 0
        res$copy_atom <- list(prior_mass = copy_atom_mass,
                              posterior_mass = sum(wv[av == 0]))
    }
    res$arm_layout  <- backend$layout(arms, prior)
    res$prior       <- prior
    res$responses   <- responses
    res             <- .joint_attach_observations(res, responses)
    res$copy        <- copy
    res$cell_coupling      <- cell_coupling
    res$adaptive_grid_info <- refine_info
    # Every axis's declared nodes, whatever subset of the tensor a pilot
    # thinning or the prune screen kept: the placement mode-find searches an
    # axis over these (`.nl_placement_mode()`).
    res$outer_axis_nodes <- lapply(
        stats::setNames(nm = colnames(theta_grid_init)),
        function(a) sort(unique(as.numeric(theta_grid_init[, a]))))
    tm$mark("postproc")
    # Outer Pareto-k-hat: re-evaluate the inner joint marginal at hyperparameters
    # drawn from the integrator's Gaussian proposal (reusing the generic
    # `kernel_fn` + `hp_fn` so no kernel-call machinery is duplicated) and
    # PSIS-smooth. Declines (NA -> quad-ESS) when an axis has unguessable
    # support (e.g. CAR_proper's rho_car).
    res <- .joint_attach_pareto_k_single(res, kernel_fn, hp_fn,
                                         max_iter   = max_iter,
                                         diagnose_k = diagnose_k,
                                         diagnose_draws = diagnose_draws,
                                         n_threads_outer = pareto_k_threads,
                                         pareto_k_by_arm = pareto_k_by_arm,
                                         k_bootstrap = k_bootstrap,
                                         k_tail_points = k_tail_points,
                                         k_conf_bands = k_conf_bands,
                                         placement_axes = placement_axes)
    res <- .nlj_inner_skew_at_theta(res, kernel_fn, skew_idx,
                                    compute = diagnose_skew)
    fixed <- .joint_fixed_layout(responses)
    res <- .nl_skew_correction_attach(res, fixed$n_fixed, skew_correct)
    tm$mark("diagnostics")
    # Per-cell fixed-effect mode + precision for the grid marginalization,
    # taken after every refinement pass has settled the grid so the cells it
    # reads are the cells the weights describe.
    res <- .joint_finalize_grid_fixed(res, fixed$n_fixed, keep_grid_hessians,
                                      keep_joint_prec)
    # `store_Q` was turned on for the closure, not by the caller: hand back the
    # one assembled precision, not the per-cell CSC scratch the caller never
    # asked for and would have to carry through every downstream object.
    if (keep_joint_prec && !isTRUE(control$store_Q)) {
      res$Q_csc_p_per_grid <- NULL
      res$Q_csc_i_per_grid <- NULL
      res$Q_csc_x_per_grid <- NULL
      res$Q_csc_n          <- NULL
    }
    # Subspace debias, after the per-cell fixed-effect pieces
    # are in place: the correction reports draws recombined from them, so it has
    # to see the settled grid the weights describe.
    res <- .nl_subspace_debias_attach(
        res, .subspace_debias_config(control$subspace_debias),
        redispatch = function(req) .joint_with_quiet_opts(
            kernel_fn(res$theta_grid, debias = req)),
        p_fixed = fixed$n_fixed, beta_names = fixed$names)
    # Corrected integrated Laplace, over the same settled
    # grid and through the same kernel closure.
    res <- .nl_cila_attach(
        res, .cila_config(control$cila),
        redispatch = function(req) .joint_with_quiet_opts(
            kernel_fn(res$theta_grid, cila = req)),
        p_fixed = fixed$n_fixed, beta_names = fixed$names,
        remoments = function(r) .nl_posterior_moments(
            r, paste0("joint_", type), within = within_cell))
    res$timing <- tm$timing()
    res <- .joint_attach_diagnose_cost(res, diagnose_k, diagnose_draws)
    .finalize_fit(res, backend = "nested_laplace_joint",
                  n_fixed = fixed$n_fixed, fixed_names = fixed$names,
                  extra_class = c("tulpa_nested_laplace_joint",
                                  "tulpa_nested_laplace", "list"))
}

# Inner-Laplace skewness diagnostic for the joint driver:
# the single-block counterpart of .nl_inner_skew_at_theta() (R/laplace_diagnostics.R),
# reusing `kernel_fn` (the SAME closure the outer Pareto-k diagnostic and the
# adaptive-grid refinement already call) so no separate re-dispatch machinery
# is duplicated -- one extra deterministic Newton solve at the fitted MAP grid
# cell, `compute_skew = TRUE`.
#
# Default probe scope is every arm's fixed-effects coefficients, read from
# `res$arm_layout$beta_start` / `$p` (.joint_layout(), already attached by the
# caller): the union of `(beta_start[k]+1):(beta_start[k]+p[k])` across arms
# k, generalizing the single-arm default (probe 1:p_fixed) to "every arm's
# beta", the direct multi-arm analogue.
#
# A coupled arm (cell_coupling != "separable" on that arm) has no per-obs
# oracle to score (build_joint_curvature3_fns excludes it, see
# laplace_newton_joint.h) -- its indices come back NaN, not a silently wrong
# 0. The joint MULTI-block path (nested_laplace_joint_multi.R, used when a
# fit carries a per-group RE / trend / arm-specific field block) has its own
# counterpart, `.nlj_multi_inner_skew_at_theta()`.
.nlj_inner_skew_at_theta <- function(res, kernel_fn, skew_idx = NULL,
                                     compute = TRUE) {
    res <- .inner_skew_decline(res, "not_requested")
    if (!compute) return(res)

    modal_theta <- .joint_modal_theta(res)
    if (is.null(modal_theta)) return(.inner_skew_decline(res, "backend_unsupported"))

    probe_idx <- skew_idx
    if (is.null(probe_idx)) {
        al <- res$arm_layout
        if (is.null(al) || is.null(al$beta_start) || is.null(al$p)) {
            return(.inner_skew_decline(res, "no_probe_indices"))
        }
        probe_idx <- unlist(lapply(seq_len(al$n_arms), function(k) {
            if (al$p[k] <= 0L) return(integer(0))
            (al$beta_start[k] + 1L):(al$beta_start[k] + al$p[k])
        }))
    }
    probe_idx <- as.integer(probe_idx)
    if (length(probe_idx) == 0L) return(.inner_skew_decline(res, "no_probe_indices"))
    res <- .inner_skew_decline(res, "solve_failed")

    warm     <- .joint_modal_mode(res)
    warm_arg <- if (is.null(warm)) NULL else list(mode = warm)
    # .joint_grids_from_cells() reads colnames(new_cells) to rebuild the named
    # grids list, so the 1-row matrix must carry them (.joint_modal_theta()
    # strips to a bare numeric vector).
    theta_mat <- matrix(modal_theta, nrow = 1L,
                        dimnames = list(NULL, colnames(res$theta_grid)))

    out <- tryCatch(
        .joint_with_quiet_opts(kernel_fn(theta_mat, warm_start = warm_arg,
                                         compute_skew = TRUE,
                                         skew_idx = probe_idx)),
        error = function(e) NULL)

    .inner_skew_attach_probe(res, out)
}

# Does a one-arm joint fit fill `fitted_eta_var`? The switch
# `tulpa_nested_laplace_joint(control$fitted_var)` sets, read in ONE place: the
# single-species kernel call below and the fused batch (`.cpp_joint_multi_batch()`)
# both reach it here, so a batched species cannot integrate under a different
# answer from the single-species fit it reproduces.
.nl_want_fitted_var <- function() isTRUE(getOption("tulpa.nl_fitted_var", TRUE))

# Thin wrapper over cpp_nested_laplace_joint_multi that injects the fit-scoped
# knobs set by tulpa_nested_laplace_joint: the outer-grid progress reporter
# (`tulpa.nl_progress`), the grid-cell checkpoint file (`tulpa.nl_checkpoint`),
# the cheap-screen depth (`tulpa.nl_screen_iters`) and the per-row predictive
# variance switch (`tulpa.nl_fitted_var`). Every backend / refinement call site
# routes through here, so these reach the cpp entry without threading scalars
# through the polymorphic backend interface. Options unset -> progress = FALSE,
# no checkpoint, the engine's own screening depth, and the per-row variance on.
# A fit running under `tulpa_joint_grid_batch()` has its main grid solve
# captured or served here (`.joint_grid_batch_intercept()`).
.cpp_joint_multi <- function(...) {
  served <- .joint_grid_batch_intercept(list(...))
  if (!is.null(served)) return(served)
  p <- getOption("tulpa.nl_progress", NULL)
  if (is.null(p)) p <- .nl_progress_args(list(progress = FALSE))
  cp <- getOption("tulpa.nl_checkpoint", NULL)
  checkpoint_path <- if (is.list(cp)) as.character(cp$path) else ""
  si <- getOption("tulpa.nl_screen_iters", NULL)
  if (is.null(si)) si <- .nl_screen("iters")
  do.call(cpp_nested_laplace_joint_multi,
          c(.joint_phi_args_to_kernel(list(...)),
            list(progress          = isTRUE(p$progress),
                 progress_every    = as.integer(p$progress_every),
                 progress_throttle = as.numeric(p$progress_throttle),
                 progress_file     = as.character(p$progress_file),
                 checkpoint_path   = checkpoint_path,
                 screen_iters      = as.integer(si),
                 compute_fitted_var = .nl_want_fitted_var())))
}