NEWS.md

tulpa 0.7.5

tulpa 0.7.4

tulpa 0.7.3

tulpa 0.7.2

tulpa 0.7.1

tulpa 0.7.0

This release collects 0.6.11 and 0.6.12; their entries below carry the detail.

On that Calluna fit (LiSC, 32 threads) the engine defaults take 25.3 min and 20 full solves in the final grid besides the 79-evaluation mode-find; the engine behind the 78 EVA fits took 198.85 min on a 22-cell grid for an Ambrosia fit of the same design. Every hyperparameter and fixed effect agrees with the unscreened grid to within 0.012 posterior SD.

tulpa 0.6.12

Behaviour changes

tulpa 0.6.11

Performance

tulpa 0.6.10

Performance

tulpa 0.6.9

Performance

tulpa 0.6.8

Performance

tulpa 0.6.7

Bug fixes

tulpa 0.6.6

Bug fixes

tulpa 0.6.5

Bug fixes

tulpa 0.6.4

WBIC

VI stops at a stationary point

tulpa 0.6.3

A large outer grid warns instead of being refused

tulpa 0.6.2

A declined adaptive lattice says why

The grid cap counts the dispersion cells crossed on top

tulpa 0.6.1

A one-level categorical predictor is named in the error

The grid-cap refusal names the axes that produced the count

integration = "grid_adaptive" stays on the lattice at four or more axes

Wrong answers returned with no warning

Accessors and labels

Input validation and messages

Performance

Documentation

family = Gamma() fits the inverse link it names

warmup >= n_iter is refused before any sampler runs

mode = "laplace" reads its control knobs

Laplace and EB predictive draws carry the random effects jointly

The default fixed-effect prior reaches the Laplace-family backends

VarCorr() reports the sigma_re a conditioning fit actually used

ranef() on a sampler fit reports the random effects, not their z

VarCorr() reports the correlation a sampler fit sampled

tulpa 0.6.0

The hyperparameter share reaches both reporting doors

The fixed-effect interval is gated with a spatial field present

A nested fit reports how much of its marginal the hyperparameter carried

A flagged inner layer names the correction it did not run

The subspace-debias coupling closure works on a grid fit

A collapsed outer grid carries its reading, not only its code

tulpa 0.5.4

The hyperparameter copula holds on coarse grids and curved posteriors

tulpa 0.5.3

Hyperparameter draws keep the posterior's correlation between axes

tulpa 0.5.2

A refined axis reads each cell's box from its own row

tulpa 0.5.1

The copy axis is declared at nine slab nodes, and the consistency pass bisects where the mass is

A random-effect term's prior is an exported header

tulpa 0.5.0

A copy amplitude was reported, and drawn, below zero

So the level was given a cell, half a node spacing wide, reaching below the axis's own support, and two things followed from that one geometry: the level's posterior mass was reproduced by no draw, and draws left the support. On a 25-cell ICAR joint fit whose alpha nodes are 0, 0.25, 0.5 and 1, with 0.460 of the posterior on the level, tulpa_hyper_draws() put 0.000 of 50000 draws on it and 22.8% of them below zero, down to -0.125; the fit's own reported 2.5% bound for the amplitude was -0.111. Each now reads 0.461, 0, and 0.

The split is one declaration, so every read takes it: the reported interval (theta_ci_lo / theta_median / theta_ci_hi) under either within-cell construction, the moment-matched interval a CCD design is summarized with, tulpa_hyper_draws() and the "theta" attribute of tulpa_posterior_draws(), and the per-axis hyperparameter summaries. The copy axis's reported resolution (outer_grid_cell_width / outer_grid_h_over_sd) is now its continuum's, measured in the coordinate the continuum is laid out in; the level owns no cell and so contributes no width. Downstream, a consumer deriving a field SD as alpha * sigma no longer receives a negative one.

A copy fit whose alpha grid carries no zero node also moves, by the same correction and in the third decimal: its outermost cells are now closed half a node step out in log, the coordinate the axis is spaced and integrated in, rather than in the value. Nothing inside the outermost nodes changes, and no axis that was already reported on its own coordinate moves at all.

The batched joint driver returned 15 of the 17 fields it promises

A default axis had only the user's pin to declare itself with

Documentation and packaging

tulpa 0.4.12

A GP lengthscale started five times the spread of its own data

Only the centered parameterization evaluates that log-determinant in the target, which is why the non-centered default never showed it: there the Cholesky enters only through f = L z. Both arms of both doors now pass the check. init_gp_lengthscales() is the one starter, shared with the GP TVC gcol33/tulpa#847 added it for.

Restricted spatial regression on a continuous field

Measured on a confounded fixture (150 points, binomial, 5 seeds): the reported field is orthogonal to the restricted design to 6.2e-15 of its own scale, where an unrestricted fit's is 14.3.

tulpa 0.4.11

A coefficient can evolve as a continuous-time GP

Measured on 24 irregularly-spaced instants (exponential gaps, poisson, n = 288) against a smooth true trajectory: the GP recovers it at cor 0.999 / RMSE 0.042, the rw1 that assumes the spacing away at 0.989 / 0.114. temporal_tvc() gains cov / nu / period / scale_coords, and sigma_prior_U / sigma_prior_alpha are read on this structure too -- the same anchor pair, whether the field samples a log-precision or a log-variance.

A GP TVC lengthscale started five times the spread of its own data

A GP TVC now starts at 0.2 * sd(time) (init_tvc_gp_lengthscale()), and all four kernels pass the check. temporal_gp(parameterization = "centered") starts at the midpoint still and reproduces those deviations to the last digit, which is what attributes them; it is gcol33/tulpa#851.

tulpa 0.4.10

A joint fit integrated a dispersion it could not predict with

cpp_nested_laplace_joint_multi() now attaches fitted_eta by replaying the driver's own eta accumulator at each cell's stored mode, so the arm's offset, every block kind and each cell's block scaling are the ones the inner solve used. A cell the cheap screen pruned reads NA rather than a predictor evaluated at the zero row the grid runner left it -- the predicate is the cell's own log_marginal, not the mode row's contents. fitted_eta_var, the within-cell spread a replicate is drawn with, follows control$fitted_var (new on this front door, TRUE by default) rather than being pinned off, and is requested only at one arm, where the mean it belongs to is stored.

tulpa 0.4.9

An SPDE field beside a random intercept pinned the RE SD at 1

The generic nested driver already carries an spde block type, so the field now goes there as a block beside the RE term's own iid block and both SDs sit on one outer grid. theta_grid gains a third axis, and tulpa_hyper_draws() reports a posterior for it. Measured over sd_re in {0.2, 0.8, 2.0} x 8 seeds, the integrated posterior mean comes out 0.268 / 1.175 / 2.727 -- monotone in the truth and the right order of magnitude, against the flat 1 the conditioned path used whatever the data said. The fixed-effect slope barely moves (mean |error| 0.063 vs 0.076 at sd_re = 0.2, 0.072 vs 0.092 at 2.0, inside a seed-to-seed spread of ~0.07), so this is about the variance component being estimated at all.

Two constructors took arguments that no backend could fit

tulpa 0.4.8

Fractional-nu SPDE took minutes where integer nu takes seconds

The search is now compiled (src/brasil.h), 49x faster per call, with R/brasil.R kept as the reference oracle and pinned to it by test-brasil-cpp-oracle.R. n = 120, the issue's headline case: nu = 1.5 647 s -> 45.9 s, nu = 0.5 600 s -> 51.0 s; n = 40 nu = 1.5 321 s -> 32.2 s. The fitted coefficients are unchanged (n = 40, nu = 1.5: 0.313, 0.630, the values the issue logged), and the two implementations agree on the roots to 1.5e-07 relative across orders 1-4, five beta, four spectrum ratios.

The VI stopping rule was scaled by a quantity with no scale

The rule now compares the ELBO gain ACROSS the patience window -- the mean of its second half minus the mean of its first -- against tol_rel_elbo * (best - worst), the span the run has covered. Both sides are differences, so the additive constant cancels; test-vi-stopping-rule.R pins that the stop does not move when a whole run is shifted by up to 1e5 or rescaled by 100x.

A tier-2 test read a draw column by a name the fit stopped using

The R family registry had no test tying it to the compiled kernels

tulpa 0.4.7

tulpa_posterior_draws() sampled hyperparameters as bare grid-node atoms

tulpa_hyper_draws(fit, cells) is the new door, and tulpa_posterior_draws() attaches it as attr(., "theta") so a consumer reading the attribute is correct by default. It is not a new construction. Each within-cell read the engine ships already defines a density on an axis and each is a per-cell mixture with a closed-form component, so what a draw samples is that component conditional on its own cell: Uniform(e_c, e_{c+1}) over the cell's own box under box_uniform, and an equal mixture of Uniform(v_{c-1}, v_c) / Uniform(v_c, v_{c+1}) under chord. Summed over cells each reproduces the read it came from exactly, so the draws' own quantiles reproduce the fit's theta_ci_lo / theta_median / theta_ci_hi -- measured to 5e-3 at 4e5 draws on a two-axis fixture and to 2e-3 on an end-to-end BYM2 fit, against numbers produced by CDF inversion in a separate code path. Neither conditional reads the weights, so the continuization cannot disagree with the sampler that chose the cell.

The attachment runs with the RNG restored: the latent draws are formed before it and are bit-for-bit what they were, and an existing set.seed() script's downstream numbers do not shift because an extra attribute appeared. An explicit tulpa_hyper_draws() call consumes the stream, as a draw should.

An axis whose support does not admit the box read falls back to chord and one that reaches no CDF at all (a CCD moment rule) keeps its node value, both recorded per axis in attr(., "within_cell") / attr(., "within_cell_declined") from the vocabulary .nl_summary_quantile_read() already uses.

A truth-correlated PIT is not on its own evidence of this defect, which #823's own reading assumed. The same measurement with an exact-posterior arm gives Spearman cor(truth, pit) of +0.490 exact against +0.472 atom at n = 5, +0.281 against +0.286 at n = 20, and -0.029 against -0.036 at n = 200: the correlation is a property of how informative the data are about the hyperparameter, and a weakly identified one produces it under the EXACT posterior too. What the atom provably breaks is uniformity. Any re-read of the occu_cover SBC evidence needs a reference arm before a quartile-bucketed rate is scored against a flat 0.10.

tulpa 0.4.6

Four audit findings: a continuous field + RE combination auto could not carry, the field missing from every non-nested spatial fit's linear predictor, an inline field fit's response stored as a list, and a singular MCAR fixed-effect covariance

tulpa 0.4.4

Five audit findings: dropped offsets, EP quadrature, tweedie accessors, SPDE NUTS

tulpa 0.4.3

Polya-Gamma Gibbs routes sample their stated posterior

A free-covariance block's prior design is declared, not inferred

The HSGP-ST interaction refuses the non-centred flag

Output accessors follow their contracts

An output-contract sweep over every tulpa() backend (34 fits) found the accessors below returning shapes, names or values their generics do not specify. Re-run on this release: 82 checks that failed now pass, none that passed now fails, and every remaining refusal names the fit class and reason.

The spatiotemporal interaction reads its AR1 time margin

Refinement slices carry the whole default hyperprior

Per-thread workspaces are built in place

tulpa 0.4.2

A refined axis's SD reads each cell over its own row's box

The default range prior's dimension is the one the coordinates span

Tests

tulpa 0.4.1

tulpa() states the outer prior, under one name at every door

Fixes

Open at this release: #744 (the outer-grid measurement fixture fits a residual variance four times smaller than it simulates).

tulpa 0.4.0

The nested-Laplace doors state their outer prior

A grid-batched species fit is the fit its own call returns

The batched joint driver weights its grid the way the multi-block driver does

The slow tier re-read under the proper default hyperpriors

Free-covariance blocks and two dispersions carry a proper prior

The cheap screen ranks what the grid integrates

A refined outer grid is measured cell by cell

WAIC, LOO and posterior_predict() read the whole linear predictor

logLik() on a nested fit is the evidence of its outer grid

Every outer axis carries a proper hyperprior, and the evidence reads it

logLik() names its quantity from what the fit estimated

The nested-Laplace route carries offset()

Vignettes compare fits on the evidence logLik() reports

A field-SD axis's decline survives a dispersion placement

Progress knobs are read by exact key

The control surface: one check per door, one home per default

One log-marginal convention, and one entry bundle

Counters that counted the wrong thing, and a rescue that erased what it rescued

Dead code that could still be revived wrong

Reporting and comments

Front-door arguments the door documented and did not honour

Spec constructors: one check, one anchor

Build

What made the failure look machine-specific was the repo's own .Rprofile, which sets options(pkg.build_extra_flags = FALSE) and so keeps load_all() on R's -O2. R reads a project .Rprofile only when the session STARTS in the package root and startup files are not skipped, so the same clone with the same toolchain built from an interactive session in the root and failed from a driver launched elsewhere or under --vanilla.

The outer k-hat's tail is resolved once per fit, and every backend records it

mode = "auto" never picks a backend that errors on the call that selected it

A reported diagnostic says what it is, and what it could not establish

Accessors that returned the wrong shape, or the wrong quantity

Reachability: a family, a backend argument, and a shipped C callable

One covariance code names one kernel, on every path

A periodic temporal GP's period is stated in the units the kernel sees

The two HSGP-ST prior defects the Kronecker branch's fixes did not reach

Every per-observation input passes the finite guard

No NUTS entry can start a trajectory on an unverified gradient

A covariate transform is carried through the formula strippers unchanged

The chain checkpoint fingerprint covers the data, not only its shape

The temporal node index is checked where the areal one is

One spatial gate at the front door, whichever way the spec arrives

tulpa 0.3.2

A placement rescue no longer deletes the checkpoint its own fit wrote

tulpa 0.3.1

A per-arm dispersion axis is placed like a prior block's scale axis

.joint_phi_grid_rescue() now moves it. Everything the placement needs was already there -- the transform registry tags a phi_<arm> column "log", and the FD stencil's re-evaluation already varies it, since .joint_grids_from_cells() hands grids$phi_<arm> straight back to .joint_phi_grid_per_arm() -- so what is new is a slot source, a provenance read and a write target, not a second placement machine. It fires on the axis's OWN sizing (railed, or h / sd past .NL_RECENTER$resolve_mult) rather than on the whole grid's collapsed_edge regime, because a dispersion axis is crossed onto the tensor independently of the field's geometry.

Measured on a BYM2 + Gaussian copy-arm fixture (spatial hyperparameters pinned at truth, four nodes over residual SD 0.02 to 2, true phi 0.09): pinned, h / sd 98.2 and a posterior mean of 0.1745; marked, a posterior mean of 0.0898. The placement is the estimate, not only the report.

Reporting

tulpa 0.3.0

Closes every open engine issue. 0.2.14 closed seven of eleven; the remaining four are below.

The CCD curvature-reuse path is measured, and removed

One dispersion convention, in the fixtures too

One axis marginal, three measures

tulpa 0.2.14

The outer grid says what it did

Placing a design is now budgeted, and screening can no longer hide a collapse

?tulpa_nested_laplace claimed a silently-wrong pruned posterior was impossible. It was not, and that sentence is gone: the gate bounds a mis-ranking among the cells that were solved, and a cell it discarded and never solved is outside what it can see.

A bounded axis has a domain, and a stated axis is a bound

This moves outer weights on any fit carrying a correlation axis whose top cell reached past 1: rho_car's outermost prior weight falls 0.0606 to 0.0388, BYM2 rho's 0.00944 to 0.00527. Point estimates and modes are unaffected; the grid-marginalised posterior moves. An evenly spaced grid whose top node crowds the boundary no longer gets equal weights, which is deliberate -- its top cell cannot extend past 1.

Also fixed here: rho names four different parameters (BYM2 mixing, proper-CAR, AR1, and the multi-output cross-field correlation), so the name-keyed c(0, 1) domain was a false claim that made the shipped hsgp_mo default grid unconstructible. Axis bounds now resolve through the one registry.

This narrows some existing fits. Any fit that stated an axis and relied on adaptive_grid = TRUE extending it now integrates the stated range. It also reaches the var-of-means consistency pass, which runs by default on every joint fit and not only under adaptive_grid = TRUE; and since the engine has no default for phi_grid, every phi_* axis counts as stated. Restore the old behaviour per axis with axis_refine = c(alpha = "extend").

Hyperpriors reach every block they name

This changes results for a direct tulpa_nested_laplace_joint() caller who passed one spec to a grid carrying that axis on several blocks. Broadcast rather than refusal is what keeps the engine's own .joint_multi_sigma_grid_rescue() coherent: it sets a scalar prior_sigma on a multi-block grid and refits, and fires exactly on the multi-copy-block collapsed-edge fits this issue is about.

Fixed in passing

tulpa 0.2.13

tulpa 0.2.12

tulpa 0.2.11

The engine is unchanged across this. Built clean at 66e5f60 and at the fixed fixtures and run on the same outer grid, every cell whose coordinates match exactly returns a bit-identical log_marginal; the cells that differ are the adaptive-refinement nodes, whose coordinates differ because the phi axis does.

tulpa 0.2.10

tulpa 0.2.9

Its two comparators had meanwhile stayed on an inline outside = "clamp", so the comparison varied the collapse construction and the edge policy together when it exists to isolate the first. All three arms now go through the one dispatch at the one within-cell construction. Scored on a common policy the gate's verdict holds either way and the margin is wider under the shipped one (lo: mixed 3.14443 against 4.79910 and 4.48729, where under clamp it is 2.69889 against 3.92249 and 3.83155), so what moved is the margin and not the finding. Bisected across 205 commits, and the pins are re-measured rather than regenerated: the comment records both the flip and the alignment.

tulpa 0.2.8

tulpa 0.2.7

tulpa 0.2.6

tulpa 0.2.5

tulpa 0.2.4

tulpa 0.2.3

tulpa 0.2.2

tulpa 0.2.1

tulpa 0.2.0

Measured, paired arms on the poisson SVC fixture differing only in the constraint, 4 chains x 400 iterations: the penalty was NOT misbehaving on the centered path. beta_x 1.661 (penalty) against 1.654 (centring) at n = 80 and 1.632 against 1.598 at n = 120, truths 1.650 and 1.647; Rhat 1.04 / 1.02 and 1.03 / 1.04; no divergences either way. The change is one construction for one alias, not a repair of a measured defect.

s2z_centre_blocks() (tulpa/sum_to_zero.h) is the single centring the four sites that had hand-written the mean subtraction now share, and that header carries the proper-vs-intrinsic derivation, including why no constant makes a ridge right: matched to the field's own prior on the sum it is a term the field prior already carries, and anything else is a second, unstated prior on the level. The unreachable plain-double compute_svc_eta and svc_sum_to_zero_penalty in hmc_svc.h are deleted.

tulpa 0.1.23

tulpa 0.1.22

tulpa 0.1.21

tulpa 0.1.20

The same branch has a second side: a payload declaring one of these families but carrying no bound fell past the callbacks' bound branches into the ladder, which does not carry them and raises from inside an OpenMP reduction -- std::terminate rather than an R error. prepare() now checks the bound is present on the calling thread, beside the tweedie phi2 check that is there for the same reason. Tests: test-bound-family-joint-arm.R, covering both sides of the branch.

tulpa 0.1.18

The convention this follows is stated once above linkinv: a bounded link's tail is built, never recovered. logit branches on the sign of eta, cloglog uses -expm1 (#454), cauchit uses the arctangent identity. The stable form differs per link, so there is nothing to factor out -- the rule is what is shared. Tests: test-family-link.R, reading the ladder against pcauchy on relative tolerance and pinning that the cancelling form it replaces is wrong by orders of magnitude more, and worse the further out it is read.

check_arg_length is the one place either entry raises that error, so the two no longer carry their own copies of the message. Tests: test-laplace-arg-lengths.R, which also pins that each argument is still accepted at its own length and that supplying all of them reproduces the base fit.

tulpa 0.1.17

Collecting the bundle is one argument-free macro rather than a per-site aggregate initializer. The package builds at C++17, where an aggregate initializer of sixteen members -- several adjacent ints, several adjacent doubles -- can be mis-ordered at one of eleven sites without the compiler seeing it, and the resulting fit would run at the wrong iteration budget or the wrong tolerance. Every entry declares these arguments under identical names, so one member-by-member token sequence binds them at all eleven sites and cannot bind max_iter to n_threads.

Behaviour is unchanged: the entries' existing test files are the gate, and test-nl-entry-forwarding.R is the new arbiter for the failure the duplication produced. It drives all eleven entries and asserts each shared argument is observable in the fit -- max_iter caps the iteration count, x_init moves a capped solve's mode, store_Q / compute_skew / debias / cila each add their fields, and the twelve that enter the checkpoint fingerprint each make a resume against a written checkpoint refuse. A build with one bundle member dropped fails it at all eleven entries.

tulpa 0.1.16

tulpa 0.1.15

tulpa 0.1.14

The checkpoint payload carries pd_conditioned, so CheckpointLog::MAGIC is bumped: a file written by the previous layout would be replayed field by field into the new one and mis-parsed.

tulpa 0.1.13

tulpa 0.1.12

tulpa 0.1.11

The fixture is now ogd_fixture_sim() / ogd_fixture_fit() in helper-outer-grid-dump.R, with within_cell an argument defaulting to the shipped read, and test-outer-grid-dump.R asserts that default still equals .NL_DIAG$within_cell -- so the next flip fails there, naming the files to re-measure, rather than silently changing what they score.

Re-measured under the shipped read, the finding that moves is which PART of the read a rule is shown to reach. The two constructions place the same mass in the same cells and differ by half a cell in where inside one they place it, which on these coarse grids is the scale the location is resolved at, so the floor narrows on the widths (0.2158 -> 0.0882 at five levels) and widens on the median (0.0024 -> 0.0184). The box-mass rule's reach moves from the location to the spread; the barycentre placement's five-level median margin falls from 22x the floor to 2.0x while its widths margin rises to 3.1x. The rules are unchanged -- the box multiplier's per-cell values are identical under both reads. Write-up: dev_notes/issue599/RESULTS599.md.

tulpa 0.1.10

The penalty now carries a third family, each site's squared linear trend, at the precision the s2z contract implies: s2z_precision(n) holds the field MEAN -- the coefficient of the constant direction -- at sd = kappa, and s2z_precision_weighted(u'u) holds the coefficient of an arbitrary direction u' phi / u' u at the same sd, reducing to s2z_precision(n) exactly at u = 1. One predicate, st_needs_trend_pin() (src/st_null_space.h), answers where it applies, and the density, the sparse matrix form and the precision-informed mass override all read it. TYPE_II takes the same term by the same derivation; TYPE_III's kernel is already spanned by the column sums and TYPE_I is proper. A CYCLIC RW2 gets nothing -- a ramp is not periodic, so rw2_rank reports T - 1 there -- and an RW1 fit is unchanged bit for bit.

Measured against the numerical Hessian of the engine's own log posterior, cond(M^-1 Q) under mass_matrix = "gmrf_margin" grew LINEARLY in tau without bound (13.6 / 18.6 / 75.2 / 529 / 3900 over log_tau 0 to 6) where RW1 saturated at 17.4; it now saturates at 25.4. On the #598 harness the RW2 fixture's divergences fall from 43.5 to 8.4 per fit and, per EFFECTIVE sample, from 5.48 to 0.627 against the adapted dense metric's 0.691 -- closing the one finding gcol33/tulpa#598 could not account for. What remains is the tau funnel every arm shares. AUTO still does not select gmrf_margin.

The low-rank mass storage (hmc_mass_lowrank.h) is now generic over WEIGHTED group sums rather than indicator groups, since a trend group covers the same coordinates a row group does and differs only in its weights. Consumer-package paths only: nothing in tulpa sets ModelData::has_spatiotemporal.

tulpa 0.1.9

tulpa 0.1.8

G1 diag(1/C0) G1 carries neither kappa nor tau, so init() builds it once and each proposal costs a scaled add in place of two sparse products. The rational path collapses the same way: its per-pole terms differ in two scalars only, sum_k w_k (kappa^2 + r_k)^2 and sum_k w_k (kappa^2 + r_k). Both hyper derivatives become exact rather than resting on C0 D = I: d(kappa^4)/dlog_kappa = 4 kappa^4 and d(2 kappa^2)/dlog_kappa = 4 kappa^2 give dQ/dlog_kappa = 4 kappa^2 tau^2 K termwise, which is the closed form backward() and forward_with_tangent() already applied.

tulpa 0.1.7

The storage (src/hmc_mass_lowrank.h) is generic over GROUP-SUM directions rather than over the Type-IV margins, so a single sum-to-zero on an ICAR / RW1 / RW2 block is the one-group case of the same term. The inverse is Woodbury on a k x k inner matrix, keeping the per-step cost O(n + nnz(U) + k^2) against a dense metric's O(n^2), and the momentum draw is a sum of two independent Gaussians rather than a factorization of M.

Scored on the #585 paired design -- 6 configurations x 8 seeds, arms sharing the data and the chain seed -- at adapt_delta = 0.95: leapfrog steps per effective sample come out at a pooled geometric-mean ratio of 0.040 against the adapted diagonal (46 of 48 pairs, sign test p = 8.4e-12), raw sampling leapfrog steps at 0.053 on 48 of 48 (p = 7.1e-15), with ess_min rising at the same time and max-treedepth saturation -- 26% to 98% of iterations on the diagonal arms -- disappearing entirely. The "gmrf" arm of the same run reproduces #585's null (pooled 1.08, p = 0.67).

Two limits are on the record rather than smoothed over. At adapt_delta = 0.8 the non-centered configuration is unusable (250 divergences per 1000 iterations at an adapted step size of 2.06), which the 0.95 target repairs; and at 0.95 the three small 3x3 fixtures still show divergences the diagonal does not, pois_3x3_T4_rw2 worst at 43.5 per fit, while the two largest configurations are clean in every arm. Whether that is new pathology or newly visible pathology in chains the diagonal metric was keeping stuck is not settled here and is carried as #598. Write-up: dev_notes/issue597/RESULTS597.md.

tulpa 0.1.6

Scored on 96 paired NUTS fits (6 configurations x 8 seeds, the two arms sharing the simulated data and the chain seed), leapfrog steps per effective sample come out at a pooled geometric-mean ratio of 1.03 on the worst parameter, 0.91 on the interaction block and 1.09 on the hyperparameters, with every per-configuration sign test above p = 0.47 and per-pair ratios spanning 0.18 to 9.2. The conditioning says why, without a sampler: cond(Q) is 1.3e5 to 2.1e5 and neither diagonal read moves it (diag(Q^-1) reaches 1.1e5, 1/diag(Q) changes nothing), while deleting the S + T soft sum-to-zero margin directions takes the same matrices to 11.5 to 51. Those directions are 1_S (x) a and b (x) 1_T, so no diagonal metric reaches them. A mass carrying them as a rank-(S+T) term reaches 6.8 to 18.6, which is #597.

tulpa 0.1.5

Behaviour changes, in the order they are most likely to affect a fit.

tulpa 0.1.4

tulpa 0.1.3

tulpa 0.1.2

tulpa 0.1.1

tulpa 0.1.0

First CRAN release. The engine's surface is unchanged from 0.0.198; this version number marks the release rather than a feature change.

Measured on test-spde-ccd.R's analytic fixture across 1e-4 to 5e-2: every step up to 2.5e-2 returns convergence 0 on Linux and Windows alike and agrees bit for bit, while 5e-2 at factr = 1e5 is the one cell that converges on one platform and aborts on the other, at a mode both platforms agree on to 11 significant digits. On the real inner-Laplace marginal the smaller step reaches the same mode in the same 13 evaluations at a lower objective (505.675740 against 505.675808), so the wider step was costing accuracy at the optimum.

The values move to .NL_MODE_FIND (R/settings.R), which records why the consumers differ, and the stats::optim() invocation itself moves to .nl_lbfgsb_mode_find() (R/mode_find.R), so fit_spde_nested_ccd() and fit_st_nested()'s auto-grid rescue share one call site instead of carrying a copy each. tulpa_pathfinder() reads its factr from the same table without sharing the helper: it is unbounded, takes an analytic gradient when the caller supplies one, and carries maxit / pgtol as its own arguments.

Behaviour is unchanged on every path. The spatiotemporal rescue previously inherited optim()'s own ndeps default, and test-mode-find.R pins that the value now written out is that default by running both calls and comparing par, value, convergence and hessian; the same file checks the helper against the call it replaced for each consumer, and test-settings.R lints the sources so a new factr or ndeps literal outside settings.R fails.

Three consequences, all of them measured on the package's own four-axis joint fixture at axes_at(5), where b1.alpha holds 57.7% of its axis mass on its lowest node (2.88x the uniform share). The whole-grid verdict was read off the axes that happened to score, so resolved could come back TRUE with an unscored axis present; the note named the coarsest SCORED axis, which on that fixture is a different, healthy one, and told the reader to add nodes there; and .nl_railed_axes() -- which finds b1.alpha:lower instantly from stored weights -- was attached only inside the registry grid rescue, which does not run on a caller-pinned grid, so the placement the engine leaves alone by construction was the one that never reported it.

Now: .NL_AXIS_SD_REASONS is a closed vocabulary (too_few_nodes, mode_at_edge, coord_not_finite, stencil_degenerate, curvature_not_negative) carried on the NA as an attribute, so every existing is.finite() caller is unaffected; it is stamped through as outer_grid_resolution_declined, per axis, beside the existing theta_within_cell_declined / theta_cell_edge_declined; resolved requires every axis to have scored; the rail report is attached wherever the resolution is; and the note reports unscored axes and railing BEFORE the coarsest scored one. A fit where nothing scored reports that nothing could be, rather than returning NULL as if it carried no resolution at all.

tulpa 0.0.198

tulpa 0.0.197

tulpa 0.0.196

tulpa 0.0.195

At nominal 0.50: the cap is EXACTLY inert at nine nodes (1.5, 3 and 6 give identical coverage to three decimals in 7 of 8 cells); at the shipped five nodes, 3 is nearer nominal than 1.5 at span = 4 (0.450 / 0.415 against 0.370) and ties at span = 2.5; and 6 reaches 0.520 only on a doubled width. Under the shipped sd_clamp_policy = "decline" a lower ceiling also is not free -- dropping to 1.5 abandons the placement on 32-39% of fits. span and n_pts are kept with it: nine nodes moves 50% coverage FURTHER from nominal while costing proportionally more inner solves.

tulpa 0.0.194

gpu_cuda.h is now included from exactly one place, and the NNGP kernels reach it through gpu_backend.h like everything else. Compiling it in needs no CUDA SDK at build time and no GPU at run time: the driver, cuBLAS and cuSOLVER entry points resolve dynamically and every entry returns FALSE when absent, which is what makes "use CUDA if available" expressible as a default. TULPA_DISABLE_CUDA builds the stubs instead, and that is now a whole-program choice.

cpp_gpu_backend_kind() reports "cuda" or "stub", so which implementation was compiled is observable rather than inferred -- a silent either/or is what let this sit. It is deliberately separate from cpp_gpu_available(), which asks whether a usable device is present at run time.

Two closed investigations should be re-read in this light: gcol33/tulpa#283 fixed a column-major/row-major bug in this path, and gcol33/tulpa#389 ruled the GPU dispatch out as a cause -- a conclusion consistent with the stub having won the link. The per-matrix batch verification added in 0.0.192 (gcol33/tulpa#392) is what now guards the path on machines where it does run, and the gcol33/tulpa#389 determinism arbiter passes with the CUDA backend live on a device-equipped machine.

tulpa 0.0.193

tulpa 0.0.192

Every matrix is now checked against the factorization's own defining identity on the diagonal, sum_k L[j][k]^2 == C[j][j], plus a positive finite pivot. That is O(n_nb^2) per matrix against the O(n_nb^3) the factorization costs, so full coverage sits an order below the work it verifies. The exact one-matrix CPU cross-check stays alongside it as the layout check.

tulpa 0.0.191

The 2-D path is unchanged to the bit: the same neighbour indices and distances, asserted at tolerance = 0 rather than approximately, and a CONSTANT extra column reproduces the 2-D graph exactly.

tulpa 0.0.190

MEASURED before the fix, the same seeds fitted twice in ONE process at n_threads = 1: 10 of 20 fits differed at 49 locations with log_marginal moving 4.09, and 20 of 20 at 50 and 60. Which SIZES broke moved between sessions -- 30 and 60 broke and 120 was clean in a run where the original report had 30 and 60 clean and 120 broken -- so the size-dependence the first reading proposed was the allocation history, not a buffer edge. It was never the batched-Cholesky dispatch either: 49 locations is on the CPU path and broke.

The three loops now share one helper, tulpa_linalg::coords_dist() (src/linalg_fast.h), which sums over EVERY column the coordinate matrix carries. So the coordinate dimension is whatever the caller supplied, and a 1-column matrix is a 1-D domain -- a transect, a depth profile, a time axis -- rather than a 2-D one with a column missing. Three coordinate columns work on the nested-Laplace path for the same reason. The arbiter is that a CONSTANT extra column, which cancels in the distance, now reproduces the 1-D fit bit-for-bit at two different constants.

tulpa 0.0.189

MEASURED over 200 fixed-truth seeds on each of six configurations x two placement policies, arms paired seed by seed and differing only in this setting. Summed |coverage - nominal| over nominal 0.95 / 0.80 / 0.50 at the shipped placement, and the paired 95%-level record against "clamp": "decline" 0.1393 against 0.1464, changing 35 trials, winning 7 and losing none (sign test p = 0.0078) at an unchanged width (ratio 1.0000). A third policy, "relative" -- cap the re-placed span by the INCOMING axis's own span rather than by an absolute bound -- scores 0.1536 and loses 7 trials to 1, and is kept as a selectable arm rather than shipped.

tulpa 0.0.188

The construction has been selectable since 0.0.175 and the default was held back by two things, both of which the placement work in 0.0.187 changed. Until gcol33/tulpa#361 the default axes were laid without reference to the posterior, so every earlier measurement of this choice was taken on a grid pinned coarser than any a user now gets. Re-measured at the placement the engine ships, summed |coverage - nominal| over nominal 0.95 / 0.80 / 0.50, chord against box-uniform: 0.2900 against 0.1233 on gcol33/tulpa#337's own pre-registered fixed-truth instrument (300 seeds), 0.2004 against 0.0361 over the 4680 truth-swept fits of the same fixture whose axis contained the truth, and 0.2467 against 0.1572 over nine (config, axis) rows spanning seven families at 200 seeds each. Box-uniform is nearer nominal on 6 of those 9 rows and at all three levels of the other two, at 0.69 to 1.08x the width.

The conditional-coverage sensitivity that held the default back is a property of ANY within-cell reconstruction, and at the shipped placement it no longer separates the two: box-uniform's realized 95% coverage swings 0.110 across a cell where it swung 0.415 on the coarse pinned grid it was measured on, and at nominal 0.50 the two reads swing 0.238 and 0.231. A resolution-conditional default was scored rather than assumed and is dominated by the fixed rule (best threshold 0.1733 against 0.1572).

The regime it is weakest in is recorded rather than left to be rediscovered: on a coarse grid PINNED by the caller with four crossed blocks and no placement pass, box-uniform wins one resolution and loses the next (per-axis summed deviation 0.5083 against the chord read's 0.6500 at four levels, 0.6875 against 0.5917 at five), because four axes give four independent box positions. That is the regime every measurement before 0.0.187 was taken in, and it is not what a default fit gets.

Two consequences worth knowing. A locally CCD-refined grid is not a partition that tiles, so it declines the new default and reports the chord read with theta_within_cell_declined = "support_mixed" -- a difference a performance knob introduces, which is why it is recorded per axis rather than silent. And .re_cov_derived_summary() is deliberately pinned to the chord read: its values are derived quantities at the nodes rather than the design's own cell coordinates on the reported axis, so a box partition of them would assert a property nothing measured. The RE-covariance backends are unchanged.

tulpa 0.0.187

tulpa 0.0.186

The three things a default-on had to be safe against, all measured on current main after the two changes in 0.0.185.

The flip survives the shipped gate. Scored by paired CRPS against the read a default-OFF fit gives -- the gcol33/tulpa#336 grid mixture -- over 400 prior-predictive replicates each: the rare-event binomial-logit intercept t = -1.895, and the small-group Bernoulli design's two coefficients t = -3.765 and t = -3.201. The pre-gcol33/tulpa#354 read that had no location term scored +3.54 / +6.12 / +4.64 on the same fixtures.

Coverage holds across model classes. Twelve configurations -- the six built-in families on the single-block driver, a rare-event small-group binomial, a small-group Poisson, the same data on the joint driver, and three crossed groupings at outer dimension 3 -- read off ONE solve per seed by the shipped recov_sweep(), so the corrected and uncorrected arms are paired and differ only in the marginal read. Pooled over 960 trials at nominal 0.95: 0.9510 -> 0.9542, standard error 0.0070, with every configuration inside the 3-standard-error acceptance the shipped gates use and the gaussian family identical to the bit. Two small-sample classes move, in opposite directions, and are the whole of the movement: at 200 seeds and three levels the correction takes the small-group Poisson design from 0.8950 / 0.7050 / 0.4200 to 0.9400 / 0.7950 / 0.4650 at nominal 0.95 / 0.80 / 0.50, and the rare-event binomial from 0.9650 / 0.8050 / 0.4900 to 0.9175 / 0.7550 / 0.4700. Summed distance from nominal over the nine cells: 0.295 uncorrected, 0.175 corrected.

The rare-event drop is the exact answer. Coverage at a FIXED truth cannot say whether an interval moved away from the right one -- a credible interval attains its nominal rate averaged over the prior, not at one parameter value. Fixture A's posterior is exact by one-dimensional quadrature, so the same rare-event regime runs at five fixed truths with the exact posterior as a third arm (2000 fits): pooled, exact 0.9470 / 0.8650 / 0.5630, corrected 0.9290 / 0.8650 / 0.5630, Gaussian 0.9625 / 0.8210 / 0.4165. The corrected interval reproduces what the exact posterior does at 0.80 and 0.50 and is 0.018 from it at 0.95, where the Gaussian is 0.044 and 0.147 away at the two lower levels. At beta = -2 and level 0.50 the Gaussian interval contains the truth on 0 of 400 replicates, the exact posterior's on 367, and the corrected one on 367.

The decline paths are no-ops. A coupled fit (every arm a multi-eta unit, so no location term), a coefficient the importance k-hat flags, a coefficient past the shape band and a non-nested fit all report bounds identical to the correction-off fit, to 0.000e+00, while an eligible coefficient on the same fit moves by 0.397. That required gcol33/tulpa#386 in 0.0.185; without it every one of those classes moved.

tulpa 0.0.185

The mechanism is why. m_i = (1/2) sum_j c_j rho_ij and gamma_3(i) = sum_j c_j rho_ij^3 are the same weighted sum at the first and third powers, so a large centre carrying a small gamma_3 is uniformly WEAK correlation -- the well-behaved incidental-parameter regime -- rather than a strong direction being extrapolated. The band was anti-correlated with the pathology it was imagined for.

Measured on those stored coefficient-seeds: 376 rows change and every one is centre_unreliable -> eligible, with the 2620 the shape band and the combined inner band decline untouched. On those 376 the correction they now take scores -0.09665 against the Gaussian they used to report (se 0.00676, t -14.31) where the exact posterior scores -0.09721 -- 99.4% of the achievable gain. A live re-fit of the small-group Poisson fixture agrees through the engine's own record: 21 of 80 rows move, max admitted |m| 2.404, and every one still banded good by gamma_3 and by the combined inner band.

Nothing else changes. cornish_fisher_in_band() is still the one predicate behind both the eligibility record and the quantile path, centre_unreliable is still in .SKEW_CORRECT_REASONS and in the precedence, and the cutoff is still an argument -- .nl_skew_correction_attach(max_abs_centre =) now takes it -- so a finite value restores the band on every path at once, which is what the decline tests drive. gamma1_not_computable is what guards an unformed location term and is untouched.

tulpa 0.0.184

The arithmetic already existed: BuiltinFamilyResponse::weights scales each row's log-density, score and Fisher information alike, and the third-derivative oracle behind gamma_3 reads it too. What was missing was the argument. All six spatial kernels now take it -- cpp_laplace_fit_spatial / _bym2 / _car_proper / _hsgp through build_spec_family_inputs(), and cpp_laplace_fit_gp / _spde / _spde_precomputed through a weights field on JointArm that make_single_arm() fills and the joint driver hands to its built-in-family response. None needed the error fallback. The fractional-nu SPDE marginal (cpp_spde_fractional_logmarginal) is weighted with them, so a weighted mode is never paired with an unweighted log-marginal; its gaussian branch uses the |I + K B| form, which stays finite at a zero weight where the scaled-variance form does not.

Measured rather than asserted. The weighted mode reproduces an R Newton solve on the penalized objective written from its definition to 8.5e-15, against 0.140 for the unweighted one; the weighted score at the reported mode falls from 14.86 to 1.4e-06 (the loop's own step-based floor); H_beta and the Schur at the reported mode now agree to 4.8e-17 with the same w. An integer weight equals row replication on every kernel, which reads the weighted path against the engine's own unweighted solver on a different data set. Nothing moves at weights = NULL: 25 fits spanning every kernel through both the raw C++ entry and the front door (including H_beta) are identical() to the pre-fix build, worst numeric difference 0. No ABI change -- every struct touched lives in src/, none in inst/include/tulpa/.

Evidence in dev_notes/issue385/.

tulpa 0.0.183

tulpa 0.0.182

tulpa 0.0.181

tulpa 0.0.180

tulpa 0.0.179

The convex form (1 - t) y0 + t y1 cannot overflow there, but it is NOT substituted wholesale: the two are not the same double, and over 840000 randomized reads 21.49% of them move between the forms, by up to 4.07e-10 relatively where the difference cancels. Three fixes on this path are pinned on identical() against what approx returns. So .nl_interp_repair() is reached only where the straight read already failed to return a double, and a bracket whose knots are not both finite declines rather than inventing a number. Measured over 1620000 reads (60000 node sets straddling the double range in both signs x five declarations x three outside policies): non-finite bounds 9504 -> 0, finite pre-fix reads that moved 0, repaired bounds outside their own bracket 0, mis-ordered reads 0.

tulpa 0.0.178

Such a mirror now falls back to the extreme coordinates -- representable by construction, and inside whatever support the axis has, since they are coordinates of it -- and records mirrored_edge_not_representable, a fourth entry in the closed .NL_EDGE_DECLINED vocabulary (gcol33/tulpa#293: a silent-disable path needs a reason). It is a second name rather than mirrored_edge_outside_domain because there is no domain here to be outside of, and it takes precedence over a nodes_outside_declared_domain / unknown_domain already in hand: those name a declaration that was set aside, after which the guess ran, and this one names what the edges ARE.

The guess itself is not otherwise narrowed. An undeclared all-positive axis whose log mirror underflows keeps its 0 lower edge, which is gcol33/tulpa#377's boundary from the other side -- with no declaration there is no support to measure an edge against, and car_proper's rho_car keeps exactly the edge it has.

Measured over 595010 reads (59501 randomized node sets spanning the double range in both signs x five declarations x both within-cell constructions): non-finite partition edges 69930 -> 0, non-finite reported bounds attributable to the mirror 251746 -> 0, mis-ordered bounds 0 -> 0, finite bounds outside a declared containing support 0 -> 0, and no read gains a non-finite bound. Byte-identical elsewhere: over 540 (node set, domain, weighting) cases and 4320 quantile comparisons, plus 118405 randomized node sets, everything whose pre-fix partition had finite edges is identical() on the edges, the boxes, the coordinate, the reason and every quantile read; 0 of 118405 moved despite a finite pre-fix edge, and all 960 replaced partitions are finite, bracketing and tile.

Only the linear coordinate reaches it -- log / qlogis / atanh land inside about +/- 745 -- and no default .NL_GRID axis is near it; it is reachable through a pinned grid and through hyper_axis_spec(bounds = c(-Inf, Inf)). test-within-cell-box-uniform.R's boxes_do_not_tile fixture moves for the second time, off the undeclared c(1, 1e300, .Machine$double.xmax) -- which now tiles -- and onto the unit axis at the subnormal floor, whose decline is a property of the node set against its coordinate and which no edge guard changes.

tulpa 0.0.177

Only the LINEAR coordinate reaches it -- log / qlogis / atanh land inside about +/- 745, where no sum of two coordinates can overflow -- so it is the unbounded declaration and the undeclared axis whose values are not all positive. No default .NL_GRID axis is near it; it is reachable through a pinned grid and through hyper_axis_spec(bounds = c(-Inf, Inf)).

a / 2 + b / 2 and NOT a + (b - a) / 2, measured rather than pattern-matched (dev_notes/issue378/midpoint378.R, 4e6 randomized pairs spanning the double range in both signs). The offset form overflows on b - a for opposite-sign extremes -- reachable, since the linear coordinate is signed -- taking 1117 pairs to Inf that a / 2 + b / 2 returns finite, and it rounds twice on ORDINARY operands, moving 52536 of the 3998914 pairs where (a + b) / 2 is finite. a / 2 + b / 2 produced 0 non-finite results, and stays inside [a, b] on 2e7 subnormal-heavy ordered pairs.

BYTE-identical wherever the sum form worked at all. Halving a normal double decrements the exponent and leaves the significand alone, so a / 2 and b / 2 are exact and their sum carries ONE rounding of the same real the sum form rounds once -- rounding to nearest commutes with scaling by a power of two. Asserted with identical() against the pre-fix file sourced from git into its own environment over the namespace (dev_notes/issue378/identity378.R): of 135 (node set, domain) cases over 4 weightings x 4 support kinds x 2 within-cell constructions x 11 probabilities, 133 are identical on the cell edges, the box partition and every quantile read, and the 2 that move were both pre-fix DECLINES that now tile. Over 197374 randomized node sets across all four domains, the fix removes 24045 declines and introduces 0, with 172860 tiled under both and 0 of them moved; all 23913 recovered partitions are finite, strictly increasing, of the right length, contain their own coordinates and stay inside the declared support. The 3621143 pairs where both halvings are exact carry 0 differences; the two forms part only where a halving is subnormal (|x| below about 2.2e-308), by at most 1.0e-320.

The decline is kept, not removed. test-within-cell-box-uniform.R's boxes_do_not_tile pin moves to the same node set UNDECLARED, where the decline is a property of the node set rather than of the arithmetic: the partition mirrors the extreme half-spacing in the value itself and double.xmax + (double.xmax - 1e300) / 2 is past the double range in any form, so there is no better midpoint to compute. A second pin covers a node set the coordinate map cannot separate (a unit axis at the subnormal floor, where plogis(qlogis(1e-320)) underflows to 0), alongside the zero-width interior box already pinned. A new randomized identical() sweep over all four domains holds the bisector to the sum form's own numbers.

tulpa 0.0.176

Measured on the issue's own node set, c(1e-320, 1e-310, 1) declared positive: exp(log(1e-320) - 0.5 (log(1e-310) - log(1e-320))) underflows to exactly 0, and the fit reported a lower bound of 5.93e-323 for a quantity whose declared support is x > 0. From the other end of the double range, c(1, 1e300, .Machine$double.xmax) overflowed to Inf, fell past the log guess to the LINEAR mirror, and reported -4.97e+299 and Inf for the same declared support. It is gcol33/tulpa#369's defect one step down: there a unit axis reported 1.0028 because the domain never reached the edge; here it reaches it and is overruled.

THE PRECEDENCE, not a patch on the symptom. A declared domain that CONTAINS every coordinate is authoritative: its mirrored edge is taken when finite and in-domain, and otherwise the partition DECLINES to the extreme coordinates themselves -- which that same containment test has already placed inside the support, and which are the conservative answer -- recording mirrored_edge_outside_domain. The guess is RESTRICTED, not removed: it runs where there is no declaration to contradict (no domain named, nodes_outside_declared_domain, unknown_domain), and car_proper's rho_car on the adjacency eigenvalue interval keeps exactly the edge it had. A node set the declaration does not contain cannot be honoured by any edge at all -- the edges bracket the coordinates -- so that case is the guess's by construction and says so rather than silently taking it.

The reason field is gcol33/tulpa#293's requirement and rides the hook gcol33/tulpa#357 already built: .nl_cell_partition() reports the COORDINATE it settled on, so it is also where the decline lives, and both readers take it from the same return. It travels out through .nl_summary_quantile_read() and .nl_axis_quantiles() to the fit as theta_cell_edge_coord / theta_cell_edge_declined, one per axis, and diagnostics() says which axis reports a conservative bound and why. .NL_EDGE_DECLINED is the closed vocabulary.

WHAT MOVES AND WHAT DOES NOT, proven with identical() rather than a tolerance (dev_notes/issue377/identity377.R, identity377.out): the pre-fix R/nested_laplace_moments.R is sourced from git into its own environment over the tulpa namespace, and 22 node sets x 5 domain declarations x 4 weightings x 4 support kinds x 2 within-cell constructions are compared call for call. Of 110 (node set, domain) cases, 100 are byte-identical on the edges, on the box partition and on every quantile read; the 10 that move are exactly the ones whose partition reports mirrored_edge_outside_domain, and no such case arises on an undeclared axis. Over 10324 node sets inside their declared domain -- the fixtures plus a 20000-draw randomized sweep across all four domains -- 174 produced an edge outside that domain before the fix and 0 after.

The invariant asserted alongside gcol33/tulpa#369's unit / correlation regressions in test-nl-interval-support.R is the general one: whenever every coordinate is inside a declared domain, both edges are finite and inside it too, for every entry of .NL_DOMAIN_TRANSFORM. test-within-cell-box-uniform.R's boxes_do_not_tile fixture moves to the unbounded declaration, where an interior midpoint's own sum still overflows before it is halved; declared positive, that same node set now has a finite in-support partition and the box read no longer declines on it.

tulpa 0.0.175

THE DEFAULT DOES NOT MOVE, and that is shown rather than asserted: seven fits spanning the single-block registry, the joint single-block driver, the multi-block tensor, CCD and adaptive paths, a bounded bym2 rho axis and the gcol33/tulpa#337 fixture were run against HEAD's R/ and against this tree with the same DLL, data and control lists, and every reported theta_median / theta_ci_lo / theta_ci_hi / theta_mean / theta_sd / log_marginal / weights and per-block moment is identical() (dev_notes/issue357/identity357.R, identity357.out).

It is a second FIELD on each .NL_SUPPORT entry rather than a fifth kind: outside is a fact about what the producer left behind, within is a caller's choice about how to read it, and the two are orthogonal. Only a "density" support admits box-uniform -- a CCD design is a moment rule, a locally refined grid's replacement clouds sit inside one base cell so a Voronoi partition of the node set is not the design's own boxes, and a posterior sample's values are order statistics. .nl_cell_edges() is now .nl_cell_partition()$edges, which also reports the COORDINATE the edges came out of, because the box read has to bisect the interior spacings in the same coordinate the outer half-cells were mirrored in or the partition it tiles the axis with is not one partition. That reconciles the construction with gcol33/tulpa#369: on a unit axis the boxes are bisected in logit, so a BYM2 mixing weight with a node at 0.99972 (gcol33/tulpa#361 extended bym2_rho to six nodes including 0.999) stays inside (0, 1) where the log guess reaches 1.0015. dev_notes/issue337/recon.R's slice_masses() STOPS when the boxes do not tile; a reported interval cannot take that behaviour, so a partition that is not finite and strictly increasing DECLINES with a recorded reason ("support_<kind>", "single_node", "boxes_do_not_tile", "no_usable_node") and that axis falls back to the chord read on its own.

THE PARTITION COMES FROM THE GRID AND THE MASSES FROM THE WEIGHTS, which is what gcol33/tulpa#337's own "keep the masses, tile the axis" says and what recon.R does (axis_geometry() on the whole grid, slice_masses() on the positive weights). A cell whose integration weight underflows to exactly 0 still SITS on the axis, and its coordinate is what fixes its neighbour's box edge; filtering the coordinates by weight -- which is right for the chord read, whose knots ARE the positive-weight coordinates -- shrinks that neighbour's box to nothing. On the coarsest rung of the re-scored ladder (2 cells at 400 groups) the softmax underflows one of two cells on 43 of 150 seeds, so that filter was the difference between measuring the construction and measuring the chord read under another name. An empty interior box is a FLAT segment of the CDF: the quantile is located on the cumulative mass and evaluated inside the box it lands in, so an empty box is stepped over rather than interpolated across.

tulpa 0.0.174

tulpa 0.0.173

tulpa 0.0.172

tulpa 0.0.171

tulpa 0.0.170

tulpa_parallel_sum() (src/omp_threads.h) now cuts the range into team contiguous chunks by index arithmetic, sums each chunk left to right into its own slot, and adds the slots in chunk order. Nothing about the answer is left to the runtime: not the combination order, and not the partition either, so a fit reproduces itself even where the runtime hands back a smaller team than the one requested. The parallel loop runs over the chunks rather than the observations, which keeps the combined parallel for construct. Because gcol33/tulpa#373 had already funnelled every one of these sums through that one helper, the change is one function and reaches the joint and non-joint paths together.

What this buys is reproducibility AT a thread count, not across thread counts: chunking the range imposes its own association, so a fit at n_threads = 4 still agrees with the same fit at one thread only to floating-point tolerance (measured at 6e-14 on log_marginal over four families), which is the invariant ?tulpa_nested_laplace_joint now states.

The one-thread route is untouched: it keeps its plain loop and allocates nothing, and the same 90 fields are identical() before and after the change. The same capture taken twice at n_threads = 4 after it moves 0 of 90. Cost, with the three reduction shapes alternated inside one process so they share the machine's load: 1.00 / 1.03 / 1.02 times the old cost at 48 observations with a team of 2 / 4 / 8, and indistinguishable at 960 and 12000 observations. Rows where all three arms run the same code move by up to 17% in the same capture, so the 3% is at the edge of what the instrument resolves.

tulpa 0.0.169

Sites moved onto the helpers: the single-arm spec eta assembly (compute_eta_spec, which had no guard at all and so entered a one-thread region unconditionally on every non-joint nested-Laplace objective evaluation), both compute_total_log_lik overloads for built-in families, tulpa_linalg::matvec (reached per leapfrog step through precompute_generic_fixed_eta), the five per-observation regions in pg_gibbs_core_step and the ten per-sweep regions across the four Polya-Gamma samplers, and the NNGP non-centered reverse sweep (nngp_nc_backward), whose per-thread workspace slot and accumulators are now parameters of a single-sourced row body rather than captures.

A serialised region body runs in index order, which is what the plain loop does, so both routes agree bit for bit, reduction(+:) included. Verified over 292 captured fields and 68968 doubles -- non-joint nested Laplace at three families and three field types, a corrected-integrated-Laplace fit, cpp_laplace_fit, six Polya-Gamma samplers, generic NUTS, NNGP exact NUTS and single-arm joint fits -- at one thread and at four, and again under OMP_NUM_THREADS=1. Every field is identical() before and after except a handful on the joint fitter at n_threads = 4, which a control run showed is not reproducible run to run under a FIXED build in either direction: a num_threads(4) reduction combines its per-thread partial sums in an unspecified order.

Cost off a corrected-integrated-Laplace M ladder on the non-joint fitter (7-cell sigma grid, one thread, variant = "qmc", two runs per arm): the per-objective-evaluation fixed cost falls 3.39 -> 1.31 us binomial, 3.66 -> 1.34 us poisson and 3.25 -> 1.39 us gaussian, and the per-observation cost falls about 8% with it, since a plain loop keeps the eta assembly in the caller's frame instead of an outlined region body. Per auxiliary draw at N = 48: 9.96 -> 7.27, 8.97 -> 6.58 and 9.91 -> 7.77 us.

tulpa 0.0.168

BuiltinFamilyResponse::prepare() (and FamilyLogLik::prepare()) resolves the family code once into a FamilyKind plus a parsed FamilyLink, and fills a per-observation vector of the constant. log_lik_for_family and grad_hess_for_family gain an entry point taking both, so the objective adds a looked-up double where it called lgamma. Each ladder is written once, in a core that the std::string overload forwards into after an allocation-free classification, so a caller still passing a string does exactly the work it did before and no path regresses. A response that was never prepared reports so and takes the string route, so a construction site that forgets the call loses the saving, never the constant.

The split densities are kernel + const in that order at every call site, so a caller holding a precomputed constant reproduces them bit for bit. Measured over 23 fits (ten families x single-arm and joint, plus three corrected-integrated-Laplace fits) and 253 captured fields -- log-marginals, modes, per-cell modes and Hessians, grid weights, coefficients, intervals, vcov, logLik -- every field identical() before and after, no tolerance.

Cost, read off a corrected-integrated-Laplace M ladder (7-cell sigma grid, one thread, variant = "qmc"), the two arms alternated so they share the machine's load: per-observation cost of the objective 0.190 -> 0.149 us binomial, 0.230 -> 0.138 us gaussian, 0.189 -> 0.118 us poisson. Per auxiliary draw at N = 48 that is 9.9 -> 6.4, 11.3 -> 6.4 and 10.3 -> 5.7 us.

tulpa 0.0.167

tulpa 0.0.166

tulpa 0.0.165

tulpa 0.0.164

MEASURED on gcol33/tulpa#341's fixture B -- a small-group Bernoulli GLMM over a 7-point hyperparameter grid, which is the shape of a real consumer fit -- at 400 prior-predictive replicates scored by tests/testthat/helper-sbc.R against a machine-accurate quadrature. Paired CRPS against that quadrature on the intercept:

| arm | CRPS - exact | t | SBC p_unif | |---|---|---|---| | shipped nested read | +0.0447 | 4.44 | 1.6e-03 | | + subspace debias | +0.0191 | 3.17 | 0.017 | | + cila, M = 1024 | +0.00090 | 1.10 | 0.15 | | + cila, M = 4096 | +0.000033 | 0.05 | 0.80 |

The slope behaves the same way (subspace +0.0138 at t = 2.38, cila at M = 1024 +0.00027 at t = 0.30), and the hyperparameter, already correct before, stays correct. All three variants sit inside k < 0.7 on every one of the 400 replicates, median 0.42 / 0.51 / 0.51 for qmc / is / rqmc.

COST on the same fixture, medians over 8 seeds at 3 repetitions: M = 1024 is 10.2x a plain fit and M = 4096 is 37.6x, against the subspace debias's 26.1x. The correction is one extra pass over the settled grid plus M evaluations of the Newton loop's own penalized objective per cell, and that closure takes one latent vector at a time, so the M evaluations dominate and are what a batched density interface would address.

tulpa 0.0.163

tulpa 0.0.162

tulpa 0.0.161

tulpa 0.0.160

tulpa 0.0.159

The two PD policies now live in one place (src/joint_pd_step.h) and both loops are backends for them: pd_lm_escalate() is the smallest diagonal load making the factorization succeed (Nocedal & Wright, Numerical Optimization 2nd ed., Alg. 3.3), pd_eigen_clamp_solve() the clamped spectrum. With H already PD the first attempt succeeds, no load is added, and the step IS the plain Newton step, so every fit that factorized before is unchanged -- the dense and sparse loops agree to 4e-16 on the mode and exactly on log_marginal for both an indefinite and a PD start. Measured over the issue's own 80-replicate sweep, non-convergence goes from 32.5 / 22.5 / 13.8 / 7.5 / 3.8% (the five affected configurations) to 0% in all ten, each converged fit landing within 1.9e-06 of an independently optimized mode in a median of 4-7 Newton iterations.

Three reports named the wrong layer for it, the gcol33/tulpa#293 shape one layer over. inner_skew_declined said backend_unsupported, which is false -- the same backend computes gamma_3 on every converged solve of the same fixture; grid_fixed_declined said block_not_extracted, a retention step downstream of the real cause. Both now say not_converged, which is also a reason in the shared Pareto-k vocabulary, and the three probe drivers settle it in one place (.inner_skew_attach_probe()) rather than each writing its own tail. $modes is kept -- it is the solver's record of where it stopped, and the warm-start chain and the refinement passes read it -- but the per-cell-mode AVERAGE every coefficient report falls back on now reads only cells that reached a mode, so a stalled fit reports NA with interval_declined = "not_converged" instead of its own start vector as an estimate. A non-PD Hessian at the returned point also withholds the stored precision and the fixed-effect block, whose inverse is not a covariance there.

New test-coupled-indefinite-start.R carries the sparse-detection fixture the suite lacked: every coupled fixture sat at b_occ = 0.2, b_det = -0.5, where the start Hessian is positive definite, which is why nothing reached the path. It establishes that it reaches the regime -- lambda_min < -5 read off the R log posterior, with the old setting asserted PD as the control -- rather than reaching it by luck.

tulpa 0.0.158

tulpa 0.0.157

tulpa 0.0.156

tulpa 0.0.155

tulpa 0.0.154

tulpa 0.0.153

tulpa 0.0.152

tulpa 0.0.151

tulpa 0.0.150

tulpa 0.0.149

tulpa 0.0.148

tulpa 0.0.147

tulpa 0.0.146

tulpa 0.0.145

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".

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 gcol33/tulpa#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 gcol33/tulpa#318 gate above the quantile scores 24.59641 there).

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.

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.

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.

tulpa 0.0.144

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.

tulpa 0.0.143

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.

tulpa 0.0.142

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 gcol33/tulpa#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).

tulpa 0.0.141

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); gcol33/tulpa#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.

tulpa 0.0.140

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.

tulpa 0.0.139

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.

tulpa 0.0.138

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.

tulpa 0.0.137

tulpa 0.0.136

No parallelism is given up for this. The chunks are dispatched as OpenMP tasks, so however many threads are genuinely idle still drain them -- only the partition is pinned, never the number of workers executing it. In the bulk of the grid the budget is 1 exactly as before, so those cells stay serial and allocate no partial buffers.

tulpa 0.0.135

tulpa 0.0.134

Pinning is also PER AXIS rather than all-or-nothing: tau_lower / tau_upper hold the two precision axes (they build both), n_grid_spatial / n_grid_temporal one each, and n_grid_rho / rho_lower / rho_upper the ar1 autocorrelation axis. A pinned axis keeps exactly the nodes its knobs built and is named in the new outer_grid_pinned_axes; the rest are recentred as usual. Only when EVERY axis is pinned does the rescue decline outright.

tulpa 0.0.133

tulpa 0.0.132

Provenance is now explicit, in one predicate (.nl_axis_is_pinned(), R/nested_laplace_auto_grid.R) that all four rescues share -- the joint single-block, the joint multi-block copy block, the standalone registry (icar tau_grid / bym2 sigma_grid) and fit_st_nested(). An axis counts as a DEFAULT (recentre-able) when it is absent, when it is marked with the new auto_grid(), or when its node set is exactly the engine's own default axis for that field -- a grid identical to the default carries no information a pin would add. Anything else is a pin and is never moved.

Contract change: handing the engine's own default axis in explicitly used to hold the grid fixed; it now recentres like the defaulted axis it is (test-nested-laplace-registry-auto-grid.R updated, and that fit's k-hat now agrees with the defaulted-grid fit's on a collapsed grid). Use control$auto_recenter = FALSE to hold a grid where it is.

tulpa 0.0.131

tulpa 0.0.130

The fixed grid is now a starting axis, not a ceiling: a fit that collapses onto a boundary re-centers via the mode-Hessian its own outer Pareto-k diagnostic already computes (R/nested_laplace_auto_grid.R) and refits, reusing that curvature rather than running a second optimizer. An explicit sigma_grid / tau_grid always wins -- auto-recenter only engages when left at its default. Wired through: - the joint driver's single-block backends (bym2/icar/car_proper) and multi-block copy blocks, with a second attempt composing a light default PC(U=3, alpha=0.01) prior on sigma for a genuinely unidentified (near-separation) mode that geometry alone cannot settle; - the standalone .NL_REGISTRY path (icar's tau_grid, bym2's sigma_grid), one recenter attempt, reusing the joint path's generic axis-tagging and FD-Hessian machinery; - fit_spde()'s explicit method = "grid" path, which now attaches the same pareto_k_regime diagnostic (visibility only -- fit_spde()'s default control$method is already "ccd", the mode-Hessian path, so "grid" is a deliberate opt-in the fix respects rather than overrides).

Byte-stable when the mode already sits inside the old default axis (a no-op branch, exercised by regression tests in test-nested-laplace-joint-auto-grid.R and test-nested-laplace-registry-auto-grid.R). car_proper (its rho axis is unguessable, same limitation the existing outer-k-hat diagnostic already has) and MCAR (log-Cholesky axis geometry, a materially different recentering problem) are out of scope; fit_st_nested()'s spatiotemporal grid got the diagnostic only, since it has no existing mode-find machinery to reuse -- tracked as #291.

tulpa 0.0.129

tulpa(y ~ x, temporal = temporal_gp("t"), mode = "hmc") now fits: the spec is validated at the front door, build_sampler_model_inputs() accepts type = "gp" and fills TemporalGPData from the unique time instants, and the field's two hyperparameters are sampled jointly with it. They are named too -- log_sigma2_temporal_gp / logit_phi_temporal_gp rather than param[3] / param[4]. The field is sampler-path only (there is no nested-Laplace kernel laying a grid over a dense T x T Gaussian), and it cannot yet share a fit with a spatial or latent() block; both now say so.

New src/temporal_gp_kernel.h holds the covariance templated over the scalar type, so the sampled (sigma2, phi) can be autodiff variables -- which is why the old plain-double kernels could never have been wired here. The exponential kernel (equivalently Matern nu = 0.5) is an Ornstein-Uhlenbeck process, so it keeps the exact O(T) Markov recursion and its numbers are unchanged; Matern 3/2 and 5/2, Gaussian and periodic have no finite-dimensional state-space form and are evaluated from a dense T x T Cholesky, in both the centered and non-centered parameterizations.

Matern is offered at nu in {0.5, 1.5, 2.5} only -- the smoothnesses with a closed form -- and anything between them is now rejected at construction rather than quietly run as exponential. test-temporal-gp-frontdoor.R asserts the five kernel configurations DISAGREE (a test asserting they agree would have passed before this), that Matern nu = 0.5 reproduces the exponential fit to the bit, and that the periodic kernel tracks its period.

tulpa 0.0.128

The structural fingerprint is the one that punished a mistake quietly: it keys the grid checkpoint, so a copied block that folds the wrong structure produces a checkpoint that MATCHES across runs it should not, and a resumed run then reuses cells computed under different inputs. New tulpa::NlFieldIdentity (nested_laplace_checkpoint.h) names each structural group once -- areal(), nngp(), hsgp(), temporal() -- and each entry point chains the groups it carries:

cpp const std::uint64_t struct_seed = tulpa::NlFieldIdentity("st_icar") .areal(n_spatial_units, adj_row_ptr, adj_col_idx) .temporal(temporal_type, n_times, cyclic, temporal_idx) .seed();

Fold order is part of the fingerprint, so the optional members (the BYM2 mixing scale, the standalone temporal field's group count) sit in the slot they have always occupied and every seed is unchanged. New test-nl-field-identity.R checks that bit for bit against the folds the entry points used to write by hand, for all eleven field models, and pins that each structural input still moves the seed.

The 1-based-to-0-based skew_idx conversion had nineteen copies across src/; all now call the tulpa::unwrap_skew_idx that already existed for it in laplace_spec_fit.h. The five spatiotemporal entries returned their temporal axes through a repeated pair of lines carrying an ar1-only conditional; that is nl_attach_temporal_grids now.

No behaviour change: the entry points shed 225 lines of plumbing, and the fingerprint values and returned lists are identical.

tulpa 0.0.127

Deleted whole: hmc_latent_grad.h (a closed component superseded by hmc_latent.h's apply_first_zero), hmc_tvc_autodiff.h and log_post_car_proper_det.h (never #included anywhere), and hmc_temporal_gp.h. That last one WAS included, by hmc_sampler.h, but every symbol in it was unreachable: its namespace tulpa_temporal_gp is named nowhere but its own opening and closing brace. The live temporal GP is tulpa_priors_temporal.h, which is templated for autodiff and writes the exponential-kernel state-space recursion and the non-centred transform inline; the deleted header was a plain-double implementation that could not serve the gradient modes the sampler uses. Its five temporal_cov_* kernels were a second, untested copy of covariance math the canonical templated tulpa_svc::compute_cov already holds, pinned by test-cov-kernel.R.

Deleted in place: gp_nngp_gradient_w_analytical (hmc_gp_gradients.h), multiscale_gp_log_lik (hmc_gp_log_lik.h), and hmc_tvc.h's dead hyperparameter-prior and gradient block (log_prior_tau_pc, log_prior_rho_uniform, log_prior_rho_beta, the finite-difference rw2_gradient, parse_tvc_structure). The live TVC priors are in tulpa_priors_tvc.h (PC prior on log-tau via pc_prior.h, Uniform(-1, 1) on rho) and the live gradients are the analytic ones in hmc_tvc_grad.h.

SelectedInverse::at was reported as dead and is not: implicit_diff.h:196 calls it as H_inv.at(...), which a search for the qualified name misses. It stays.

tulpa 0.0.126

There is now one implementation of each solve, templated on tulpa_linalg::TriLayout, and the layout is a required template argument rather than something a call site inherits from argument order: tri_solve_lower<TriLayout::RowMajor>(L, n, ld, b, y) and tri_solve_lower_transpose<...>. chol_factor_lower and nngp_moments_from_chol -- the producer and the consumer that must agree with the solve -- carry the same parameter, so a call site states the whole convention it is asserting. chol_log_det does not: the diagonal sits at i * ld + i under both. The layout folds at compile time, so the emitted arithmetic and its summation order are unchanged.

tri_solve_lower's old column-major body had no callers and is gone; the one caller of tri_solve_upper_transpose (the dense mass matrix's momentum draw, which reads an Eigen factor and so really is column-major) now says so. New test-tri-solve-layout.R states the contract -- each routine reads the lower triangle of the matrix its declared layout spells out of the buffer -- and checks it on matched and mismatched buffers, including the cuSOLVER-shaped one whose opposite triangle still holds the input.

tulpa 0.0.125

tulpa 0.0.124

tulpa 0.0.123

tulpa 0.0.122

tulpa 0.0.121

tulpa 0.0.120

Three changes, all on the joint nested-Laplace backend (single- and multi-block paths alike):

The pareto_k band itself is left to speak for the number it reports: the fix is the number, not a verdict override. sn_match() remains the single source of truth for the cumulant inversion; the proposal path adds only a vectorized sampler and log-density (sn_cdf() / sn_quantile() route through Owen's T by quadrature per point, which is right for a few reported quantiles and unusable for hundreds of proposal draws).

tulpa 0.0.119

tulpa 0.0.118

tulpa 0.0.117

tulpa 0.0.116

tulpa 0.0.115

The per-cell cross-Hessian scatter multiplies EVERY active dof of one coupled arm's row -- beta, RE, and any latent-block dof that row's idx / obs_indices resolves to -- against every active dof of another (or the same) coupled arm's row sharing the cell. build_joint_hessian_pattern's cross-arm section only ever registered the beta/RE part of that product; a block reached by just one side (a private random-slope block, or a field the other arm's field_coef = 0 decouples) had no pattern entry for its cross term against the other arm's beta/RE/latent dofs, even though the scatter produces a real nonzero value there whenever the coupling spec's cross-Hessian for that arm pair is nonzero. A new section walks each coupled cell's per-arm active-dof union (mirroring the scatter's own collect_coupled_row_latents resolver) and adds the missing cross entries, scoped per cell so it costs no more than the scatter already does. Two regression tests reproduce both shapes on tulpa's own bivariate test-coupling spec (test-cell-coupling-cross-hess.R), independent of tulpaObs.

tulpa 0.0.114

auto now routes any random-effect term -- intercept-only or slope -- to the exact Metropolis-within-Gibbs covariance debias (re_cov_gibbs), which already treats a scalar (1 | g) as the degenerate c = 1 covariance block, so no special-casing by term shape was needed. An explicit mode = "laplace" / "mala" / ... still conditions on sigma_re (defaulting to 1) when the caller names it directly -- only the auto default changes. The message() that reported the conditioning is now a warning(), so a script that promotes warnings sees it (same sub-issue as #265).

tulpa_eb() and tulpa_re_cov_nested() now share a hyperprior argument, "flat" (default) or "pc_lkj": "flat" matches the nested-Laplace convention everywhere else in the engine; "pc_lkj" opts into the weakly-informative PC + LKJ prior re_cov_pc_lkj_prior() builds, still available on request. tulpa_re_cov_gibbs()'s Sigma | b conjugate draw cannot go fully flat -- an improper prior is not a valid target for that step -- so it keeps its existing minimal-proper Inverse-Wishart default (prior_df = n_coefs + 1), documented as the closest analogue. See vignette("priors").

tulpa 0.0.113

Each (1 | g) now becomes an iid latent block, so its SD is one more axis of the outer grid beside the smoother's tau. The block type, its integration and its recovery already existed (the #86 coupled field + RE capability); only the front door was not using it, and .tulpa_fitter_args()'s own comment already named this as the intended treatment. Consequences:

Recovery is asserted as coverage, not as a point tolerance. The grid integrates the SD under a prior flat in log(sigma), the convention every nested scale axis uses; that does not shrink, so at G = 15 the posterior is wide and its point summary runs above the truth (mean of medians 1.02 against 0.9 over 6 seeds) while the 95% interval covers the truth 6/6. Asserting |est - truth| would encode the +0.12 as the target. Whether the engine should carry one hyperprior convention across the nested and RE-covariance paths is #268, deliberately not bundled here: a PC prior on the iid axis alone would trade the cross-path inconsistency for one inside a single fit.

tulpa 0.0.112

The warning is opt-in per redirect site (notify =), because two different things were being conflated. A smoother sending an eb / agq request to the nested kernels takes away the random-effect SD those two would have estimated, and warns. A random slope under mode = "laplace" has no scalar sigma_re to condition on, a temporal field has no conditional-Laplace kernel, and an SPDE field redirected to the spde backend is the same mode and tier reaching its own integrator -- those are documented routes for the structure, not capabilities taken away, so they are recorded on the fit without warning on every fit.

tulpa 0.0.111

The adaptive Gauss-Hermite inner marginal (control$re_cov = "aghq", or n_quad > 1) integrates each group out by quadrature and so forms no per-group posterior at all. It now says that, with the two modes that do report one, rather than returning the empty frame: ranef() errors on the stated fit$ranef_unavailable reason. ?ranef documents what each backend reports and why.

tulpa 0.0.110

tulpa 0.0.109

tulpa 0.0.108

Nothing looked wrong from R. icar_center_field still centres the field on its way into eta, so the linear predictor stayed invariant and the field's shape and the fixed effects came out right; only the level was loose. On tulpaRatio's 8-unit chain the reported mean(phi_spatial) was +200.9, -342.8 and +423.1 at three consecutive seeds against a legacy level held within 0.02 of zero, and tau_spatial sat at 16.89 against 11.90 with a seed-to-seed spread of 0.42 -- the variance component had a different posterior, not just a different origin. A chain also spends its adaptation and treedepth on a flat direction, which is what makes the field's convergence diagnostics meaningless while it does.

ModelData::set_spatial_adjacency(n_units, row_ptr, col_idx, n_neighbors) now sets the CSR arrays, the partition and the component count together; the engine's own loader uses it at both of its sites, and compute_param_layout() rejects a field whose partition does not describe its adjacency, alongside the existing PC-range-anchor guard. After the fix the same three seeds report a level of -0.003 and a tau_spatial of 12.05, inside the legacy spread.

The cause is neither a cycle nor a conditioning floor but plain local divergence. The Newton weight is the working (expected / quasi-likelihood) one wherever it differs from the observed curvature, so the error map near the mode is I - Hw^-1 Ho, which contracts only while every eigenvalue of Hw^-1 Ho is below 2. On the neg_binomial_1 fixture at phi = 6, sigma_re = 5 the largest is 2.23: the iterate walks away from the mode geometrically, losing a few parts in 1e9 of objective per step, and the stall test then reads the growing step as a floor and stops. inverse_gaussian (Fisher weight 1 / (phi mu)) does the same, so this is not one family's quirk; families whose Newton weight already is the observed curvature cannot reach it at all.

No acceptance rule reading the objective can fix that. The penalized log-posterior is stationary at the mode, so a step gains about decrement / 2 while the objective's own accumulation noise is 8 eps |obj|; on these fixtures the two cross at a joint score near 1e-6, which is exactly where every affected solve stopped. Tightening the test makes it worse -- an Armijo sufficient- decrease condition rejects genuine progress as noise and left the worst residual at 1.4e-05, three of the reported rows unimproved and one clean phi = 0.5 solve regressed from 4.9e-10 to 5.7e-07.

So the steering moves to the Newton decrement g' H^-1 g, which is already computed every iteration, is the affine-invariant distance to the mode, and keeps full relative precision exactly where objective differences are noise. Below the near-mode gate the scale the line search opens with (newton_trust_scale) halves whenever the decrement grew -- the previous step overshot -- and relaxes back toward 1 by 1.5 whenever it fell. A solve whose decrement never grows below the gate holds the scale at 1 and takes exactly the trial sequence it always did, so every converging fit, the whole #255 fixture (81 iterations at phi = 4, 323 at phi = 6) and the rational-SPDE case the stall path exists for are unchanged.

tulpa 0.0.106

What separates the two cases is whether the step is still shrinking: a converging solve sets a new shortest step every iteration however slowly, while one at its conditioning floor bounces around it (on the rational-SPDE fixture the step wanders over 3e-6 .. 4e-5 with per-iteration ratios from 0.25 to 7). The stall test now runs on that, so the affected solves are left to finish -- 81 iterations at phi = 4, 323 at phi = 6 -- and the gradient agrees to 1e-10 across every family and phi in dev_notes/probe_phi_gradient_families.R. The rational SPDE case the path exists for is unchanged (converged, n_iter 9, field correlation 0.96), and so is every solve that was already tripping the step criterion.

tulpa 0.0.105

tulpa 0.0.104

It went unnoticed because nothing downstream distinguishes the two under the canonical link: dW/deta is identically zero for gaussian_identity, so a wrong phi cancels out of the mode motion and the curvature channel. It stops cancelling for gaussian_log, whose observed-minus-working delta became nonzero in 0.0.103. Changes fitted values for existing suffixed-normal fits.

tulpa 0.0.103

gaussian and lognormal now share one set of derivatives read at their own residual -- the lognormal is the same normal density on the log scale plus a dispersion-free Jacobian. Tweedie's log-density derivative reads the mean event count off the SAME compound Poisson-gamma enumeration that produces the density, so the two cannot be taken over different truncations of the series.

An exact gradient says the optimizer walks the right surface, not that the maximizer of that surface is the generating dispersion, so each new family additionally recovers its phi from simulated truth over several seeds (test-eb-dispersion.R). Measured mean estimates against truth: lognormal 1.44 / 1.5, neg_binomial_1 2.11 / 2.0, beta 8.32 / 8.0, inverse_gaussian 0.48 / 0.5, beta_binomial 8.08 / 8.0, t 1.21 / 1.2, truncated_neg_binomial_2 2.96 / 3.0, with the random-effect SD recovered at 0.68-0.74 against 0.7 throughout. tweedie is the exception, and not for a reason of its own: the EB path never threads phi2, so it cannot be fitted there at all (gcol33/tulpa#257).

This is what unparks beta, whose derivatives had been kept unregistered since its assembled gradient landed ~1e-4 off the objective. It also makes the exact mode Jacobian and the marginal fixed-effect precision available across the registry, where before they were declined for gamma, beta, inverse_gaussian, beta_binomial, tweedie and t.

truncated_neg_binomial_2 needed the border itself fixed. It differentiates u and the phi mode motion, both formed on the true-curvature inverse, while reading Hinv -- correct only where the two weights coincide, which is what the old registration gate quietly relied on. The border now carries dH_true/dpsi, which needs the phi-derivative of the observed-minus-working correction as well as the eta-derivative the random-effect block already had; a family in that position supplies dobs_weight and .family_dphi2() withholds the whole entry without it, rather than pairing two different inverses. The second mode derivative's right-hand side reads the same correction: the score's eta-derivative is the OBSERVED curvature, and the two differ for exactly these families.

The truncated family's three weight derivatives are one chain rule over f(a, a_e) and the shape's own derivatives, the same partials curvature_deta2_for_family carries, so dweight, dweight2 and the mixed dweight_deta come from a single derivation instead of three. The bordered Hessian matches a central difference of the exact gradient to ~1e-9 for all six registered families. neg_binomial_1 keeps the differencing stencil.

tulpa 0.0.102

src/hessian_pattern_guard.h adds the drop counter, the scatter_slot() helper that single-sources every guarded slot write (14 sites across the dense-basis, indexed-cache, joint-batch and sum-to-zero scatters), and HessianPatternGuard, which snapshots the counter and raises. Detection is a counter rather than a throw at the write because the scatter runs inside OpenMP parallel regions, where an Rcpp::stop escaping the structured block is std::terminate -- the constraint that also makes LaplaceResult::start_infeasible a flag. The guard is a pure snapshot, so a driver nested inside another measures its own window without disturbing it. Placed in run_nested_laplace_grid (covering the single-block and both joint paths), the batched joint driver, spde_run_single_fit and cpp_laplace_fit_gp.

Only a NONZERO discarded contribution counts. The scatter index caches resolve whole cross products up front -- every (beta_j, RE_g) pair of an arm -- and legitimately hold -1 for pairs no observation ever touches; those slots are written with a structural zero, which changes nothing whether it lands or not. Counting at index-resolution time instead would fire on every real fit.

The shipping kernels scatter entirely inside their registered patterns: the suite runs with the check armed and zero drops, so this lands as a regression barrier rather than a bug fix. New test-hessian-pattern-guard.R covers the counter on both write paths, the raise, the zero-versus-nonzero distinction, and a real fit.

tulpa 0.0.101

tulpa 0.0.100

Warm-starting the sampler from a cheaper fit of the same model.

New:

Under a mixture the refusal now tests the BASE for zero-truncation rather than testing for a mixture at all. A hurdle's zero branch is log(pi) and carries no dispersion, so the base family's derivatives are the mixture's once masked at y = 0; genuine zero inflation on an untruncated base keeps the refusal, since its zero branch depends on phi through P(Y = 0) as well. Verified to ~1e-8 against the two-process log-marginal on a fixture with 35% zeros, and phi_hat is consistent -- 2.921, 2.871, 2.998 as the count arm grows 750, 2269, 5057 (12 seeds each, standard error 0.144, 0.091, 0.065).

The coupled y = 0 branch differentiates D = pi + (1 - pi) P(Y = 0) a fourth time, which is the only genuinely new input: it needs P(Y = 0)'s fourth log-derivative and so the second eta-derivative of the observed curvature, which has_zi_curvature_2nd_derivative() gates on. That gate is no longer narrower than the gradient's.

The five fields are checked along every route to them: most are reachable by differentiating two or three DIFFERENT third-order fields, in different directions, and all routes agree to ~1e-11 -- which is what tests the -Hess(log D) identification rather than the arithmetic. Against a central difference of the exact gradient the assembled Hessian reaches 1.5e-9 (poisson, 53% zeros), 4.1e-9 (neg_binomial_2, 57%), 5.2e-9 (binomial, 60%), and 3.3e-11 under a hurdle -- the last unchanged to every printed digit by the generalization.

The sixteen fields that needs are sixteen CALLS on one engine, not sixteen formulas: the mixture density is q(z) B(eta, phi) + C(z) in every branch, so every mixed partial of log D follows from the multivariate cumulant-from-moment recursion over the alphabet (eta, z, phi), exact at any order and any mix of directions. The base registry covers the rows the mixture leaves additively separable and the engine covers the coupled y = 0 rows; the two are summed over disjoint rows, and the engine returns exact zeros wherever the registry applies (every y != 0 row, a hurdle's y = 0 rows, a family with no free dispersion).

Anchored on zi_loglik / zi_score_eta / zi_neg_hessian and the curvature engines, all written independently of this one: 40 field checks agree to 1.9e-10, with the two shared fourth-order names each reached from BOTH members of their pair. The assembled dm/dlog_phi matches the two-process log-marginal to 5.6e-09 and the bordered Hessian a difference of the exact gradient to 2.1e-10, on a fixture with 57% zeros. The one-process path is not a second assembly: it is this one with the six z-direction fields at zero, and it expands term for term to the expression it replaces.

The same missing derivative was also what kept the family off the mixture curvature gate, so that gate now asks for the property it needs -- the observed curvature's eta-derivative, being the registered working-weight one plus the correction -- rather than naming families. beta_binomial, t and tweedie stay out, having no observed form at all. Under zero inflation neg_binomial_1's third- and fourth-order mixture fields reach 2.2e-11 and 6.3e-11, and its closed Hessian 2.5e-11.

Still on the stencil, for a stated reason rather than by omission: estimate_phi for neg_binomial_1 in any model, zero-inflated or not, since the family has no .FAMILY_DPHI entry at all -- a separate gap from the one closed here.

The returned fit carries "tulpa_glmm" / "tulpa_gam" ahead of "tulpa_fit", so print() adds the section that model class is read for -- the random-effect covariance (with VarCorr()'s estimated / sampled / conditioned label) or the smoother table -- and plot() on a GAM draws the fitted smooths. A door's print() composes the backend's own report rather than displacing it, so a nested-Laplace GAM still shows its integrated hyperparameters, grid size and outer Pareto-k.

Doors are registry-driven (.TULPA_DOORS): a new one is an entry plus a stub, and test-doors.R asserts each signature stays exactly tulpa()'s minus its withheld arguments, so a new statistical argument on tulpa() cannot leave a door silently stale.

This is what makes a Gaussian GLMM comparable to lme4::lmer(REML = FALSE). With the residual variance held at its default the outer maximization absorbs the mismatch into the random-effect scales: on a nested (1 | g1/g2) fit the inner term came out at 0.10 against lmer's 0.31. With the dispersion free the same fit reports 0.32, a residual SD of 0.5266 against lmer's 0.5260, and fixed effects within 0.006 of an lmer standard error.

Conditioned on the reference's own dispersion and given a diffuse fixed-effect prior, the Laplace mode is the MLE of the same likelihood, so the tests assert agreement to a thousandth of a reference standard error rather than to a fraction of one. The dispersion half of each likelihood is checked separately by where the profile log-marginal peaks. A parameter-recovery test against a simulated truth cannot separate a likelihood bug from sampling noise at finite N; these can.

The exact outer gradient follows the mixture. With two linear predictors the per-observation curvature is a 2 x 2 block over the count predictor and the zero predictor, and both move with the mode, so d log|H| / d theta picks up all six of its partials rather than the single dw/deta the one-process path uses. laplace_family_zi_curvature.h supplies them in closed form, contracted against the three linear-predictor (co)variances. cpp_family_has_zi_curvature_derivative() gates it, and the gate is narrower than zero inflation itself: an untruncated y = 0 branch differentiates P(Y = 0) a third time, which needs the eta-derivative of the OBSERVED curvature, and only the families whose working weight already is that have it registered. A hurdle base is zero-truncated, so its y = 0 branch is flat in the count predictor and needs no third derivative.

Refused rather than fitted as a different model: n_quad > 1 (the adaptive Gauss-Hermite inner marginal runs through a single-predictor compiled oracle), estimate_phi = TRUE (the registered dispersion derivative is the base family's, while the mixture's zero branch depends on phi through P(Y = 0) as well, so the two describe different objectives), and mode = "re_cov_gibbs" (its conditional carries no second process).

Fixed:

The same transform is now available on the NNGP spatially-varying-coefficient (spatial_svc(parameterization = "noncentered"), per term) and multi-scale (spatial_multiscale(sampler = "noncentered"), per scale) fields, which the three blocks reach through one shared per-term applier rather than a copy each. The HSGP variants of all three are already spectrally non-centered and are unchanged.

SVC keeps the CENTERED default, unlike GP. The issue that prompted this work expected SVC to share the funnel but recorded that as untested; measured on the test-svc-nuts-frontdoor.R recovery fit (Poisson, n = 120, same data / seed / budget for both arms), the centered path is already funnel-free -- 0/700 divergent, field correlation 0.939, sd ratio 0.662, mean treedepth 6.3 -- and it fits in 163s where non-centered needs well over ten times that for the same answer. The funnel that motivated the GP flip was measured on a weakly identified field, which is where non-centered earns its cost; a well-identified response does not pay it. Reach for "noncentered" when the field is weakly identified or the centered fit reports divergences.

Correspondingly there is no new SVC amplitude test: test-svc-nuts-frontdoor.R already asserts the same sd-ratio band plus a divergence guard on well-identified data, which is the correct design for that quantity -- on weakly identified data a posterior mean shrunk toward zero is the right answer, not evidence of a funnel, so a sd-ratio contrast there measures shrinkage rather than geometry.

The multiscale transform is verified by finite differences against the analytic backward on both scales, but its end-to-end amplitude test is skipped against gcol33/tulpa#244: that block's range prior is still a Uniform behind a hard -INFINITY wall (the defect #144 fixed for GP and SVC, which multiscale escaped only by being unreachable), and it leaves 82-88% of post-warmup draws divergent whichever parameterization is used.

spatial_multiscale() reaches exact NUTS through the tulpa() front door for the first time: the prior and its parameter layout existed, but no spatial-type branch ever built the sampler inputs for it, so the path was unreachable. Both scales' fields are reported as gp_local[i] / gp_regional[i] draws alongside their own (sigma2, phi). It has no nested-Laplace kernel, so mode = "structured" still routes it to the conditional Laplace path (which rejects it); auto routes it to exact NUTS, the integrator that fits it.

An intrinsic field's constant null direction is jointly unidentified with the intercept. The spatial (ICAR / BYM2) and sampler paths identify it by augmenting the field precision during the solve, so the level sits under the field's own tau; the temporal Laplace kernels never did. They identified the level only by centring the mode afterwards and folding the removed mean into the intercept -- which fixes the point estimate and leaves that direction flat in the Hessian the marginal covariance is read from. Point estimates, field shape, log-marginals and every other coefficient were unaffected, which is why nothing else moved: on the same fit the slope SE was already correct to 2%.

The RW1 / RW2 kernels now carry the same augmentation, over the field's global constant (the between-group level contrasts of a panel field stay improper by design, matching the sampler). The pin is registered for the solver to fold in rather than stored, so an RW chain's Hessian stays tridiagonal / pentadiagonal; both storages are exact and agree to five significant figures. The augmentation itself is now one helper shared with the ICAR kernels, whose numbers are unchanged.

summary() / vcov() / confint() marginalize the field out of the joint Hessian by reconstructing the field precision in R. That reconstruction augmented the ICAR structure with a unit 11', where the kernel it is meant to mirror augments with 1_c 1_c'/J_c per component -- a factor J too much weight on the constant direction, which is exactly the direction the intercept is aliased with, and no per-component handling on a replicated field. The reconstruction now matches the kernel and follows the component split for_each_icar_component walks. Slope standard errors are unchanged to five decimals, the pin being orthogonal to covariate contrasts.

The reconstruction is now pinned against the kernel's own exported log_prior_icar rather than against a second R statement of the convention: the previous reference encoded the same unit 11', so it agreed with the helper by construction and could not have caught this.

Components are now taken from the adjacency graph itself: the connected- component partition (inst/include/tulpa/graph_components.h, one DFS computed once per fit and carried on ModelData) drives every consumer -- the Laplace gradient / Hessian / pattern / log-prior, the autodiff NUTS prior, the MCAR block, and the R marginal-SE reconstruction -- so each component's constant is pinned over its actual nodes. The rank-1 sum-to-zero fold the sparse solver uses (Woodbury step and block-Schur log-determinant) now carries a node-index list, so a component that is too large to densify folds exactly whether its nodes are contiguous or not. A connected graph and an equal-size replicate are the trivial partition and stay byte-identical. ModelData grew spatial_partition and the ABI version is 37; consumer packages that LinkingTo: tulpa must be rebuilt.

An exact analytic gradient of the joint-field Laplace log-marginal with respect to the random-effect covariances, and the outer optimizer and marginal correction rebuilt on top of it.

New:

Changed:

Notes:

Identification of intrinsic fields:

Fixed:

tulpa 0.0.97

The soft sum-to-zero constant that identifies intrinsic latent fields, moved onto the reference idiom and single-sourced (#241).

Fixed:

Internal:

tulpa 0.0.96

A settable prior on the zero-inflation coefficients, and the correctness fix that finding it turned up.

New:

Fixed:

tulpa 0.0.95

One diagnostic front door, selected by draws provenance.

New:

Deprecated:

Fixed:

tulpa 0.0.94

Empirical Bayes over random-effect covariances, and the lme4 / posterior accessor surface.

New:

Correctness:

Three symptoms were live: posterior_sample() returned "iid" and tulpa_draws_array() built a 1x1x1 array from it for every Laplace-shaped fit; laplace_diagnostics(), which exists for exactly those fits, could never reach its "no posterior draws" branch; and print() on an AGQ fit reported the random-effect standard deviation under the sigma: label, where it reads as the dispersion. The remaining collisions were latent -- masked by a branch ordered ahead of them or by a companion && !is.null(...) guard -- so no coefficient, standard error or interval changes: coef(), vcov(), confint() and summary() return exactly what they did before on every backend. The draws accessors additionally route through one .fit_draws() helper.

This reaches model packages that set class = c("<model>_fit", "tulpa_fit"): they inherit exact $ too, so any of their code that was relying on a partial match now gets NULL. That reliance was always a bug, but it will surface here. Cost is ~2 us per $ read (an S3 dispatch); coef() and summary() are ~0.5 ms, so it is not measurable at the accessor level. - offset() was silently dropped by the RE-covariance backends. tulpa() never threaded the offset into tulpa_re_cov_nested() (which had no offset argument at all), so a rate model reached under mode = "re_cov_nested" -- or via the automatic random-slope redirect off mode = "laplace" -- fitted counts instead. On a simulated rate model with exposure spanning 50x and a true intercept of -1.5, it returned +1.66. tulpa_re_cov_nested() and tulpa_eb() now take offset and thread it through the inner tulpa_laplace() solve, which always supported it. Where an offset genuinely cannot be carried it now errors instead of dropping: n_quad > 1 (the compiled per-group AGHQ oracle has no offset term) and the re_cov_gibbs backend. - The outer optimization over the random-effect covariance(s) now warns when it does not converge, in both tulpa_eb() and tulpa_re_cov_nested(). It previously returned wherever the optimizer stopped without a word -- for EB that is the estimate, and for the nested path it is the centre the integration grid is placed around. control$outer_maxit (default 500) sets the budget, and tulpa_eb() reports the code as $outer_convergence. - One-dimensional outer optimization uses Brent rather than Nelder-Mead, which R warns is unreliable there. k == 1 is the common case (a scalar (1 | g) block), so this affects tulpa_re_cov_nested() as well as EB. Brent reports success at a bracket endpoint, so a variance component pinned at the bracket now warns rather than being reported as a fitted value -- the low end is the classic empirical-Bayes collapse to sigma = 0.

tulpa 0.0.93

Zero inflation as a composition over the count families, and compiled kernels for the last three families that had none.

New:

Correctness:

tulpa 0.0.92

Audit fixes (0.0.91 review, issues #228-#239).

Correctness:

Validation & docs:

Clean-up:

tulpa 0.0.91

Audit fixes (0.0.90 review, issues #218-#227).

tulpa 0.0.90

Audit fixes (0.0.89 review, issues #207-#217).

tulpa 0.0.89

Audit fixes (0.0.88 review, issues #193-#206).

tulpa 0.0.88

tulpa 0.0.87

Bug fixes and cleanups from a second whole-repo audit (issues #176-192).

The nested-Laplace tensor-grid integrator was investigated (#179) and left unchanged: CAR_proper (tau, rho) parameter recovery confirms the marginal is already in the internal log-scale parameterization, so equal-weight grid integration is correctly calibrated and adding a user-scale Jacobian biases the scale posterior.

tulpa 0.0.86

Bug fixes from a whole-repo audit.

Warnings and cleanups.

tulpa 0.0.85

Per-block quadrature order and an optional variance-component prior on the AGHQ path.

tulpa 0.0.84

Checkpoint fix (#161).

tulpa 0.0.83

Front-door API convention cleanup (#156, fully closed) and the missing-front-door features (#158, fully closed). ABI break (TULPA_ABI_VERSION 34 -> 35) for the proper-CAR exact-NUTS eigenvalue log-determinant; downstream packages must rebuild.

tulpa 0.0.82 (2026-07-15)

The sampled spatial range gets a real prior. ABI break (TULPA_ABI_VERSION 33 -> 34); downstream packages must rebuild, and any package placing a gp() / svc() NNGP block must now supply the range anchors (see the last bullet).

tulpa 0.0.81 (2026-07-15)

Third deep-audit pass: statistical, memory-safety, and backend-consistency fixes surfaced by a fan-out code audit.

tulpa 0.0.80 (2026-07-14)

CRAN-preparation release.

tulpa 0.0.79 (2026-07-12)

tulpa 0.0.78 (2026-07-12)

tulpa 0.0.77 (2026-07-10)

tulpa 0.0.76 (2026-07-10)

tulpa 0.0.75 (2026-07-08)

tulpa 0.0.74 (2026-07-07)

tulpa 0.0.73 (2026-07-07)

tulpa 0.0.72 (2026-07-07)

tulpa 0.0.71 (2026-07-06)

tulpa 0.0.70 (2026-07-06)

tulpa 0.0.69 (2026-07-06)

tulpa 0.0.67 (2026-07-03)

tulpa 0.0.66 (2026-07-03)

tulpa 0.0.65 (2026-07-03)

tulpa 0.0.64 (2026-07-01)

tulpa 0.0.63 (2026-07-01)

tulpa 0.0.62 (2026-06-30)

tulpa 0.0.61 (2026-06-23)

tulpa 0.0.60 (2026-06-23)

tulpa 0.0.59 (2026-06-22)

tulpa 0.0.57 (2026-06-22)

tulpa 0.0.56 (2026-06-21)

tulpa 0.0.55 (2026-06-21)

tulpa 0.0.54 (2026-06-20)

tulpa 0.0.53 (2026-06-20)

tulpa 0.0.52 (2026-06-20)

tulpa 0.0.51 (2026-06-20)

tulpa 0.0.50 (2026-06-19)

tulpa 0.0.49 (2026-06-19)

tulpa 0.0.48 (2026-06-19)

tulpa 0.0.47 (2026-06-19)

tulpa 0.0.46 (2026-06-18)

tulpa 0.0.45 (2026-06-18)

tulpa 0.0.44 (2026-06-18)

tulpa 0.0.43 (2026-06-18)

tulpa 0.0.42 (2026-06-18)

tulpa 0.0.41 (2026-06-18)

tulpa 0.0.40 (2026-06-17)

tulpa 0.0.38 (2026-06-17)

tulpa 0.0.37 (2026-06-16)

tulpa 0.0.36 (2026-06-16)

tulpa 0.0.35 (2026-06-16)

tulpa 0.0.34 (2026-06-12)

tulpa 0.0.33 (2026-06-12)

tulpa 0.0.32 (2026-06-12)

tulpa 0.0.31 (2026-06-12)

tulpa 0.0.30 (2026-06-10)

tulpa 0.0.28 (2026-06-09)

tulpa 0.0.27 (2026-06-09)

tulpa 0.0.26 (2026-06-08)

tulpa 0.0.25 (2026-06-08)

tulpa 0.0.24 (2026-06-08)

tulpa 0.0.23 (2026-06-08)

tulpa 0.0.22 (2026-06-08)

tulpa 0.0.21 (2026-06-08)

tulpa 0.0.20 (2026-06-08)

tulpa 0.0.19 (2026-06-08)

tulpa 0.0.18 (2026-06-08)

tulpa 0.0.17 (2026-06-07)

tulpa 0.0.16 (2026-06-07)

tulpa 0.0.15 (2026-06-07)

tulpa 0.0.14 (2026-06-06)

tulpa 0.0.13 (2026-06-06)

tulpa 0.0.12 (2026-06-06)

tulpa 0.0.11 (2026-06-06)

tulpa 0.0.10 (2026-06-05)

tulpa 0.0.9 (2026-06-03)

tulpa 0.0.8 (2026-06-02)

tulpa 0.0.7 (2026-06-01)

tulpa 0.0.6

tulpa 0.0.5

tulpa 0.0.2

tulpa 0.0.3 (2026-05-28)

tulpa 2026-05-13 — ABI v13: Phase D — delete legacy ratio path

Closes the tulpaRatio migration tracker (gcol33/tulpa#15). After v12 gated the legacy ratio body of compute_log_post_impl<T> behind a generic-layout check, Phase D removes the body itself and every consumer of LegacyRatioData / LegacyRatioLayout.

Downstream rebuild notes: tulpaRatio already routes through the generic LikelihoodSpec path via tulpa_bridge.cpp + per-family payloads in lik_specs/ (B1+B2 of the migration); rebuild against v13 headers, no logic changes needed. tulpaObs never used the legacy ratio path; rebuild against v13. * tulpaGlmm Day-22+ already targets the generic path; rebuild against v13.

tulpa 2026-05-12 — ABI v12: generic-layout safety in compute_log_post_impl + ESS port

tulpa 2026-05-12 — ABI v11: caller-supplied inv-mass diagonal for NUTS

tulpa 2026-05-11 — NNGP Laplace: full off-diagonal precision scatter

tulpa 2026-05-11 — nested-Laplace ST family: 5 more indexed × indexed combos

tulpa 2026-05-11 — nested-Laplace joint spatial × temporal (ICAR × AR1)

tulpa 2026-05-11 — nested-Laplace HSGP returns modes + store_Q

tulpa 2026-05-11 — nested-Laplace BYM2 returns modes + store_Q

tulpa 2026-05-11 — nested-Laplace store_Q on RW1/RW2/AR1/CAR_proper

tulpa 2026-05-06 — Takahashi partial inverse as a registered C-callable

tulpa 2026-05-05 — multi-term + slope REs on the spec-Laplace path



Try the tulpa package in your browser

Any scripts or data that you put into this service are public.

tulpa documentation built on Oct. 19, 2026, 5:09 p.m.