tulpa 0.0.150

0.0.150

• The cross-cell estimator asymmetry in local-CCD refinement is documented as a
  known and instrumented property rather than left reading as an expected
  consequence (#319). Weight conservation is the statement on a flat
  integrand; the integrand refinement is selected for is peaked, so a refined
  cell's mass is re-estimated by a 25-node rule while its unrefined siblings keep
  the base grid's single atom, and the refined cell's share rises. The header now
  says so, points at log_mass_ratio as the per-cell reading, and records what
  the asymmetry costs in coverage rather than in grid-internal metrics.

0.0.149

• A CI-coverage gate now reaches local-CCD refinement (#320). The
  joint recovery fixture builds one iid block, so it has one latent axis and
  .joint_local_ccd_engage() -- which requires four -- declined silently: every
  number measured on this path over four rounds came from grid dumps rather than
  from a coverage sweep. recov_fit_joint_local_ccd() is that sweep at four
  crossed groupings on a deliberately coarse four-level base grid, run twice on
  the same seeds with the refinement on and off, and a companion test asserts the
  refinement actually engaged rather than trusting the settings.
• Measured, 150 seeds x 2 coefficients per level: pooled coverage 0.8800 with
  refinement on against 0.8800 off at nominal 0.95, and 0.7300 against 0.7267 at
  0.80, against a standard error of 0.0126. What refinement moves is the width,
  in the direction #319 predicts -- the intercept's mean 95% interval
  is 2.9% narrower -- but that buys 0 of 300 trials at 0.95 and 1 at 0.80. Where
  it pays is the hyperparameter axis: the sigma_1 interval is more than
  fourfold sharper (0.2330 against 1.0590) with half the median bias (0.0245
  against 0.0592) while still covering 149 of 150 against a nominal 0.95.
• recov_sweep() takes the nominal level to judge at (z / level, defaulting
  to the 95% every existing gate runs at) and reports mean interval width, so a
  second level and the width effect are read through the shared harness rather
  than a parallel one.

0.0.148

• A locally refined cell now reports the ratio of the two estimates of its own
  mass it already carries (#323): the coarse atom
  Delta_c exp(ell_c) the base grid gave it against the refined cloud
  Delta_c sum_j delta_j exp(ell_j) its own nodes give it, as
  log_mass_ratio = logSumExp_j(log delta_j + ell_j - ell_c) over the full node
  set including the centre. That comparison is the embedded-rule local error
  indicator classical adaptive cubature uses to decide whether a subdivided
  region's estimate is comparable to its unrefined siblings', and this grid
  evaluates both rules already, so it costs no inner solve. Recorded per cell on
  local_ccd_info with the two masses it is formed from and max_node_weight,
  the share the single largest node takes of its own cell's refined mass, for
  refined and declined cells alike.
• The centring score gains its curvature-scaled form (#324). A
  central composite design identifies a full quadratic exactly, so the whitened
  Hessian sits in the same coefficient vector the gradient does, and
  mode_gain = 0.5 g' (-H)^-1 g is the nats the quadratic model predicts the
  log-density gains by moving the expansion centre to the cell's own fitted
  peak. The plain gradient norm offset says nothing about how sharply the
  log-marginal falls away in the direction it points, so two cells with the same
  norm and curvature an order of magnitude apart are displaced by very different
  amounts; the scaled form is comparable across cells and the unscaled one is
  kept for the cases the scaled one cannot be formed in. mode_gain declines to
  NA where -H is not positive definite: a cell whose fitted quadratic is not
  concave has no interior peak to be displaced from.
• The three scores are orthogonal and none of them gates: misfit is
  non-quadraticity, offset / mode_gain are off-centring, log_mass_ratio is
  mass correction, and skew_max reads misfit and nothing else.

0.0.147

• The outer-grid dump / rebuild harness is in the test suite
  (tests/testthat/helper-outer-grid-dump.R, #322). A candidate
  construction for the outer integration weights is pure post-processing of a
  fit that already ran, so outer_grid_dump() writes the grid state
  (joint_grid, log_marginal, dnode, weight_kind, the axis tags and
  domains, the support the read was taken off, and the summary the fit shipped)
  and outer_grid_rebuild() re-reads the per-axis summary under any weight
  vector. The read goes through the engine's own .nl_axis_quantiles() ->
  .nl_summary_quantile(), never a second copy of it, and the round-trip
  assertion in test-outer-grid-dump.R -- rebuild-with-own-weights equals the
  shipped read -- is what makes an offline difference attributable to the
  weights alone. It holds exactly (0.000e+00) on a tensor grid, a global CCD and
  a locally refined grid.
• outer_grid_noise_floor() estimates the scale below which a difference
  between two reads is not resolved by the grid, as the spread of the read under
  a weight-preserving coarsening of each axis's own atom set (consecutive atoms
  merged at their weighted mean carrying their summed weight). Total mass and
  each group's first moment are exactly preserved, so only resolution is
  removed. On a one-axis dump with a Gaussian outer log-marginal the floor
  bounds the read's true error against the closed-form quantiles at every
  resolution from 9 to 81 levels.
• The joint multi-block driver records dnode on the fit beside the integration
  weights it was folded into. Recovering it afterwards is a division by
  exp(log_marginal), which loses the scale and is undefined on a cell whose
  inner solve returned no finite marginal.

0.0.146

• The local-CCD cubic misfit score now reports the whitened gradient its own
  least-squares fit already estimated, and the refinement carries it per cell as
  offset / offset_declined on local_ccd_info (#321). The score
  puts the linear term in its own design columns, so a cell whose outer
  log-marginal is a perfectly good quadratic that simply is not centred on the
  cell fits exactly and scores near zero however steep the gradient across it:
  passing it certifies that the design can represent the cell, not that the
  cell's coordinate is a representative point of it. offset is the
  standardized displacement of the cell's own peak from the cell's coordinate,
  in units of the marginal spread the whitening used, and nothing gates on it --
  a gradient across the cell is a cross-cell estimator question, orthogonal to
  the local shape skew_max reads.

0.0.145

• Local CCD refinement of the joint outer grid now keeps a refined cell's node
  cloud only where the cell's own outer log-marginal is close to the quadratic
  the cloud was placed from, and puts the cell back as its own mass atom where it
  is not (#318). The refinement was a large win on an outer target
  that is quadratic in the transformed coordinate (summed absolute endpoint error
  against closed-form axis quantiles 7.3118 against 24.1142 for not refining, 48
  configurations of an equicorrelated Gaussian) and a net loss on a skewed one
  (26.2467 against 23.1874 over 27 configurations of a Gaussian copula with
  Gamma(2) marginals, 42.8578 against 38.1609 over 48).

  The mechanism is the cell's own non-quadraticity, and it is measurable from the
  design rather than inferred. A central composite design identifies a full
  quadratic exactly, so the least-squares residual of the nodes' measured
  log-marginals against intercept + gradient + Hessian in the whitened offset is
  the part of the cell the design cannot represent; the nodes are evaluated
  whatever the residual says, so the score costs no inner solve.
  .joint_local_ccd_misfit() reports it as a standardized cubic magnitude on the
  same convention the inner-Laplace gamma_3 uses, and on the Gaussian target it
  is identically zero in all 48 configurations while on every skewed family it
  exceeds 0.08.

  The threshold is .NL_DIAG$gamma3_ok (0.5), one number for the inner band and
  this gate because both are a standardized third-order departure from the
  Gaussian the approximation was placed from. Where it belongs was measured: on
  an eight-family ladder (the Gaussian target plus Gaussian copulas with
  Gamma(1), (2), (4), (8), (16), (32) and (64) marginals, 48 configurations each
  bar 32 for Gamma(1)), 0.5 is the only threshold on the ladder 0.01 to 2 that
  improves or ties every family. Per family, gated against refining
  unconditionally: 7.3118 / 7.3118 on the Gaussian, then 36.3139 / 37.2185,
  40.8376 / 42.8578, 31.9596 / 34.3650, 34.1830 / 36.8361, 44.4462 / 47.2740,
  61.8421 / 63.1718 and 83.2492 / 83.2492 down the ladder; pooled 340.1435
  against 352.2843, with 365.0831 for never refining. Lower thresholds score
  better pooled (0.175 gives 337.4142) by regressing on the two least skewed
  families. 399 of 1196 candidate cells are declined across the ladder, none of
  them on the Gaussian target.

  On the four-axis two-block fixture the gate is measurably neutral, which is
  what it has to be: over the 20 distinct refined configurations the summed
  absolute endpoint error against the converged m = 13 reference is 2.11091
  gated against 2.12847 unconditional and 3.12915 unrefined, and the largest
  single-configuration movement is 0.00696 against that reference's own 0.01716
  endpoint noise floor. Its per-cell scores there run 0.053 to 3.822, and the
  design-dominated m = 3 configuration reads 0.126, so it keeps its cloud.
  control$local_ccd$skew_max overrides the threshold; $local_ccd_info gains
  misfit, skew_max, cells_declined, misfit_declined and
  n_cells_declined, and a refinement whose every candidate declined leaves a
  plain tensor grid that reports theta_interval_read = "density".

• A locally CCD-refined joint outer grid now says what its per-axis
  hyperparameter intervals were read off, and how much of the support underneath
  them is a quadrature design rather than posterior mass (#317). It
  is the one node set carrying both kinds at once -- a carried-over base cell
  holds the mass of its own cell, a refined cell's replacement cloud holds a
  partition-of-unity share of its cell's mass placed at the design's radius --
  and it reports integration = "grid", so .nl_node_support() read it as a
  homogeneous density grid and nothing downstream could tell. The support is now
  named "mixed", the fit carries theta_interval_read and
  theta_interval_design_mass, and diagnostic_summary() surfaces the pair.

  The reported numbers are unchanged, on the measurement. Three replacement
  reads were scored against the converged m = 13 tensor reference (28561
  cells) of the four-axis two-block fixture, whose own noise floor is 0.01716 on
  the endpoints and 0.03853 on the widths: splitting the read into a mass CDF
  plus a per-cell moment-matched Gaussian, collapsing each refined cell's design
  block to one atom at its own weighted mean, and the #308 moment
  read. Summed absolute endpoint error over seven base grids (m = 3 ... 9,
  design_mass 0.930 down to 0.092) is 0.63446 for the shipped weighted
  quantile, 0.73159 for the collapse, 0.74464 for the split and 1.20402 for the
  moment read, against 1.00289 for not refining at all; over fourteen further
  configurations reached by varying local_ccd$max_cells on the same fixture,
  1.65810 / 1.86035 / 1.89038 / 2.82062. The collapse is ahead at the single
  design_mass = 0.930 grid (0.12022 against 0.16409) and behind at every lower
  one, which is the pattern a fix has to avoid. On analytic outer targets whose
  axis quantiles are closed form the same ordering holds: at design_mass >= 0.5
  on an equicorrelated Gaussian (19 configurations) 3.18024 for the quantile
  against 15.72792 for the collapse and 17.80607 for not refining. The moment
  read wins on that target because a Gaussian moment match is its exact family
  there; on a Gaussian copula with Gamma marginals, correlated the same way, it
  is the worst of the four (39.93958 against 26.24674 over 27 configurations,
  both measured with the refinement engaged unconditionally; under the
  #318 gate above the quantile scores 24.59641 there).

• tulpa_re_cov_nested() reports the median and 95% interval of every derived
  covariance quantity (sigma_i, rho_ij, Sigma_ij, in every block) from the
  moments its integration design reproduces, instead of from a discrete weighted
  quantile over the design's node positions (#308). A central
  composite design is a moment rule: its nodes sit where they reproduce the
  integrand's first two moments and carry no probability mass of their own, so
  the cumulative design weight across them is not a CDF. The discrete quantile
  clamped an out-of-support probability to the extreme node, which made every
  reported interval exactly the design's own extent -- at k covariance
  coordinates, theta_hat +/- 1.1 sqrt(k) posterior SDs, whose coverage is
  capped at 2 Phi(1.1 sqrt(k)) - 1 = 0.729 / 0.880 / 0.943 / 0.972 for
  k = 1 ... 4 no matter how much data the model is given. The interval is now
  moment-matched on each quantity's own coordinate (log for a scale or a
  variance, atanh for a correlation, the identity for a covariance) and mapped
  back, so scale intervals stay positive and asymmetric and correlation
  intervals stay inside (-1, 1). Measured sigma_1 coverage over 200 seeds
  per arm (poisson, 60 groups x 40 observations, nominal 0.95), varying only the
  random-effect block: k = 1 0.735 -> 0.950, k = 2 0.880 -> 0.945, k = 3
  0.950 -> 0.945 (binomial SE 0.015 - 0.031). The mean and sd columns are
  unchanged, tulpa_re_cov_gibbs() was never affected, and the tensor-grid
  layout (control$integration = "grid"), whose uniform cells do discretize the
  density, keeps the weighted quantile.

• .nl_wtd_quantile() takes an explicit outside policy for a probability
  beyond the support's own cumulative range. The default "clamp" is the
  cumulative-mass convention a sample uses and is unchanged; "na" withholds
  the number where the support is a quadrature design rather than posterior
  mass, so the clamp is a stated choice instead of a silent one.

• The nested random-effect-covariance recovery gate
  (test-re-cov-recovery.R) is raised from 75% to the 85% its Gibbs sibling is
  held to. The 75% was the defect above, measured and accepted rather than
  diagnosed.

• Every remaining consumer that reads a median and interval off a CCD-integrated
  outer grid now takes them from the moments the design reproduces, the same way
  tulpa_re_cov_nested() does since #308. The three are the joint
  multi-block per-axis hyperparameter summary (theta_median / theta_ci_lo /
  theta_ci_hi, #309), the inline spatial() bar field's per-block
  sigma / rho and its MCAR Sigma summary (#310), and
  spatial_range() / temporal_corr() (#312). Confirmed before
  changing anything: on a real three-block joint CCD fit the reported upper
  endpoint equalled the node maximum on every axis at 0.000e+00, and on a real
  MCAR fit four of the six derived quantities had BOTH endpoints exactly on the
  node extent.

  Measured coverage of the reported interval, gaussian response, 40 groups per
  factor and 1200 observations, varying ONLY the number of crossed random-effect
  blocks (which is the outer dimension d), 200 seeds per arm, nominal 0.95:
  d = 3 0.9350 -> 0.9450, d = 4 0.9900 -> 0.9700, d = 5 0.9650 -> 0.9300
  (binomial SE 0.007 - 0.018). Before, the mean interval width grew with the
  design -- 0.3131 / 0.3560 / 0.3899 across the three arms, in a data regime that
  did not change -- because the endpoints were the design's 1.1 sqrt(d)
  whitened-SD reach rather than a property of the posterior; after, it is
  0.3241 / 0.3251 / 0.3261, flat to 0.6%.

  On the MCAR bar field, whose simulated Sigma gives an explicit truth (8 x 8
  lattice, 20 observations per cell, 150 seeds, nominal 0.95), every derived
  quantity's coverage improves: sigma_1 0.9067 -> 0.9133, sigma_2
  0.7467 -> 0.8333, rho_12 0.7533 -> 0.8400, Sigma_12 0.8667 -> 0.9267,
  pooled 0.8211 -> 0.8767 over 900 trials. It does not reach nominal; the
  residual is the inner Laplace's attenuation of a spatial covariance at 64
  units plus the outer Gaussian grid's fit to a skewed log-Cholesky posterior,
  which this does not address.

• Each nested-Laplace hyperparameter axis names the DOMAIN its interval is
  formed on, read from the same per-axis registry the outer Pareto-k
  unconstrains with (.joint_axis_domains()): a positive scale on log, the
  BYM2 mixing weight on logit, an unconstrained coordinate (a copy alpha, an
  MCAR log-Cholesky entry) on the identity. A proper-CAR rho_car, whose support
  is the adjacency's eigenvalue interval, has no domain the engine will guess, so
  its interval is withheld rather than reported as the design's extent.

• ranef() on a tulpa_re_cov_nested() fit reports the SAMPLED values for a
  random-effect coordinate the subspace debias selected, instead of the Gaussian
  mixture the fit had stopped using (#314). The node mixture is
  reused rather than redrawn, so the group effects and the fixed effects
  marginalize one weighted node set, and a source column says per row which
  construction produced it ("sampled", "mixture" or "mode"). A fit whose
  S contains no random effect is untouched, consumes no random number, and
  reports what it reported before.

• The inner-k-hat identity test strips the $skew_correction record alongside
  $inner_* (#313). Every entry of that record is derived from
  gamma_3 or from the switch itself, so it necessarily differs with
  control$diagnose_skew on and off; the invariant the test asserts -- that the
  diagnostic consumes no randomness and changes no non-diagnostic field -- is
  unchanged.

• Every joint outer grid reports weight_kind, one entry per cell, saying
  whether that cell carries the mass of its own cell or an in-cell design weight
  (#311). A fit integrated by one rule reports one value throughout;
  a locally CCD-refined grid is the one support carrying both, and it now says so
  per cell instead of leaving a consumer to read one kind off integration.
  local_ccd_info gains n_design_nodes and design_mass, the share of the
  integration weight sitting on design-weighted nodes.

  The per-axis median and interval keep the weighted quantile on a refined grid,
  which is a measurement rather than an omission. Scored against outer targets
  whose exact axis quantiles are known in closed form, the refined grid's
  quantile beats the unrefined grid's own in 8 of 8 configurations on a
  diagonal-Gaussian target (d = 4/5, 5/7 levels per axis, two peak
  sharpnesses), so declining it would withhold a number strictly better than the
  one the same summary reports one refinement earlier. On a target the moment
  rule cannot fit by construction (each axis carrying a Gamma(k, 1) marginal on
  a log-tagged axis, 18 configurations over d = 4/5, 5/7/9 levels,
  k = 1/2/8), the moment read beats the quantile as often on an unrefined grid
  (12 of 18) as on a refined one (11 of 18), and refining improves the moment
  read in only 2 of 18 -- so that advantage belongs to coarse grids and
  near-lognormal targets, not to the mixed weights, and moving a refined fit onto
  it is a separate question about density supports. The mixing is bounded by the
  refinement itself: each node cloud is clamped to its cell's Voronoi half-box,
  so a refined cell's mass is redistributed only inside the cell the unrefined
  grid had collapsed onto one point.

• Local CCD refinement scales each node cloud by the refined cell's MARGINAL
  spread instead of its conditional one (#316). The cloud's scale
  came from a diagonal finite-difference stencil, and 1 / sqrt(-d2_j) is the
  spread along axis j with every other axis held at the cell; the per-axis
  summary reports marginal spreads, and on a correlated outer posterior -- a
  sigma-alpha copy ridge is one -- the two differ by sqrt(H_jj (H^-1)_jj). The
  stencil now also differences the cell's CORNER grid neighbours, which a cell
  interior on every axis always has and which the tensor base already evaluated,
  so sqrt(diag((-H)^-1)) costs no extra inner solve. A cell whose corners are
  missing (a grid a previous pass spliced nodes into) or whose local -H is not
  positive definite keeps the conditional scale.

  Only the scale changes. The design stays axis-aligned, the per-axis shrink to
  the Voronoi half-box and the node clamp are untouched, and so is the
  weight-conservation argument. That is deliberate: the summary reads a weighted
  quantile over the refined grid and on design weights a cumulative sum is not a
  CDF, so what it returns is close to the design's own per-axis EXTENT. An
  axis-aligned design puts an axial node at f_0 sd_j on coordinate j; a
  design rotated by the Cholesky factor of the same covariance puts it at
  f_0 L[j, k], and measured on an equicorrelated target with unit marginal SDs
  at rho = 0.8 the per-axis extents are 2.200 / 1.760 / 1.765 / 1.918 times the
  SD -- so a rotated design reports an interval that depends on the arbitrary
  order of the axes. Rotating was built and measured and is not what shipped.

  Scored against an equicorrelated Gaussian outer target on identity axes, whose
  axis marginals are standard normal whatever the correlation is, at four axes
  and three levels per axis (design_mass above 0.99), reported 95% width as a
  fraction of the exact 3.91993:

  rho     0.0     0.5     0.7     0.8     0.9
  before  1.0052  0.8587  0.6748  0.5543  0.3942
  after   1.0052  0.9740  0.9540  0.9527  0.9569
  

  Mean absolute endpoint error at rho = 0.8 goes 0.58240 -> 0.06182. The target
  is quadratic, so the stencil is exact there: the recovered marginal SDs match
  sqrt(diag(Sigma)) to 1.1e-15 and the conditional ones 1 / sqrt(diag(Q))
  exactly.

  On the package's own 4-axis multi-block fixture, refit on axis ranges that
  bracket its posterior and scored against a converged m = 13 tensor grid
  (28561 cells, 583 s, refinement off), mean absolute endpoint error per base
  grid, with the reference's own m = 11 against m = 13 movement as the noise
  floor (0.01716 on endpoints, 0.03853 on widths):

  levels m      3       4       5       6       7       8    total
  design_mass  0.931   0.289   0.222   0.254   0.169   0.121
  before       0.1797  0.1533  0.1219  0.0616  0.0537  0.0469  0.6170
  after        0.1641  0.1518  0.1219  0.0616  0.0537  0.0469  0.6000
  

  The improvement is concentrated where #316 is: at three levels per axis the
  reported mean width goes 0.67781 -> 0.71103 against a reference 1.09089, and
  since both sit below the reference that direction does not depend on the
  reference's exact value. It is modest there because that fixture's outer
  correlations are weak (the refined cell's local correlations are -0.460 to
  +0.158) and because at four of the six base grids the estimated -H is
  indefinite, so the conditional scale is kept. Four other candidates were built
  and measured against the same reference -- removing the shrink and clamp and
  deleting absorbed cells, calibrating the scale to the cell's own mass share,
  collapsing each cloud back to a mass atom, and rotating by the Cholesky factor
  -- and none of them beats this across the sweep.

  Also measured: the shrink to the Voronoi box, which #316 names as the
  mechanism, does not bind on that fixture at three levels per axis. The
  conditional scale there is (0.2120, 0.1034, 0.1954, 0.0388) against a
  half / node_reach of (0.3197, 0.1818, 0.3197, 0.1420), and a candidate that
  removes the shrink and the clamp entirely returns bit-identical numbers.

• local_ccd_info gains cell_share, the share of the base grid's integration
  weight each refined cell held before any node was placed. It is a different
  number from design_mass, which is the share the refined region holds after:
  the replacement nodes sit nearer the peak than the cell's own coordinate did,
  so refining raises it. Reading both separates how concentrated the base grid
  already was from how much the refinement concentrated it.

• A joint multi-block fit records which outer integrator the caller ASKED for
  and why the CCD did not run (#315). $integration names the
  integrator that ran, and .nl_node_support() keys the interval construction
  off it, so "grid" could not distinguish a tensor grid the caller chose from
  one a declined CCD fell back to; the reason existed and was thrown away
  outside a verbose message. $integration_requested carries the request and
  $integration_declined the reason -- NA_character_ when nothing was
  declined, otherwise "axis_count", "unguessable_axis", "degenerate_axis",
  "modefind_ridge", "modefind_boundary", "modefind_degenerate",
  "modefind_failed", "hessian_singular" or "hessian_not_pd". The
  cell-count warning drops its "set control$integration = \"ccd\"" advice on a
  fit that already asked for one and was turned down, and names the decline
  instead.

0.0.144

• The subspace debias reaches the grid and joint nested backends
  (#306, the follow-up to #304). control$subspace_debias is now
  accepted by tulpa_nested_laplace() -- both the single-block kernels
  (icar / bym2 / car_proper / rw1 / rw2 / ar1) and the multi-block driver -- and
  by tulpa_nested_laplace_joint() on both its single-block and multi-block
  paths, with the same settings and the same meaning as on
  tulpa_re_cov_nested().

  The selector costs nothing new here. control$diagnose_skew (on by default)
  already re-dispatches the kernel at the fitted MAP cell and attaches both
  inner scores, so .subspace_bands() reads the per-index gamma_3 and inner
  Pareto-k-hat off the fit instead of paying for a probe of its own -- the one
  thing it had to learn is to prefer the stored per-index k-hat over re-fitting
  the raw importance curve, which a nested fit does not retain. The correction
  itself re-runs the settled grid once with the sampler on, because the
  corrected shape is a property of each cell and which cell is the MAP is only
  known after the first pass. A corrected fit then reports $draws -- each
  cell's Metropolis sample for the selected coordinates, the rest from the
  Gaussian conditional given them, mixed by the grid weights through the same
  .re_cov_nested_beta_draws() the RE-covariance backend uses -- instead of the
  Gaussian-mixture moments, and every coefficient-facing method reads them
  through the accessor it already used.

  Measured against the exact conditional posterior, computed by two-dimensional
  quadrature outside the engine. Both fixtures put the latent block where no
  observation reaches it, so the conditional posterior factorises into the
  fixed-effect target and an independent Gaussian and the quadrature is exact
  for the reported coefficients.

  Grid backend, rare-event binomial logit (n = 120): total absolute
  interval-endpoint error 0.5229 -> 0.1883 (mean of 5 seeds, sd 0.0672), a 64.0%
  cut; the reported centre moves -2.526 -> -2.63 against an exact posterior mean
  of -2.635. S is the intercept, selected on an inner k-hat of 0.705 with
  gamma_3 = -0.375 against an exact skewness of -0.388.

  Joint multi-block, the #300 coupled occupancy fixture: 0.6189 -> 0.163 / 0.341
  / 0.313 over three seeds, the centre -0.240 -> about -0.15 against an exact
  -0.1437 and the scale 0.389 -> about 0.447 against an exact 0.479.

  Coverage, with the exact posterior's OWN coverage as the reference rather than
  the nominal level -- the correction targets that posterior, so reproducing its
  coverage is the success condition and matching nominal is not. Grid backend,
  400 seeds, per coefficient:

  level	coef	exact	plain	corrected
  0.95	beta0	0.9050 (.0147)	0.9750 (.0078)	0.9150 (.0139)
  0.95	beta1	0.9300 (.0128)	0.9550 (.0104)	0.9225 (.0134)
  0.80	beta0	0.7600 (.0214)	0.8225 (.0191)	0.7575 (.0214)
  0.80	beta1	0.7950 (.0202)	0.8375 (.0184)	0.7800 (.0207)

  The corrected rate is within 0.52 SE of the exact posterior's on all four; the
  plain Laplace sits 2.4 to 4.8 SE above it. Its apparently better agreement
  with the nominal level is an over-wide, mis-shaped Gaussian, not accuracy.
  Joint multi-block, 200 seeds with S pinned to both coefficients, same
  reference: 0.9500 / 0.9600 / 0.7900 / 0.7400 exact against 0.9598 / 0.9598 /
  0.7940 / 0.7387 corrected and 0.9749 / 0.9548 / 0.8191 / 0.7739 plain (one
  seed's fit reported an NA bound, so n = 199 for the fits and 200 for the
  reference).

  Cost, over those sweeps: 0.0324 s (SE 0.0005) -> 0.1765 s (SE 0.0011) on the
  grid backend at mean |S| = 1.060, and 0.0662 s (SE 0.0012) -> 0.4019 s
  (SE 0.0045) on the joint one at |S| = 2, i.e. 5.45x and 6.07x. Roughly half of
  that is the second grid pass and the rest the sweeps themselves.

  An EMPTY selection leaves every backend bit-for-bit identical to the plain
  path, asserted per backend (log-marginal, weights, modes, the per-cell
  fixed-effect pieces, summary() and vcov()) in the new
  tests/testthat/test-subspace-debias-backends.R.

• The band selector under-flags on the coupled fixture, measured. On the
  #300 coupled occupancy fit both coefficients band good on both inner scores
  -- gamma_3 reads 0.256 and -0.126, the inner k-hat 0.378 and 0.329 -- while
  the exact skewness of the occupancy coefficient is 1.198, which is
  unreliable. gamma_3 recovers 0.21 of it there, below the 0.564 to 0.943
  range #304 measured on the separable fixtures, and the inner k-hat over 256
  one-dimensional draws does not separate that target from a Gaussian either. So
  on that fit the selector takes nothing and control$subspace_debias$idx is
  what pins the set; the joint numbers above are the pinned run. The ok band
  floor is unchanged -- #304 measured it, and one fixture where a lower-bound
  estimator undershoots harder than usual is a fact about the estimator, not a
  reason to move a threshold that was itself set on measurement.

• Threading. The sampler draws from R's RNG, so a debiased outer grid is
  integrated serially whatever n_threads_outer asked for, on both the
  single-arm driver and the sparse joint one. The cheap warm-start screen never
  runs the correction, and a joint cell 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, for
  the same reason.

• One request, one unwrap, one assignment point. The debias request travels
  through every nested kernel entry as ONE R list (idx plus the sweep budget)
  rather than four parallel arguments; unwrap_debias() / DebiasRequest
  (src/laplace_spec_fit.h) turn it into the solver's options at each entry,
  and run_subspace_debias() (src/subspace_debias.h) guards, runs and records
  the outcome for all three Newton loops -- the single-arm spec loop and both
  joint loops -- so the "empty index set is a no-op" contract and the result
  mapping are written once. cpp_laplace_fit_multi_re()'s four #304 arguments
  are folded into the same list; pre-release, no shim.

• control$subspace_debias is left refused at n_quad > 1 on
  tulpa_re_cov_nested(), which is the honest answer rather than a gap. The
  adaptive Gauss-Hermite inner marginal integrates each group's random effects
  out, so at the fitted point there is no conditional latent field: no joint
  precision to take Sigma e_i from, hence no Gaussian-conditional-mean surface
  and nothing for the sampler to move along. What could be corrected there is
  the outer optimum's own Laplace approximation, which is a different
  construction on a different density and would not share this machinery.

• Reporting the SAMPLED values per group for a random effect the closure pulled
  into S is split out as #314: the draws exist per node but are not
  recorded on the fit, so ranef() still reports the Gaussian mixture there.

• Fixed: tests/testthat/test-inner-pareto-k.R asserted that a fit is
  bit-for-bit identical with control$diagnose_skew on and off after stripping
  inner_* and timing, which stopped holding when #302 added the
  gamma_3-derived $skew_correction record two commits later. The record is
  diagnostic-derived, so it is stripped alongside them (#313).

0.0.143

• The joint tier's fixed-effect block is extracted inside each cell's own
  solve (#307). Filling $grid_modes / $grid_hessians on a joint
  fit (#305) read the block off the cell precision, so the joint kernels ran with
  store_Q internally and the whole outer grid's precision was alive at once
  between the kernel call and the extraction. Both joint Newton loops now take
  the request -- the leading block size plus the field sum-to-zero groups, both
  fixed by the latent layout before the first solve -- and return the block on the
  LaplaceResult re_cov contract, which the grid driver emits per cell. The
  dense loop builds one cell's CSC, extracts, and releases it; the sparse loop
  reads the builder's own CSC, so there the block costs no precision copy at all.
  store_Q is once again the caller's own knob, passed straight through.

  Measured on an ICAR chain fixture (n_fixed = 8, 40-cell grid, R heap
  high-water over three paired runs): peak 49.9 vs 51.5 MB at n_x = 408,
  53.6 vs 58.2 MB at 2008, and 62.1 vs 71.9 MB at 6008 -- the saving tracking the
  grid's precision (1.63 / 4.07 / 8.54 MB, i.e. 41.7 / 104.2 / 218.7 KB per cell)
  and growing with the field where the retained block does not (868 bytes per
  cell at every size). Fit time is unchanged: 0.27 / 0.73 / 1.73 s against
  0.31 / 0.72 / 1.75 s.

  The blocks are byte-identical to what cpp_joint_inner_vcov_blocks() returns
  for the same cells -- same bytes in, same routine -- so summary(),
  confint() and vcov() report exactly what they did in 0.0.142, and draws,
  modes, weights and log_marginal are untouched.

• Local-CCD refinement no longer costs a joint fit its intervals
  (#307). The node solves carry their own fixed-effect block through
  the splice alongside the inner modes, so a refined grid reports instead of
  recording grid_fixed_declined = "local_ccd_refined". Refinement itself is
  unchanged.

• One extraction algebra behind all of it: src/inv_block_extract.h holds the
  conditioning-by-kriging constraint correction and the diagonal-block
  extraction, both templated on a solve oracle. laplace_newton.h's
  inv_block_layout path drives it against the live Newton factor;
  extract_inner_vcov_block_cell() drives it against a factorized cell with the
  constraint; the joint loops go through the latter.

0.0.142

• A joint fit reports uncertainty on its fixed effects (#305).
  summary(), confint() and vcov() on a tulpa_nested_laplace_joint() fit
  reported the point estimates and NA for every standard error and both
  bounds, on both the single-block and the multi-block path. The grid
  marginalizer .nested_fixed_moments() reads one representation --
  $grid_modes and $grid_hessians, the per-cell fixed-effect mode and
  marginal precision -- and the joint driver stored neither, so a joint fit
  could not put an interval on any coefficient.

  Both joint drivers now fill that same pair through one shared helper,
  .joint_attach_grid_fixed(), so the two tiers reach the one marginalizer
  rather than growing a second one. Both joint layouts stack every arm's
  coefficients as a contiguous prefix of the latent vector, so the extraction
  is arm-aware by construction: it takes the whole 1:n_fixed block in one
  pass and reports each arm's coefficients under its own name. The per-cell
  block comes from cpp_joint_inner_vcov_blocks(), the joint tier's existing
  per-cell inner-covariance extraction, so the reported covariance is the
  field-constrained one the fit's own tulpa_posterior_draws() mixture is
  generated from. The reported covariance is the law of total variance over the
  outer grid and so carries both the within-cell curvature and the between-cell
  hyperparameter spread.

  Measured against references outside the engine. On the #300 coupled fixture,
  whose log posterior is written independently in R, the reported covariance
  matches the inverse numerical Hessian of that density at the mode to 1.4e-09
  relative. A one-arm joint fit and the single-block fit of the same model --
  the same data, block and grid -- now agree to 5e-08 on every coefficient and
  standard error, where the joint side previously produced NA; that
  equivalence is asserted for poisson, binomial and gaussian. Over a
  multi-cell ICAR grid the mixture matches the independent R implementation of
  the same law-of-total-covariance (.joint_mixture_moments()) to 8e-17, and
  200000 draws of the fit's own posterior mixture reproduce its standard errors
  to 0.1%. CI coverage is judged by the existing recovery sweep with the joint
  fitter substituted rather than by a second harness.

  Retention is control$keep_grid_hessians (default TRUE), and costs
  O(n_fixed^2) per cell: 860 bytes per cell at n_fixed = 8, unchanged as
  the field grows from n_x = 408 to n_x = 6008, with no measurable fit-time
  overhead (-0.1% over three paired runs). Reading the block needs the cell
  precision, which the kernels now keep during the fit and drop again unless
  control$store_Q asked for it; that transient is the existing store_Q
  peak, ~64 bytes per latent per cell, and #307 tracks removing it
  by extracting the block inside the joint Newton loop.

  Draws, modes, weights, log_marginal and the hyperparameter moments are
  bit-for-bit identical with the retention on and off -- this adds reporting,
  not inference. Where the retention cannot be trusted it declines with a
  reason on $grid_fixed_declined instead of going quiet: "not_requested"
  when switched off, and "local_ccd_refined" when local-CCD refinement
  rewrote the outer grid after the cells were stored (local-CCD keeps
  precedence, so no existing fit changes).

• .nested_fixed_moments() skips a grid cell that carries no integration
  weight. A pruned cell with no retained block previously turned the whole
  marginalized covariance into NA.

0.0.141

• The reliability band is now the debias SELECTOR: exact MCMC runs on only the
  misfit directions (#304). Escalation used to be whole-fit and
  all-or-nothing (tulpa_re_cov_nested -> tulpa_re_cov_gibbs, or a grid
  refinement), with nothing saying which directions needed exact treatment, so
  every coordinate paid the sampler's price including the ones the Gaussian
  already fits. The per-index inner diagnostics are already a map of exactly
  that. control$subspace_debias (default FALSE) bands every probed index,
  takes S = the misfit set, and corrects only x_S by Metropolis with
  x_{-S} carried at its Gaussian conditional.

  The sampled surface is x(u) = mode + V Sigma_SS^{-1} L u with
  V = Sigma E_S -- the q-dimensional generalization of the one-dimensional
  conditional-mean curve both inner diagnostics already walk, reusing
  inner_probe_column() rather than a second solve. The Gaussian restricted to
  it is exactly N(0, I) in u, so the walk is spherical and the Laplace
  shaping lives in the coordinates. Random-walk Metropolis, not NUTS: each
  evaluation is one call of the Newton loop's own penalized objective (O(N), no
  factorization, no derivative), so a gradient sampler would buy nothing and
  would need a derivative the loop does not expose along the surface.

  Selection is at the ok band, one step below the unreliable band the
  reporting layer flags on, because gamma_3 is a LOWER bound on the true
  skewness (0.564-0.943 of the exact value across the engine's own fixtures), so
  selecting at the reported boundary would leave genuinely misfit coordinates
  uncorrected. The inner importance k-hat (#303), which needs no derivative and
  does not undershoot the same way, is folded in as the worse of the two -- on
  the rare-event sweep below it is what bands 225 of 400 intercepts
  unreliable where the mean gamma_3 is only -0.597.

  MEASURED against an exact reference (Bernoulli random intercept at fixed RE
  SD, the group intercepts integrated out by Gauss-Hermite and p(beta | y)
  marginalized on a grid): the exact intercept marginal is
  mean -4.0050, 95% (-6.5082, -2.2898); the Laplace Gaussian gives
  mean -3.4025, (-5.3153, -1.4897), total endpoint error 1.9930; the
  correction on S = {intercept} gives mean -3.7675, (-6.1446, -2.0611),
  endpoint error 0.5923 -- a 70.3% reduction. Residual bias remains, as expected
  from a lower-bound skewness estimate and a Gaussian conditional.

  MEASURED against correcting EVERY coordinate, which is the question of whether
  a subspace is enough (200 seeds x 2 coefficients, Bernoulli random intercept,
  60 groups of 3, at the true RE SD, against the exact quadrature marginal,
  nominal 0.95): plain Laplace 0.9050 (se 0.0147), subspace debias 0.9275
  (0.0130) at a mean |S| of 0.945 coordinates, full-S debias over all 62
  latent coordinates 0.9225 (0.0134). Correcting about one coordinate recovers
  what correcting all 62 recovers, and costs 0.313 s against 0.461 s.

  MEASURED against the full Gibbs debias (400 seeds, rare-event binomial-logit
  with a random intercept, pooled over both coefficients): at nominal 0.95,
  plain Laplace 0.9738 (se 0.0057), subspace 0.8662 (0.0120), full Gibbs 0.8888
  (0.0111) -- subspace within 1.4 standard errors of the full debias; at nominal
  0.80, plain 0.8738, subspace 0.7175 (0.0159), full Gibbs 0.7037 (0.0161),
  within 0.6 standard errors. Cost 0.468 s against the full debias's 1.287 s,
  2.75x cheaper. On the small-group binary RE fixture every probed
  coordinate bands good (max |gamma_3| 0.236 over all 122 latent
  coordinates), S is empty, and the fit is the plain one.

  A second whole-fit sweep on a denser fixture -- Bernoulli random intercept, 60
  groups of 3 (n = 180), beta = (-2.5, 1), sigma_u = 1, 200 seeds, all three
  backends at their defaults on the same data -- does NOT reproduce that match,
  and the reason is worth stating rather than averaging away. Beta coverage
  pooled over both coefficients at nominal 0.95: plain nested 356/400 = 0.8900
  (se 0.0156), subspace 359/400 = 0.8975 (0.0152), full Gibbs 376/400 = 0.9400
  (0.0119) -- a 2.2 standard-error gap, at 1.074 s against 8.123 s. The
  correction is not what falls short there. On the SAME 200 seeds its own layer,
  conditional coverage at the true sigma, goes 0.9050 plain -> 0.9275 subspace
  against 0.9225 for correcting every one of the 62 latent coordinates, so one
  coordinate recovers what all 62 do, at 0.313 s against 0.461 s. What is left
  is outer: the nested path's sigma_1 interval covers 150/200 against the
  Gibbs sampler's 199/200 and its intercept interval is 26% narrower, and
  turning the correction on moves neither number. That is a different layer,
  and the two backends are not even integrating the same hyperprior (the
  conjugate Sigma | b draw cannot take the flat default); #308
  separates it.

  Against exact Gauss-Hermite quadrature on that fixture (24 seeds x 2
  coefficients, max grid tail mass 8.4e-14) total absolute endpoint error is
  plain 24.2040, subspace 21.2363 (-12.3%), every-coordinate 7.7709 (-67.9%).
  So the band-selected subspace recovers the full correction's COVERAGE while
  recovering about a fifth of its endpoint accuracy: it puts the interval in the
  right place without fully fixing its shape. Closing that remainder is what the
  coupling closure below would do, and only by growing S to nearly the whole
  coupled block.

  The COUPLING CLOSURE (grow S by the precision-graph neighbours whose partial
  correlation with a member exceeds a threshold) is implemented and was measured
  both ways rather than assumed, which is what the issue asked for. At the
  default threshold it changes nothing: against the exact marginal it moves the
  endpoint error 0.5923 -> 0.5651, a difference of 0.027 against a combined seed
  standard error of 0.038, and across the 400-seed sweep it fires on 163 seeds
  yet leaves coverage identical on 1572 of 1600 seed-coefficient-levels. The
  reason is that the partial correlations between a fixed effect and the random
  effects only run about 0.09 to 0.25 on these models, so a threshold in the
  usual "strong coupling" range never bites.

  Lowering it far enough to bite does move the finer metric, and that is worth
  stating precisely rather than glossing: on a 14-coordinate fixture the total
  endpoint error against the exact marginal falls 4.43 -> 1.21 only once the
  threshold reaches 0.05, at which point |S| has grown to 13.4 of 14 -- the
  full debias wearing a different name rather than a subspace one. So
  conditioning x_{-S} on the Gaussian does NOT reproduce the exact marginal
  endpoint for endpoint; it removes about 70% of the Gaussian's endpoint error
  at |S| = 1 and the rest is not cheaply recoverable by growing S.

  On the arbiter the issue actually names -- interval coverage -- that residual
  does not show: at 200 seeds the |S| = 1 correction and the all-62-coordinate
  correction cover 0.9275 and 0.9225, indistinguishable. Coverage is the coarser
  of the two metrics, and the closure is off by default because nothing measured
  here asks for it. It stays available as closure = TRUE or an explicit
  threshold.

  S is recorded on the fit as subspace_debias (selected indices, the
  per-index band table they were read from, what the closure added, and the
  per-node acceptance), so the escalation is auditable rather than implicit. An
  empty S is not a special case of anything: the sampler is never entered, no
  random number is consumed, and the fit is bit-for-bit the plain Laplace fit --
  asserted on both the solver and the front door.

• One random-walk Metropolis definition, not two (src/rwmh.h). The
  starting scale 2.4 / sqrt(d), the Roberts-Gelman-Gilks target acceptance,
  the Robbins-Monro burn-in adaptation and the accept test were written out
  inline in the covariance Gibbs sweep and would have been written out again for
  the subspace debias. They are now one set of primitives both consume.
  rw_accept() draws its uniform unconditionally so a sweep consumes exactly
  one uniform per test whatever the ratio is, which is what keeps the migrated
  Gibbs sweep's RNG stream unchanged.

0.0.140

• gamma_3 is now consumed, not only graded: the inner-Laplace marginals can
  be skew-corrected (#302). The cubic term was computed, banded
  and printed, and nothing read it -- so the inner layer was nested
  approximation with no debias, the position this engine is designed against,
  one layer in from where that argument is usually made. summary() and
  confint() on a nested-Laplace fit run with control$skew_correct = TRUE now
  report Cornish-Fisher marginal quantiles at each coefficient's own gamma_3,
  gated to the good / ok bands, and the Gaussian quantiles everywhere else.
  $skew_correction records the per-coefficient gamma_3, band and
  eligibility; a skew_applied attribute on summary() / confint() records
  what was used at the requested level. Wired through tulpa_nested_laplace()
  and both tulpa_nested_laplace_joint() paths. (At this release a joint fit
  recorded the correction without showing it, because the joint driver retained
  no per-cell fixed-effect Hessians for the grid-marginalized covariance;
  #305 supplies them in 0.0.142 and the correction applies there.)

  Rue, Martino & Chopin (2009) Sec 3.2.3 fit a skew normal here, under three
  constraints -- mean gamma^(1), variance 1, third log-density derivative at
  the mode gamma^(3). Two of those inputs exist in this engine and one does
  not: gamma^(1) comes from their denominator expansion, which is diagonal
  only in their augmented x_j == eta_j representation
  (src/inner_laplace_skew.h carries the reason). A skew normal fitted on the
  cubic term alone is therefore a different construction from theirs, and its
  attainable skewness saturates at |skewness| ~ 0.995 with the shape parameter
  diverging as that bound is approached -- inside the very band the correction
  is gated to. The Cornish-Fisher expansion is the quantile-side inverse of the
  same Edgeworth series gamma_3 is the leading term of, is linear in
  gamma_3 so it does not saturate, and returns quantiles directly.

  The correction is skewness-only and therefore partial, which is measured
  rather than asserted. Against exact quadrature quantiles of rare-event
  binomial-logit posteriors it cuts total absolute endpoint error from 2.4931 to
  1.3837 (44.5%), improving both endpoints in every case. On CI coverage over a
  small-group Bernoulli random-effect fixture (N = 48, 200 seeds x 2
  coefficients) it is directionally right and immaterial: nominal 0.95, Gaussian
  0.9650, corrected 0.9600; nominal 0.80, 0.8050 vs 0.8075; nominal 0.50, 0.4950
  vs 0.5000 -- every difference inside one standard error. Two reasons the
  coverage gain is smaller than the marginal gain: gamma_3 is a lower bound on
  the true skewness (0.875-0.943 of it on the cases above), and a biased Laplace
  mode stays biased because the location term is not computed. The correction
  is therefore OFF by default (.NL_DIAG$skew_correct); the coverage
  measurement does not justify defaulting it on. Draws, modes, weights and every
  other field the solve produced are bit-for-bit unchanged either way -- this is
  post-processing on the reported quantiles.

  New: .nl_skew_marginal(), .nl_skew_by_fixed(), .nl_skew_correction_attach()
  (R/laplace_diagnostics.R), src/cornish_fisher.h / .cpp,
  tests/testthat/test-inner-skew-correction.R, and a paired
  corrected-vs-Gaussian coverage gate in test-nested-laplace-recovery.R.

0.0.139

• gamma_3 now scores coupled multi-predictor likelihoods instead of
  declining on them (#301). The cubic Edgeworth term assumed a
  log-likelihood that is a separable sum of one-eta terms, so every unit reading
  several linear predictors at once -- a zero-inflation mixture's (count, zi)
  pair, a CellCouplingSpec cell's arms (tulpaObs's occu_cover) -- had no
  per-eta third derivative and came back NaN for good. The expansion is
  unchanged; only the contraction widens, to
  sum_units sum_{a,b,c} T^{abc} u^a u^b u^c with T the unit's third
  derivative in its linear predictors and u the eta response to Sigma e_i.
  The separable case is the one-coordinate special case of it.

  T is never materialised (src/curvature3_contract.h). Partition the unit's
  coordinates into K blocks and the contraction equals
  sum_a d/ds [u' L''(e + s u^(a)) u] at s = 0, because moving along block
  a's slice of the direction differentiates exactly that block's coordinates.
  Each term is one central difference of the Hessian the likelihood already
  returns for the Newton solve, so the whole tensor costs 2K extra evaluations
  per unit and no storage, at any block sizes. For a CellCouplingSpec that
  Hessian is the analytic CellDerivs block, so this is one finite-difference
  layer on an exact quantity, not a difference of a difference.

  The step is scaled PER BLOCK off that block's own eta magnitude, matching the
  eta-space step the scalar working-weight fallback takes. Measured against a
  five-point third derivative of the cell log-density: identical to a single
  global step while the arms share an eta scale, and 1.8x more accurate once one
  arm's |eta| is 67x the other's. The contraction is symmetrised over index
  permutations; for this block decomposition that is algebraically the plain sum
  (the three relabelings coincide), so it buys robustness at a block whose own
  quotient could not be formed rather than variance reduction.

  Verified against the exact posterior, not asserted: on the coupled two-arm
  occupancy fixture the engine's gamma_3 reproduces the same quantity computed
  independently in R -- the third derivative of the exact log posterior along
  the same conditional-mean curve -- to 8e-4 relative, and the zero-inflated
  Poisson to 5e-4. Held against the two-dimensional quadrature of the same
  posteriors it has the right sign and undershoots, closely where the skewness
  is small (0.86 and 0.93 of the exact value at |skew| ~ 0.11-0.13) and by
  about half where it is moderate (0.299 of an exact 0.530). That last case is
  pinned in the suite because it has a consequence: gamma_3 is a LOWER BOUND
  on the skewness, and there it bands "good" where the exact value bands "ok".
  A coupled Gaussian cell, whose Hessian is constant, reads exactly 0.

  The scalar single-coordinate path is untouched: byte-identical across seven
  fixtures on both the family-enum and the spec entry (identical(), max
  absolute difference exactly 0), verified against a build of the preceding
  commit.

• "coupled_likelihood" is retired from the inner-skew decline vocabulary and
  from .INNER_SKEW_STRUCTURAL -- coupling several processes in one likelihood
  no longer describes anything permanently unscorable. What remains is
  "curvature3_unavailable" (a spec that ships no way to reach a third
  derivative) and "coupled_arm" (a coupled fit for which no cell tensor could
  be built at all). Every decline still returns NaN; one unreadable cell takes
  the whole contraction to NaN rather than silently understating the sum.

• New cpp_cell_coupling_curvature3() exposes the contraction at one cell,
  outside any solve, so a registered spec's tensor can be checked against a
  direct numerical third derivative of its own log-density and the step policy
  measured rather than asserted (tests/testthat/test-cell-curvature3.R).

0.0.138

• The inner Laplace layer now has a likelihood-agnostic reliability number
  (#303). gamma_3 scores the inner Gaussian by expanding the joint
  log density along the Gaussian conditional-mean curve at a probed latent
  index, which needs a per-observation third derivative -- so a coupled
  multi-process likelihood (a ZI mixture, tulpaObs's occu_cover) declines
  permanently and the fit has only the outer k-hat, which scores a different
  layer. inner_pareto_k walks the SAME curve and simply evaluates the joint
  density along it: the inner Gaussian is an importance proposal for the exact
  conditional posterior, and the Pareto-smoothed shape of that ratio scores the
  approximation directly. No likelihood derivative anywhere, so it answers
  wherever a mode was found.

  It runs on the probed subspace, not the field. Importance sampling degrades
  with dimension on its own, so a k-hat over all n_x coordinates would report
  n_x rather than the approximation; one dimension per probed index keeps
  every sampling problem 1-D and makes the number directly comparable to the
  gamma_3 for the same index. The engine returns the draws and the joint log
  density at them (src/inner_laplace_is.h); the Pareto fit is the existing
  shared .nested_is_pareto_k() core, which now accepts an injected draw
  matrix, so there is one importance-sampling k-hat in the package rather than
  two. The conditional-curve solve v_i = Sigma e_i is extracted to
  src/inner_laplace_probe.h and shared with the cubic term; neither
  refactorizes.

  A Pareto shape index is scale-free -- it describes the SHAPE of the
  importance-weight tail and says nothing about its size. Measured on the
  engine's own fixtures at 256 draws: a gaussian-family coefficient, where the
  inner Laplace is EXACT and gamma_3 is exactly 0, reads k-hat 0.19 / 0.26 at
  importance efficiency 1.000, and a balanced binomial intercept (N = 500,
  S = 230, gamma_3 = -0.007) reads 0.640 at efficiency 0.99998. Both are noise
  on a proposal that needs no correction. The k-hat is therefore banded only on
  probed indices whose realized efficiency falls below
  .NL_DIAG$inner_k_material_ess (0.995); the raw shape is reported either way,
  and inner_pareto_k_uniform records that no index carried a correction worth
  describing.

  The two inner scores agree where both compute. Across a binomial-intercept
  skewness ladder ((N, S) = (500, 230), (500, 60), (100, 3), (20, 2), (15, 1)),
  |gamma_3| runs 0.007 to 0.897 and the importance efficiency falls
  monotonically with it (0.99998, 0.9962, 0.850, 0.807, 0.634 -- Spearman 1.00);
  the tail shape follows at Spearman 0.90, and the band verdicts agree rung by
  rung. On the coupled fixture, where gamma_3 is NaN for every index, the arm
  with the larger exact posterior skewness (0.53 vs 0.13 by direct quadrature)
  is the arm with the lower efficiency (0.983 vs 0.997).

  .tulpa_combined_reliability() folds the inner layer's two scores into one
  band -- the worse of them where both computed, the one that did where only one
  did -- so a fully coupled fit reads "reliable (both layers good)" instead of
  "inner Laplace not assessed". Reported through diagnostics(),
  print.laplace_diagnostics() and diagnostic_summary(); declines carry a
  reason from the same closed vocabulary the outer k-hat uses.

  The draws are engine-owned and deterministic rather than taken from R's
  stream, so requesting the diagnostic leaves a fit bit-for-bit unchanged and
  the reported k-hat does not flap with the seed. Cost is one joint-density
  evaluation per draw per probed index -- no factorization -- which is why the
  budget is a fixed engine constant rather than the outer diagnostic's
  k_samples, whose draws each cost a full inner Laplace solve.

• Fixed: .tulpa_inner_k_reliability() reads its fields with [[. On a
  declined fit the only field carrying the inner_pareto_k prefix is the reason
  string, which $ would partial-match into the k-hat.

0.0.137

• The engine can now test its own coupled likelihood paths (#300).
  CellCouplingSpec has been virtual-dispatched per cell since the joint driver
  gained a coupled branch, but every genuinely non-separable implementation lived
  downstream in tulpaObs, so the cross-arm scatter, the dense-pair allocation and
  the per-cell derivative contract were only ever exercised by a consumer. A
  minimal coupled likelihood is now registered here as a test fixture:
  test_occupancy_mixture (src/test_cell_coupling_occupancy_mixture.h), a
  two-arm occupancy mixture whose cell density is
  psi prod_v Bern(y_v | p_v) + (1 - psi) 1{no detection}. A cell with a
  detection factorises; a cell with none puts the occupancy state and every visit
  inside one logarithm, so d^2 log p_cell / d eta_occ d eta_det and the
  cross-visit second derivatives are nonzero. It writes both dense cross blocks
  ((occ, det) and the (det, det) self block) rather than taking the rank-1
  self-cross shortcut, so a third-derivative tensor has an explicit Hessian to
  difference, and it declares those two through dense_cross_pairs() while
  omitting the one-row occupancy self block.

• cpp_cell_coupling_evaluate() exposes what a spec actually writes. The
  inner Newton chains each spec's eta-space derivatives through the design and
  scatters them immediately, so nothing a spec computes was visible from R. This
  export drives any registered spec at one cell and returns the cell log-density,
  the per-arm gradient, the per-arm negative-Hessian diagonal and every dense
  cross block, with the same buffer-allocation policy the kernel applies (pairs
  read from the spec's own dense_cross_pairs(), rank-1 descriptor supplied).
  It is the surface a finite-difference check of a spec's analytic derivatives
  runs on.

• The exact-quadrature ground truth reaches the coupled case.
  test-inner-skew.R held the separable scalar reference: integrate the exact
  posterior on a grid and hold gamma_3 against its central moments. The same
  construction is now carried to two dimensions over the coupled fixture's
  intercept-only conditional posterior, with three things asserted rather than
  assumed -- the two-dimensional quadrature reproduces the trusted scalar
  reference on a product posterior, the R density agrees cell by cell with what
  the compiled spec evaluates, and the grid is converged under widening and
  refinement. The fixture's exact marginal skewness is 0.53 on the occupancy
  intercept and -0.13 on the detection intercept, so a coupled cubic term
  (#301) checked against it has something to be wrong about. The
  joint kernel's current behaviour on it is pinned alongside: every probed index
  returns NaN with the reason "coupled_arm", never a silently-wrong 0.

• New tests: tests/testthat/test-cell-coupling-occupancy-mixture.R (the
  per-cell contract at tier 1 -- value against the closed form, gradient against
  a difference of the value, the full cross-arm Hessian against a difference of
  the gradient, the coupled/factorising branch split, the declared dense pairs,
  the grad-only path; then at tier 2 an end-to-end joint fit landing on the exact
  mode of the posterior it claims to solve, a spatial ICAR fit whose cross-arm
  curvature is measured nonzero at its own fitted mode, and dense-versus-sparse
  agreement) and four blocks in tests/testthat/test-inner-skew.R. Shared
  scaffolding is in tests/testthat/helper-coupled-fixture.R.