latent_factor() gives each process its own loading on each factor (#134).
The loading matrix is lower triangular with a positive diagonal, so one
factor adds a rank-one and two factors a full residual covariance between
the numerator and denominator predictors; more factors than processes are
refused, since the extra ones are not identified. Before, every factor
loaded with one sigma on both processes, so K factors were one effect and
n_factors = 2 gave about 80 divergences. The factor scores carry their
N(0, 1) prior on the sampled scale, which keeps the direction the
sum-to-zero constraint removes proper. Draws are named
lambda_latent[process,factor]; the scale argument, which nothing read,
is gone.
A per-group temporal walk (temporal_rw1(group_var = ),
temporal_rw2(group_var = )) sums to zero within each group (#135), so it
carries the group's trend and a group intercept such as (1 | site) its
level. Before, only the global mean was removed and the walk levels lay
along the intercept's direction: 25 divergences at 1000 iterations with
(1 | site), 95 in the non-centred coordinate; now 0 to 9 and 0.
spatial_multiscale() represents both scales on one Hilbert-space basis
(#136). Two nearest-neighbour GP fields over the same locations leave the
split of their variance to the priors, and NUTS diverged on most draws there
in both coordinates (100 to 190 of 200 on the 40-location fixture, also
with the variances sampled as a total and a share); the basis samples the
same fixture without divergences. approx, nn_local, nn_regional,
parameterization, cov and nu are gone (the basis is squared
exponential, and cov was ignored on it), and sampler = "lbfgs" now
reaches the field.
The PC prior on an HSGP field's sd is proper (#146). On the log-variance
coordinate it lacked the d sigma / d log sigma^2 Jacobian, so its density
tended to a constant as sigma went to zero (single-scale and multi-scale
HSGP).
A multi-scale field combined with anything beyond a temporal walk (a space-time interaction, latent factors, a varying coefficient) receives its gradient (#147); the composite gradient left all of its parameters at zero.
spatiotemporal(parameterization = "collapsed") integrates a GP space-time
interaction (type = "separable" or "nonsep_gp") out with a Laplace
approximation inside HMC (#143; Margossian et al. 2020): the sampler moves
the interaction's three hyperparameters and every other block, its gradient
follows from the implicit function theorem with one adjoint solve, and the
interaction is drawn from its Gaussian approximation at each saved
iteration. The soft margin constraints enter as a low-rank term through a
quasi-definite augmented system, so one sparse factorization serves every
solve. On the spatial-temporal vignette's model (25 regions x 4 years,
Gneiting, nn = 15) NUTS takes about 20 leapfrog steps per draw instead of
about 390, and 400 iterations of one chain take about 2 minutes where 200
centred ones took 10. For a likelihood Gaussian in the interaction the
approximation is exact, and a test checks the density against the centred
one there to 1e-8. It needs a log-concave family with canonical links and
no inflation.
The composite and latent-factor gradients read an observation's residuals and dispersion derivatives from the log posterior's own per-observation density (#148). Each carried its own family algebra: on gamma-gamma and lognormal both were off by order one, beta-binomial dispersions got no gradient, and neither carried zero or one inflation. A collapsed GP or areal field refuses the families and inflation its inner Laplace is not written for (gamma-gamma, lognormal, beta-binomial; negative binomial-gamma for the GP). The gradient sweep now runs every hand-coded family and inflation structure through these paths.
An intrinsic CAR or BYM2 field on a disconnected map (#141) is centred on each connected component rather than over the whole graph, an island gets a unit-variance prior instead of none, and BYM2 carries tulpa's per-component scaling (Freni-Sterrantino, Ventrucci and Rue 2018). HMC, ESS, SGHMC, SGLD, VI and the Laplace backend fit it; the Gibbs and Polya-Gamma samplers, which centre over the whole graph, and the collapsed parameterization refuse it by name, and auto mode sends it to HMC.
ESS, SGHMC, SGLD and VI fit BYM2 again (#145). Their parameter names are now written at the indices the layout assigns; before, a BYM2 field stopped with "named 23 of 22 parameters". ESS also updates the unstructured BYM2 component, which it had left at its starting value.
coef() and vcov() of a single-process Laplace fit (binomial,
beta-binomial) return the fixed effects (#144). Its draws were named by
design column rather than beta_num[j], so both came back empty.
A Gibbs fit whose numerator has no dispersion parameter (e.g.,
ratiod_poisson_gamma()) returns its draws; it stopped on the missing
numerator dispersion name, which broke the BYM2 chunk of the
spatial-temporal vignette.
The BYM2 field has one scale on every backend (#140): the structured
component is multiplied by tulpa's Riebler et al. (2016) scale, 1 / sqrt(s) with s the geometric mean of the ICAR marginal variances. The
HMC-family backends computed their own factor from eigenvalues and the Gibbs,
Polya-Gamma and Laplace backends multiplied by the eigenvalues' geometric
mean itself, so sigma and rho were different parameters on each: on a
5 x 5 grid HMC reported sigma 8.1 and Gibbs 0.89 for the same field.
The Gibbs sampler fits BYM2 again (#133). The fixed effects and both field components move as one Gaussian-approximation block, and (sigma, rho) jointly with them (Knorr-Held and Rue 2002); on the 25-site test field Gibbs and HMC now agree on every hyperparameter to within 0.1 posterior SD. Auto mode sends BYM2 models to Gibbs, as it does ICAR ones.
spatiotemporal(parameterization = "noncentered") reaches the GP
interaction types (#132): z ~ N(0, I) through the nearest-neighbour GP
autoregression, delta_i = m_i + sqrt(d_i) z_i, sampled with the arena
gradient. Stored draws are the interaction itself. The centred coordinate
stays the default: on an 8 x 4 grid it mixed as well or better and ran
faster, for an interaction SD of 0.5 and of 0.05 alike.
Families take the links their likelihoods now implement (#117): "probit"
and "cloglog" beside "logit" for a success, zero-inflation,
one-inflation or hurdle probability, and "identity" and "inverse" beside
"log" for a gamma mean. One inverse-link helper carries every link through
the log posterior, and a model with a non-canonical link takes the arena
gradient of it. fitted(), predict(), ratio(), loo()/waic(),
posterior_predict() and prior_predict() apply the family's links. The
Gibbs, Polya-Gamma and Laplace samplers and the collapsed spatial fields keep
the canonical links and name the restriction.
A known exposure (#118): ratiod_zipois(), ratiod_zinegbin() and the
hurdle families take denom_family = "fixed", and ratiod_lognormal()
denom_fixed = TRUE. The exposure scales the numerator mean, has no
coefficients, and the ratio is the rate. ratiod_zipois() and
ratiod_hurdle_pois() also take denom_family = "negbin", a Poisson
numerator over a negative binomial denominator.
ESS, SGHMC and SGLD no longer fit every family outside binomial, negbin-negbin and Poisson-gamma with the negbin-negbin likelihood, and VI no longer refuses them; a gamma-gamma or log-normal model with a GP field no longer reads its continuous numerator out of bounds (#139). Every sampler entry point reads the model type, links and responses through one builder and names its parameters off the shared layout.
latent_factors() and spatiotemporal_effects() are tulpa's verbs,
re-exported, instead of tulpaRatio copies that masked tulpa's on attach
(#100). The spatiotemporal summary and long formats now label region s and
time t correctly when the numbers of regions and times differ: the old
copy read the s-fastest draw layout t-fastest. format = "summary" returns
a tulpa_st_summary (with print() and a plot() that adds
type = "spatial_map"), and the array and long formats are a plain array
and data frame. Requires tulpa >= 0.7.5.
spatial_gp() fits no longer run away. Outside the range prior's bounds
the hand-coded GP and multi-scale GP gradients returned a zero gradient
without reporting the log-posterior, so NUTS read an out-of-support point as
carrying the previous point's density; once the range crossed its upper bound
the chain random-walked to beta near 1e6 with phi_gp = Inf (#102).
NNGP fits are several times faster: the centred gradient computes the field's
density in the same pass and the fused log-posterior reuses it, the
non-centred density no longer rebuilds a field it does not read, and every
neighbour covariance is filled once per pair on the cached distances. A
60-location spatial_gp() fit drops from 69 to 13 s. The SVC terms use the
spatial GP's NNGP density and gradient instead of their own copies, and the
block-diagonal warmup metric shrinks its covariance as the diagonal one does
(#101).
predict() builds new-data design matrices from the training terms, factor
levels and contrasts, so a factor predictor, or new data holding only some of
its levels, predicts (#103). predict(coords.0 = ) krigs spatial_gp() fits
from the stored field draws and evaluates spatial_hsgp() fits on the fit's
own sine eigenbasis, both on the training coordinate scaling (#104).
loo(), waic(), ratiod_compare() and ratiod_average() evaluate each
observation's likelihood at the full linear predictor, random and
structured effects included; non-HMC draw backends refuse models with such
terms (#108). loo() and waic() are re-exported.
summary() returns a summary.ratiod_fit object whose print() method lays
it out; coef(summary(fit)) is its parameter table. coef(), vcov(),
nobs() and formula() methods are new for fits, and as.data.frame()
returns the posterior draws again (#105).
The HMC sampler records each iteration's Hamiltonian, so plot_energy() and
diagnostic_summary() report E-BFMI, per chain through
tulpa::energy_bfmi() (#107). diagnostic_summary() returns its object
instead of printing it, without terminal colour codes.
latent_factors() reads a tratio() fit's factor scores (#106), and
validate_mode() accepts the backend a fit ran on ("hmc") as well as its
tier.
A priors = ratiod_priors(...) argument reaches the samplers: every MCMC
backend translates it through one function, passing Normal(0, sd) on fixed
effects, Half-Cauchy on random-effect SDs and Gamma on the dispersion, and
refusing other families by name instead of fitting under its own defaults
as before (#109). A component left unspecified takes the default the
samplers already used (priors_default() lists them per family), so default
fits are unchanged; prior_predict() and sim_ratiod_sbc() now draw from
those same defaults. ratiod_priors_legacy() is removed.
summary() reports every main parameter; it stopped at the twelfth.
print() on a fit is a five-line description.
ratio_contrast() compares the levels of a factor, character or logical
variable in level order and refuses a numeric covariate.
(1 | a:b) groups by the interaction of a and b, as the second term of
(1 | a/b) does.
tratio() is quiet by default: control = list(verbose = TRUE) prints the
progress messages and sampler adaptation traces.
Help pages: one Title Case style, print/summary/plot methods documented on the page of the function returning their class, and every exported topic with an example that runs.
The HMC entry points run their OpenMP regions on the threads cores
allows, so concurrent chains no longer each start a full team; a
zero-inflated test fit drops from 87 to 1 s (#110).
Correlated random slopes take their correlation through canonical partial correlations with an LKJ(2) prior and its Jacobian, and draws unpack with three or more correlated coefficients. Correlated-slope posteriors change (#111).
Inline transforms such as log(x) and I(x^2) combine with random-effect
terms: the random part is removed from the formula as a parse tree, not as
text (#112).
loo(), waic() and posterior_predict() support every zero-inflated,
one-inflated, zero-and-one-inflated and hurdle family; the zero-and-one
inflated binomial is the extended mixture its density in the sampler uses
(#113).
Hurdle fits no longer gain mass where a small mean makes P(Y = 0) round to
one: the zero-truncation uses a stable log(1 - exp(x)) instead of a
-1e10 sentinel the truncation subtracted (#114).
mode = "ess" targets the posterior (the block's reference density was
counted twice and the slice threshold was stale), runs the requested
chains, and returns natural-scale draws that fitted(), predict() and
loo() read as they read an HMC fit's. Its reference starts at a Laplace
approximation and is refitted from warmup draws as a multivariate t, as in
generalised elliptical slice sampling; on a random-intercept model its bulk
ESS now exceeds HMC's at the same length (#115). ESS refuses random slopes
and several grouping factors instead of dropping them.
Families accept only the links and denominators the samplers fit; other
links, denominator_known = FALSE, denom_fixed = TRUE and several
denom_family values used to be stored and silently fitted as the default
model. denominator_known and denom_fixed are removed;
ratiod_zinegbin() and ratiod_hurdle_negbin() gain
denom_family = "gamma" (#116).
Fixed-effect draws are named num:<term> and denom:<term> in coef(),
vcov() and the printed summary() tables (#127); the unused
latent_mean[k] draws are gone.
The Gibbs sampler fits ICAR models only. Its BYM2 sweep shifted the
structured field's level against the intercept, off the field's
sum-to-zero constraint, so the intercept wandered under its own prior
(#119); without that move it does not reach the posterior HMC samples, so
a BYM2 model in mode = "gibbs" is refused and mode = "auto" fits it
with HMC (#133). Gibbs fits report natural-scale draws under the HMC names
and work with ratio(), fitted() and loo() (#120).
The PC prior on the space-time interaction's precision was the density of a
variance, which made it improper toward large precision; with little
interaction in the data tau_st ran to 1e9 (#130). Space-time results
change. Non-centred Type IV fits rebuild the interaction from its scale for
fitted values, ratios and loo() (#131).
Collapsed ICAR, BYM2 and GP fields reach the fit for every chain count and
entry point; multi-chain collapsed fits lost them, and their fitted values,
ratios and loo() missed the spatial term (#125).
Correlated random slopes report their correlations, cor_re[t,a:b],
instead of the sampler's raw coordinates (#122).
The pre-sampling gradient check compares against a five-point finite difference. The two-point difference was too coarse for a covariate on a large scale, so the check fired and the fit switched to slow finite-difference gradients (#124).
Divergent transitions raise a warning of class tulpaRatio_divergence
whether or not verbose is set (#123), and print(summary(fit)) flags
parameters with R-hat above 1.01 or bulk ESS below 100 per chain (#126).
spatiotemporal(parameterization = "noncentered") covers Knorr-Held Types
I to III as well as IV: the interaction is sampled as z with
delta = z / sqrt(tau), and its sum-to-zero penalty acts on z, so a weak
interaction no longer leaves tau in a funnel (#132).
The centred space-time interaction's normalizer counts only the directions
its sum-to-zero constraints leave; Types I to III and the GP types counted
the constrained directions as well and overstated tau (Type I by about two
posterior SDs on a 25 x 6 grid). Space-time results change (#138).
A NUTS leaf whose step leaves the posterior's support counts as a
divergence with acceptance zero; std::min(1, exp(-NaN)) had reported it as
perfect acceptance, so the step size adaptation could grow epsilon to its
cap and freeze a chain, most often under longer warmups (#137).
ratio(newdata = ) works (#121); print() lists every grouping factor
(#128); trace plots drop overlapping axis labels (#129).
spatial_svc() samples each NNGP term non-centred by default (#96). The
new parameterization argument takes "noncentered" (default) or
"centered". Non-centred, each term's block holds z ~ N(0, I) and the
field is w = L(sigma2, phi) z through the NNGP autoregression the spatial
GP already uses, so the variance reaches the field through the transform.
Centred, the field's level is identified by the NNGP prior alone and wanders
with sigma2: on the 80-row recovery fixture, five chain seeds at 4 x 3000
gave intercept R-hat up to 1.117 and a sigma2 bulk ESS of 12, against at
most 1.020 and at least 167 non-centred. Both coordinates place the same
prior on w; a spec built before the argument existed runs centred.
An NNGP SVC range no longer flags most transitions divergent (#99). The
range carried its uniform prior as a hard bound on log(phi), so a proposal
past either end of range met a log density of -Inf, which the sampler
records as a divergence; with the posterior near the upper end, 442 of 500
transitions were flagged on an 80-row fixture. The range is now sampled on
the logit scale of its interval, which places the same uniform prior with no
edge in the sampled space: 2 of 500 centred and 0 of 500 non-centred on the
same fixture and seed. The draws still report phi_svc itself.
A stored SVC draw is the field the likelihood reads (#98). eta reads each
term centred over locations, but the sampler stored w uncentred, so
svc() and the fit's eta carried a level the likelihood removes: 0.73 on
average on an 80-row fixture with a true slope of 0.3. The draw store now
writes the centred field, in either coordinate.
spatial_rsr() is carried, or refused, by every backend (#92). HMC
received the restriction only on its centred NNGP branch, where it projected
the field and then evaluated the NNGP prior on the projected values, with an
N x N projection applied to a field held at the unique locations; on an areal
field the entry point cleared the flag, so spatial_rsr(spatial_car(...))
gave draws bit-identical to the unrestricted fit. The log posterior now
builds the spatial field's value at each observation once (areal, NNGP,
HSGP or multi-scale GP), and under RSR the predictors read (I - Q Q') f
with Q the orthonormal basis of the restricted covariates; every field
prior reads f. The gradient of a restricted model is the arena gradient of
that density, and fitted() applies the same projection. A collapsed
parameterization is refused under RSR. The PG sampler drew the field from
the unrestricted ICAR conditional and read a corner of the projection; it now
draws it from its exact Gaussian conditional, with precision
tau Q + A' Omega A for A = (I - Q Q') S, and its R wrapper no longer
errors before sampling. ESS, Gibbs, SGHMC, SGLD and VI, which never read the
restriction, refuse it, and PG refuses it on anything but an intrinsic CAR
field. predict() refuses a restricted field at new rows.
predict() on the fitted rows reproduces fitted() (#93, #94, #95).
Three separate defects pulled the two apart:
log(p) and took plogis() of that, which is not the predictor; the field
now enters the linear predictors and the response is taken once through
hmc_response_draws(), as fitted() takes it (#94).s paired
the parameters of iteration s with the predictor of iteration s - 1; the
predictor is refreshed at the save (#95). The field-free and
random-effect-only fits run tulpa's own sampler, which stores the saved
draw's predictor from gcol33/tulpa#941 on.The log posterior no longer depends on the thread count on Windows arm64
(#83). On aarch64-w64-mingw (Rtools45 clang 19) an OpenMP
reduction(+:double) over a 4-thread team returned the serial value in 0 of
500 repeats of a plain Poisson sum on a 4-core runner, while 1 and 2 threads
and a hand-rolled per-thread sum were exact. The observation sum in the
density, the Laplace log-likelihood and the test kernels now go through
ratiod_omp::sum_range() (src/omp_sum.h): each thread sums one chunk into
its own slot and the slots are added in thread order. The HMC gradient's
observation loop summed the zero- and one-inflation coefficients, the
temporal and spatial fields and the random slopes with omp atomic, and on
the same platform the analytic beta_zi gradient missed finite differences
by up to a factor of four; those sums now share the loop's per-thread slots
with the other blocks.
The Laplace tier runs on the nested-Laplace engine (#85, #84). The
in-tree mode-finder (laplace_core.cpp) picked one value per hyperparameter
and conditioned the fixed effects on it, so intervals omitted hyperparameter
uncertainty. The tier now hands the model to tulpa::tulpa_nested_laplace_joint():
one likelihood arm per process, one latent block per random-effect term,
areal field (ICAR, BYM2, proper CAR), Gaussian-process field (single and two
scale, on a Hilbert-space basis) and multi-scale temporal component. The
engine integrates every block hyperparameter and every dispersion on its outer
grid. summary() reads the fixed effects off the grid mixture (mean,
covariance, skew-aware interval), so the rows no longer move with a draw seed;
hyperparameters are reported as posteriors instead of constant columns; the
stored draws are drawn from the mixture.
Consequences of the move:
poisson_gamma, negbin_gamma, negbin_negbin,
gamma_gamma) are fitted as two arms with their own fixed effects and
dispersions, shared latent blocks entering both. The old tier fitted the
numerator alone, so ratio() of such a fit was a numerator rate (#89).predict() works for fits without a spatial or temporal field (it always
errored, #87), and fitted() reads each random effect at the group each row
belongs to (#90).spatial_hsgp() are fitted.(I - Q Q') S z) or as a projected
Hilbert-space basis; the projected effect is orthogonal to the restricted
covariates. Spatially varying coefficients have no block on the engine and
are refused with mode = "hmc" as the way out.summary() lists every fixed effect of a one-process fit; it listed the
intercept only (#91).The SoftAbs divergence retry is replaced by a WALNUTS transition (#86).
control$riemannian used to resolve to an automatic setting that switched
the retry on for BYM2 and ICAR fits with a dense mass matrix. The retry
replaced a divergent NUTS trajectory with a second one under a
Hessian-based metric, so which kernel moved the chain depended on the
current state and no acceptance correction covered the choice. On the
engine's copy of the scheme (gcol33/tulpa#937, Neal's funnel) a seed settled
in the neck with 451 divergent transitions against 6 without it. The retry,
its finite-difference Hessian and control$riemannian are removed.
control = list(walnuts = TRUE) runs each trajectory as the engine's
WALNUTS transition (<tulpa/walnuts.h>, the code tulpa's own sampler runs):
every macro step takes the coarsest leapfrog subdivision that keeps the
joint log density within tolerance, and only if the reversed step picks the
same one, so the kernel is exact. On a BYM2 field over 12 units with one
observation each at adapt_delta = 0.5, NUTS diverged 2 to 122 times over
5 seeds and WALNUTS never, and both sat within 0.06 reference SDs of a long
adapt_delta = 0.95 NUTS run on every fixed effect and hyperparameter, with
SD ratios 0.93 to 1.05 under WALNUTS. NUTS stays the default, and its draws
are bit-identical to those before the change. Needs tulpa 0.7.1 or later,
the first version that ships <tulpa/walnuts.h> (gcol33/tulpa@7fddf507).
mcmc_diagnostics() is gone; use diagnostics() (#82). It was an
exported lifecycle::deprecate_warn() shim forwarding to diagnostics(),
and the engine's shim of the same name was deleted in gcol33/tulpa#773, so
two packages exported the name and whichever attached second masked the
other. The package is pre-release, where the convention is to delete the old
name rather than carry a shim. lifecycle was reachable from nothing else
and leaves Imports with it.
The standalone intrinsic walks take the same coordinate (#79).
temporal_rw1() and temporal_rw2() take parameterization, reading the
same transform the multi-scale arms do. temporal_effects() and
temporal_gmrf_prior_grad() are the two places it lives -- the one chokepoint
every eta assembly reads and the one all eight gradient functions reach -- so
the coordinate is a branch in each rather than a change to any caller.
Measured on a 500-observation fixture with no temporal signal, centred ->
non-centred: an RW1's bulk ESS for tau_temporal goes 409 to 805 at 0 and 1
divergences, and an RW2 picks up 23 divergences of 2000 for no ESS gain. Both
return the same posterior, the intercept and slope agreeing to four decimals.
Non-centred helps where the variance is shrunk toward zero and can cost where
the data pin the field, which is why it is an argument rather than a default.
The non-centred walks take group_var (#81). The prior augments ONE
global constant whatever the number of groups, so G - 1 group contrasts stay
improper: a per-group level shifts eta only for that group's observations and
the likelihood identifies it. Running the single-group transform per group
would augment each group's own constant and shrink those contrasts, which is
a different model, so the transform re-mixes instead. Each group's coordinate
splits into its increments and one sum coordinate; the G sum coordinates
split in turn into the global direction the augmentation fills and G - 1
contrasts, which stay unscaled because they are flat and scaling them would
couple the group levels to tau. The cyclic branch reads its mean coordinate
the same way, so one map covers RW1, RW2 and the seasonal walk.
cpp_ms_temporal_nc_coordinate_spread() now runs over groups, and
cpp_ms_temporal_nc_group_contrast_flat() moves the contrast directions
alone: it asks that the prior stay where it is and that group g's field move
by exactly the unscaled offset. A finite difference cannot separate a
transform that augments per group -- that density is self-consistent and its
gradient matches its own difference quotient -- and neither can a flatness
check on its own, which a transform that dropped the contrasts would pass.
Every gradient fixture is checked to build the structure its name claims.
make_model() selects a structure by a chain of flags, each naming the fields
it covers, so a name added to KNOWN_FIELDS but to no flag falls through every
branch and yields a model with no structured field: the gradient check then
passes on a two-parameter regression while reporting the structure's name.
rw1_nc did exactly that, coming back with 2 parameters against the 13 its
siblings carry. test-gradient-correctness.R now asserts every field builds
more than the structureless model.
A non-centred coordinate for the multi-scale temporal arms (#79).
temporal_multiscale() takes parameterization, and "noncentered" samples
the trend and seasonal arms as z ~ N(0, I), reaching the effects through
sigma * A^-1 z. A stacks the arm's difference operator on 1'/sqrt(n),
which is exactly the augmented precision's factor Q + 11'/n = A'A, so the
arm's quadratic form and the transform's Jacobian cancel term for term and
what is left is N(0, I) with the variance reaching the linear predictor
through the transform alone. A^-1 is a cumulative sum, once for RW1 and
twice for RW2, so the transform and its adjoint are O(n).
A non-cyclic RW2 annihilates a linear ramp as well as the constant and only
the constant is augmented, so one coordinate stays a free direction under the
flat prior it already had. Any right inverse of A reproduces the density,
but not any right inverse is well conditioned: the double integral's own
particular solution carries a large ramp that the free coordinate can undo,
which measured 406 divergences of 2000 before the ramp was projected out and
3 after.
Measured on the multiscale RW2 recovery fixture, whose true trend is zero so
the whole posterior sits in the funnel's neck, centred -> non-centred: the
divergence count goes 30 of 2000 to 3, and by decile of log sigma2_trend the
neck goes from 12.5% divergent to 0%. The bulk ESS of sigma2_trend goes 21
to 310 and its rhat 1.126 to 1.013, and the intercept's ESS 596 to 1776. The
centred run reports a posterior median of 1.6e-04 against 4.2e-05, biased
upward by the part of the neck it never reaches; only the non-centred run has
converged on that parameter. The intercept and slope are unmoved.
test-gradient-correctness.R carries the three walks in either coordinate
across all four gradient modes, and
cpp_ms_temporal_nc_coordinate_spread() checks the invariant a finite
difference cannot see: a transform that builds the wrong field leaves the
density self-consistent and its gradient still matching its own difference
quotient, but at fixed hyperparameters the two coordinates' densities have to
differ by one number for every z. Measured spread 1e-13.
A proper CAR field's mean is removed from the linear predictor (#80).
Q(rho) = D - rho*W is full rank at rho < 1, so the prior identifies the
field's mean, and the field was left uncentred on that reasoning. The
likelihood cannot see that mean apart from the intercept though: eta reads
beta_0 + phi_s, and moving the intercept by c against every phi by -c
leaves eta where it was. Measured on a 36-unit fixture, the two correlated
at -0.977, each carried a posterior sd near 0.06 while their sum carried
0.013, and the intercept's bulk ESS was 71 of 2000 at rhat 1.053. Longer
warmup made it worse, not better: 1.053 -> 1.112 and ESS 71 -> 32 going from
1500 to 3500, because adaptation tunes for the bulk and the ridge is not in
it.
spatial_car() now takes center, TRUE by default, and an SVC term's
treatment is what it applies: eta reads the centred field and the prior
still reads the raw one. Since Q(rho) does not annihilate the constant, the
two are different quadratic forms, and the density and every gradient now
take the prior's view of the field separately from eta's.
The identified posterior does not move. On that fixture the level comes back at 0.5200 against 0.5204, the slope at 0.3003 against 0.3004, with rho and tau inside their own uncertainty. What changes is that the level is the intercept: ESS 71 -> 1953 and rhat 1.053 -> 1.001.
center = FALSE keeps the previous parameterization and requires
proper = TRUE; an intrinsic field's mean is unidentified until it is
removed, so there is no model on the other side of that switch.
car_proper_raw covers it in test-gradient-correctness.R.
A proper CAR field reaches only gradients written for it. GF_AREAL
covered a proper CAR field as much as an intrinsic one, so
spatial_car(proper = TRUE) alongside temporal_multiscale() or a
spatiotemporal interaction selected
compute_gradient_ms_temporal_handcoded or
compute_gradient_spatiotemporal_handcoded. Neither writes proper CAR's
rho gradient or its log-determinant, and both centred a field the log
posterior leaves alone, so the gradient they returned was the gradient of a
different density. Measured against central differences of the log posterior
on a 400-observation, 25-unit fixture: 210 on the regression coefficients,
where every field in the sweep agrees to better than 1e-04 relative.
Proper CAR now carries GF_AREAL_PROPER, so a function whose covers() mask
does not name it falls through to the composite -- the path
spatial_car(proper = TRUE) on its own already took, and the one that
carries the rho block. The same combinations now agree with finite
differences to 2e-07 and 3e-06. car_proper_ms and car_proper_st join
test-gradient-correctness.R, which is what nothing reached before.
One predicate decides whether the spatial block is centred. The rule was
written out at five sites, and three spellings had drifted from the log
posterior's. spatial_block_is_centred() (hmc_sampler.h) is now the only
one, so a gradient and the density it is evaluated against cannot disagree
about which field the likelihood saw.
Recovery assertions read the quantity their model identifies. The SVC recovery test compared a single 80-observation realization's posterior mean against the value the data were simulated from, where the maximum likelihood slope carries a standard error near 0.05: seed 11's sample implies 0.2533 and seed 22's implies 0.2062, so the assertion failed on sampling noise. It now runs the three seeds the RW1 test runs, asserts the posterior mean against the same sample's unpenalized MLE -- measured agreement within 0.002, an order of magnitude sharper than a tolerance against 0.3 can be -- and keeps recovery of the simulated values as interval coverage across seeds.
The proper CAR test asserted rhat < 1.05 on the intercept, which that model
does not identify on its own: the intercept and the field's mean correlate at
-0.977 and their sum carries a fifth of either one's posterior sd. It now
asserts on the sum, which converges (ESS 1970 of 2000, rhat 1.002), and on the
intercept's interval coverage.
The multiscale RW2 fixture's divergences are pinned where they belong.
They sit in the neck of sigma2_trend, whose true value on that fixture is
zero: 12.5% of draws in the lowest decile of log sigma2_trend diverge and
none above the third, and the divergent draws sit 1.96 posterior sd below the
clean ones. test-recovery-constrained-fields.R asserts that location rather
than a bare count, which is what separates the known centred-parameterisation
funnel from a fault elsewhere. The trend arm's non-centred coordinate is #79.
The SVC range prior a spec declares is the one the fit runs (#78).
prepare_svc_for_hmc() wrote phi_prior_lower = 0.3, phi_prior_upper = 30
as literals into each of its three branches, so every SVC fit ran on that
interval and nothing a caller set could change it. spatial_svc() carried no
field for the bounds at all, and the %||% 0.01 / %||% 10.0 fallbacks in
fit_hmc() -- the numbers ?spatial_svc and the ModelData defaults both
named -- were unreachable, because svc_info always supplied a value for the
%||% to pass through.
spatial_svc() now takes range = c(lower, upper), validated as
0 < lower < upper the way spatial_multiscale() validates range_local
and range_regional, and prepare_svc_for_hmc() reads it off the spec.
svc_default_range() is the one place the default lives; the ModelData
defaults carry the same pair, and the print method reports the interval in
force. Fits that set nothing are unchanged: the default is the [0.3, 30]
those fits already ran under, and the draws are bit-identical.
The bounds apply to approx = "nngp". The HSGP arm parameterises a
lengthscale under LogNormal(0, 1), which this uniform does not set, so a
range passed alongside it is an error rather than a value nothing reads.
ess_sampler.cpp and vi_sampler.cpp each wrote a third and fourth pair of
bounds into a block they declare unsupported (has_svc = false); both are
gone.
Intrinsic fields identified by augment-and-centre (#22). The temporal
random walks and the multi-scale trend and seasonal arms were identified by a
soft penalty on the field's sum, -0.5 * lambda * (sum phi)^2 at
tulpa::s2z_precision(n). That pins the constant direction without removing
it: the direction stays in the linear predictor at precision `1 / (kappa*n)^2
, five orders of magnitude above the tau` the field's own directions
carry, and a step size that holds it stable is the one the whole field then
has to move at.These fields now use the engine's construction (tulpa/sum_to_zero.h,
gcol33/tulpa#241): the precision is augmented, Q -> Q + 11'/n, so the
constant carries the field's own tau, and the effects are centred on their
way into eta, so the likelihood never sees it. The normalizer takes the
augmented rank -- one more pinned direction, except that a non-cyclic RW2
keeps its linear null direction and stays deficient by one. Measured on a
500-observation binomial fixture, the adapted step size rises from 0.101 to
0.351 for temporal_rw1() and from 0.105 to 0.255 for temporal_rw2(), and
the leapfrog steps per iteration fall from 35 to 12 and from 42 to 20; the
multi-scale arms move from 0.066 to 0.124 and from 0.040 to 0.056.
The reported field changes with it. fitted(), ratio() and predict() add
the stored field back into eta, so an intrinsic field's draws are now stored
centred -- the series the likelihood saw, rather than one differing from it by
a level. temporal_ar1(), temporal_iid() and the temporal GP are proper,
identify their own level, and are unchanged.
spatial_multiscale(parameterization = "noncentered") (#75). The
multi-scale field was the one continuous field in the package with a single
coordinate: spatial_gp() offers "centered", "noncentered" and
"collapsed", and temporal_gp() and temporal_ar1() each offer a second
coordinate, while the multi-scale block carried two variance parameters over
the same locations with no alternative to the effects themselves. #36 removed
the name that selected nothing; it now selects a coordinate.
Non-centred, each scale's block holds z ~ N(0, I) and its field is
w = L(sigma2, phi) z through the NNGP autoregression the single-scale field
already uses, so all four hyperparameters reach the linear predictor through
the transforms as well as through the priors. Stored draws are the effects in
either coordinate. approx = "hsgp" ignores the argument: its basis
coefficients are already that coordinate.
MultiscaleGPFieldT is the reader and MultiscaleGPView::accumulate the
adjoint, one of each rather than a branch at each of the six sites that read
the field -- the templated density, the three gradient functions that write
the block, and the sample-storage transform. Both scales enter the linear
predictor through the same sum, so the likelihood's derivative with respect
to either is the one vector, and the adjoint carries it onto whichever
coordinate is being sampled. multiscale_gp_log_lik_t() built its own copy
of the two GPData views beside make_msgp_gp_views(); it now calls it.
msgp_nc joins the finite-difference sweep in all four gradient modes. A
finite difference cannot see a transform that builds the wrong field -- the
density would stay self-consistent and its gradient would still match its own
difference quotient -- so cpp_msgp_nc_coordinate_spread() reads the
invariant that does: at fixed hyperparameters the two densities differ by
log|det L|, one number for every z. It measures 2e-13.
A collapsed workspace no longer reads a previous model's structure (#76).
The collapsed GP's per-thread workspace caches the NNGP coefficients for the
last (sigma2, phi), the location-to-observation map, and a Newton mode the
next call warm-starts from. Its cache key was the number of GP locations, so
a thread that had evaluated one model and was then handed another with the
same number of locations kept the first model's coefficients and map. Two
collapsed-GP fits in one session over a field of the same size was enough:
the second fit's chain threads started from the first fit's structure, and
the same call with the same seed depended on what had run before it.
CollapsedICARWorkspace was keyed the same way on (S, bym2, N).
Both workspaces now carry the id of the ModelData their structure was built
for and drop every cached quantity when handed a different one. The Newton
warm start goes with the rest: the mode of the previous leapfrog step is
worth starting from, another model's mode is not, and the search is capped at
20 iterations, so a bad enough start is not guaranteed to recover.
cpp_gradient_race() gains prime_other_model, which hands the calling
thread a second model of the same dimensions at the same point before the
race, and clears the Newton warm start before every evaluation so a collapsed
gradient is a function of the point alone. The probe read 4.40e-01 for
gp_collapsed and 8.49e-01 for gp_collapsed_st4 against the old key and
reads 0 against the new one; the eight collapsed fields are now swept by
test-omp-gradient-race.R alongside every other field.
compute_log_post is compute_log_post_impl<double> (#28). The two
densities returned the same number on every structure the templated one
expressed, so keeping them apart bought nothing and cost the standing risk
every bug in #11, #14 and #18 came from: a term in one copy and not the
other, invisible because the runtime gradient check differences whichever
density the active mode differentiates. The analytic body is deleted; what
used to separate the two now sits inside the template.
The structures no autodiff scalar can express -- the collapsed ICAR/BYM2 and
collapsed GP marginals, which locate a mode with Newton and add a Laplace
correction, and proper CAR's dense log-determinant -- are if constexpr
arms, so log_post_impl_gap() names exactly the blocks the other
instantiations leave out, as before. The observation loop is one lambda,
summed by an OpenMP reduction in the double instantiation and serially in
the others, so the H path keeps its parallel likelihood without a second copy
of the 250-line body. skip_obs_loop and the temporal GP's precomputed block
prior become a LogPostOptions struct defaulted to the full density;
precomputed_st_log_prior goes, having never had a caller.
Writing the collapsed fields into the template found the guard they needed: the spatial hyperparameter priors were gated on the field having a parameter slot, so a marginalized field would have dropped the prior on its own scale. Nothing reached that -- the gap kept collapsed models out of the templated density entirely -- but the merge could not have kept it.
The double twins the deleted body was the last caller of go with it:
gp_nngp_log_lik, multiscale_gp_log_lik, log_prior_sigma2_pc,
log_prior_phi_uniform and gp_gradient_w in hmc_gp.h, hsgp_evaluate
and hsgp_log_prior_beta in hmc_hsgp.h, ar1_nc_log_prior in
hmc_temporal.h and log_prior_temporal_sigma2_pc in hmc_temporal_gp.h,
each of which had a templated counterpart that is now the only definition.
So do compute_gradient_gp_autodiff, compute_gradient_msgp_autodiff and
compute_gradient_gp_temporal_autodiff, three hand-written ad::Var
transcriptions of the density that resolve_gradient_fn never returned and
nothing else named -- a copy no test could reach and no assertion could hold
to the original. 2260 lines net.
What the harness asserts changes with it: comparing compute_log_post to
compute_log_post_impl<double> is now comparing a function with itself, so
that check is replaced by two that still have content. The density must not
depend on how many threads sum it, which is what the new n_threads argument
to cpp_gradient_check reaches; and the value each mode reports must be the
density its gradient describes, which holds the fused skip_obs_loop path
and the arena::Var instantiation to the double one.
Across twenty-two structures every entry of the analytic gradient is bit-identical to before, and the density moves by at most 3.7e-16 relative -- one to two units in the last place, from summing the same terms in the template's order. Wall-clock on a 4000-row binomial fit, 400 iterations, before and after: plain GLM 0.59 s / 0.60 s, ICAR + RW1 16.5 s / 15.7 s, collapsed ICAR + RW1 320 s / 281 s, crossed random effects 46.7 s / 44.1 s.
test-temporal.R's AR1 fit asserts a correlation its fixture identifies.
It read mean(rho_samples) > 0 on ten time points, where the posterior for
rho is nearly its prior: sd 0.32 around a mean of 0.07, with a 90% interval
of [-0.54, +0.69] even at 8000 iterations. The sign of the mean over 100
retained draws sat inside its own Monte Carlo error, and across five seeds it
came out anywhere from -0.03 to +0.40, so what the assertion turned on was
which trajectory the sampler happened to take. The fixture now carries eighty
time points, where the posterior mean is +0.47 or better across eight seeds
against 0 for a chain that left rho at its prior, and the threshold says so.
The strengthened fixture passes on the code either side of #28.
The templated log posterior expresses the HSGP multi-scale GP and the
non-centred GP (#26). These were the two structures log_post_impl_gap()
named as expressible but unwritten, and both were reachable from R:
spatial_multiscale(approx = "hsgp") and
spatial_gp(parameterization = "noncentered"), the latter the default
coordinate for a GP fit. Declared rather than written, they cost the
A_r / A / A_t gradient modes those two models, since the front door
refuses a mode whose density is missing a block. The multi-scale case was a
missing block of the kind #18 fixed for the plain HSGP field -- the
parameters would carry no prior and contribute nothing to the linear
predictor -- and the non-centred case was worse than missing: the templated
density read the parameter slot as w when it holds z, so it evaluated a
different model rather than an incomplete one.
Both are now written against the same primitives the analytic density uses,
rather than transcribed beside them. hsgp_evaluate_t() is one spectral
expansion for every scalar type, with the Eigen matvec as its double arm,
and the plain HSGP field, the SVC terms and both multi-scale scales now read
it instead of carrying three inline copies of f = Phi (sqrt(S) . beta).
NNGPNCWorkspace and nngp_nc_forward() are templated, so the z -> w
recursion has one definition rather than one per scalar type, and the
conditional-variance floor is ratiod_cov::nngp_floor_cond_var on both.
Measured against central differences of compute_log_post, all three
autodiff modes now agree to 5.7e-09 relative on gp_nc and 2.0e-09 on
msgp_hsgp, with no identically-zero gradient entry, and the two densities
agree to 5e-13. Both fields move from the gap list into the harness's kernel
sweep, where they are checked under handcoded and arena and included in
the density-equality and fused-value assertions.
temporal_iid() reaches the standalone temporal field, and
temporal_ar1(parameterization = "noncentered") reaches the AR1 one (#40).
TemporalType::IID had a branch in temporal_log_prior() and in its
templated twin, but no R constructor and no "iid" case in the C++ parser,
so nothing could select it; the H-mode and A/A_r branches had no IID case at
all and would have contributed zero. A complete non-centred AR1
reparameterization -- ar1_nc_forward, ar1_nc_log_prior, ar1_nc_gradient
-- sat in hmc_temporal.h with no callers anywhere. Both are now wired:
temporal_effects() reconstructs the effects wherever a density or a
gradient reads the block, and temporal_gmrf_prior_grad() writes the
gradient in whichever coordinate the sampler moves in. Wiring
ar1_nc_gradient up showed it was wrong: its log_tau and logit_rho sums
paired the backward adjoint adj[t], which carries the AR1 recursion
backward, with a forward derivative dphi[t]/dtheta that already carries the
same recursion, counting it twice. Against central differences of the log
posterior the two hyperparameters were off by 12% and 8% while all ten field
entries matched to 4e-10; the correct pairing is the direct partial
grad_phi_lik[t], which the function's own T == 1 branch already used.
All four gradient modes now agree with finite differences on temporal_iid
and temporal_ar1_nc, both added to test_gradient_check.cpp.
spatial_multiscale() under mode = "laplace" fits instead of reporting a
missing group_var (#39). cpp_laplace_fit_multiscale_gp was a registered
export with no caller: fit_laplace_spatial() had no "multiscale" branch,
so the model fell through to the areal path and failed on a message about
areal grouping. spatial_gp() failed the same way -- the areal preparation
ran ahead of the GP dispatch -- and would have failed again on
result$hessian, which cpp_laplace_fit_gp never returned. The two
coordinate-based fields now dispatch before the areal preparation, and both
run through one mode-finder: laplace_mode_nngp() takes a list of NNGP
blocks, so a single-scale GP passes one and a multi-scale GP passes two,
replacing two near-identical Newton loops. Three defects went with the
duplication: the multi-scale copy ignored cov_type and hard-coded an
exponential covariance, and wrote a conditional mean from raw covariances
rather than the Kriging weights; both copies assembled only the diagonal of
the NNGP precision, so Newton converged somewhere other than the mode the
reported log-marginal describes; and both mapped observation i to location
i, dropping the spatial term entirely for every observation past the last
unique location. obs_to_loc is now passed in, the Hessian is assembled at
the converged point and returned, and hsgp and svc name themselves in
BACKEND_STRUCTURE_SUPPORT rather than reaching a dispatch that would fit
the model without the field.
The Laplace NNGP path also read R's 1-based nn_order as 0-based, writing one
past the end of the field block on the last location. Nothing reached it
before; it now converts at the call site, as the sampler entry does, and
add_nngp_block_laplace() refuses an ordering outside 0..n-1 rather than
corrupting the heap.
Multi-scale HSGP is reachable (#38). data.msgp_is_hsgp selected a
complete C++ branch that swapped the two neighbour sets for one shared
Hilbert-space basis, and no R function could set it -- spatial_multiscale()
had no approx argument, and the internal validator the dispatch resolved to
raised an error saying so. spatial_multiscale(approx = "hsgp", m =, c_boundary =) now selects it, the neighbour computation is skipped on that
path, and msgp_hsgp joins the finite-difference harness. On 90 observations
with m = 8 it fits in 3.2 s against 14.3 s for the neighbour sets. The
templated log posterior still declares the gap it declared before, and now
also skips the block rather than reading neighbour sets sized for the
locations against a parameter block sized for the basis.
spatial_multiscale(sampler =) offers the two strategies that exist
(#36). Five of its seven names -- "noncentered", "centered",
"interweaved", "adaptive", "riemannian" -- produced byte-identical
sampling: the field is centred at every site that reads it for the density or
the gradient, and the only enumerator read anywhere selects L-BFGS warmup
adaptation. The argument now takes "auto" and "lbfgs", which is what it
selects between, and the enum's documented default no longer claims a
parameterization the code does not carry. A non-centred coordinate for the
multi-scale field is filed separately.
make_msgp_gp_views() existed twice, and the copy feeding the gradients left
nn_neighbor_dist empty. Its current callers do not read that field, but the
non-centred NNGP helpers do, so the two are now one builder that fills every
field GPData declares.
spatial_gp(solver =, cg_tol =, cg_maxiter =) are gone (#37). The CG /
PCG / GPU dispatch in hmc_gp_cg.h was complete and had no callers -- the
header was not #included anywhere, so it was never compiled. Routing the
sampler through it would have made the package slower: the systems it solves
are the k x k neighbour covariances, and CG needs about k iterations on a
dense covariance rather than converging early. Measured against the Cholesky
already running, on exponential-covariance neighbourhoods in the unit square,
10,000 solves each: k=5 3.1x slower, k=10 2.3x, k=15 (the nn default) 3.1x,
k=30 (the nn_regional default) 2.4x, k=50 2.8x, returning a 1e-6 solution
rather than a machine-accurate one. The N-scaling of an NNGP is O(N k^3)
and the inner solver does not touch N, so "cg"'s documented case -- "better
for large N (> 2000)" -- does not hold. "gpu" had no path either: nothing in
the package builds CUDA, and its own caller commented the GPU branch out. The
three arguments and the dead header are removed; gpu_available() and
gpu_info() are unaffected.
An areal spatial field paired with temporal_multiscale() now gets its
multi-scale block (#74). cpp_hmc_fit set has_multiscale_temporal from
the temporal bundle and then cleared it about a hundred lines later, among
the feature flags that belong to cpp_hmc_fit_gp. compute_param_layout()
reads that field to allocate the trend, seasonal and short-term blocks, so it
allocated none of them while R went on naming their columns. The combination
is reachable because use_gp_sampler is
is_gp_spatial(spatial) || (is_multiscale_temporal(temporal) && !has_areal_spatial):
a multi-scale temporal term with an areal spatial field is the one case that
routes to cpp_hmc_fit. Measured on num | denom ~ x with
spatial_car(level = "group") and
temporal_multiscale("year", trend = "rw1", short_term = "ar1"), 6 units by
12 times, poisson_gamma, 100 iterations: R named 39 columns, 27 of them the
multi-scale block's, and those 27 were not a posterior -- sigma2_trend and
sigma2_short spanned [0, Inf] with a standard deviation of NaN,
rho_short reached exactly -1 and 1 though it is logit-mapped onto the open
interval, and the field columns had standard deviations in the thousands.
All 27 now respect their own priors. test-temporal-multiscale.R checks the
block against its prior support the way test-gp-sampler-blocks.R does, and
fails 29 assertions without the fix.
The GP sampler entry carries the random-effect block it was handed (#72).
cpp_hmc_fit_gp set n_re_terms = 0 and has_re_slopes = false
unconditionally, and the R backend handed it re_group / n_re_groups
alone rather than the re_params bundle cpp_hmc_fit receives.
prepare_hmc_data() and hmc_param_layout() above it parsed and named the
slope and multi-term structure regardless, so the two layouts disagreed on
the size of the RE block and every column after it was read at the wrong
offset. Measured on num | denom ~ x + (1 + x | g) with
spatial_gp(~ lon + lat, nn = 3): R named 28 columns, the sampler wrote the
legacy block instead, and the column labelled phi_gp ran to 2.37e+149
against a uniform prior on (0.01, 100). It now reports 4.08 to 98.10.
Affects random slopes (correlated and uncorrelated) and crossed / nested
multi-term random effects on spatial_gp(), spatial_multiscale(),
spatial_hsgp(), and on a multiscale temporal model with no areal spatial
term.
The AR1 temporal rho prior reaches the GP sampler entry (#73). That
entry's temporal block read type, time_idx, group_idx, n_times,
n_groups, n_params, cyclic, shared, tau_shape and tau_rate and
stopped there, so the rho_prior_a / rho_prior_b the R backend put on the
list crossed the .Call boundary and were never assigned; the fit ran on
ModelData's default anchors. Measured on num | denom ~ x with
temporal_ar1("year"), 6 groups x 10 times: the posterior mean of rho_ar1
under prior_beta(2, 2) and prior_beta(50, 1) was bit-identical on a GP
main effect and moved from -0.071 to 0.972 on an areal one. The GP arm now
moves -0.103 to 0.933.
Both bundles are assembled by one builder each entry point calls
(apply_re_params / apply_temporal_params, src/hmc_model_data_blocks.h),
the way #70 did for the spatiotemporal, latent, TVC and SVC bundles, and the
R backend builds them once above the branch. A field read by one extraction
and not the other is what both of these were.
The random-effect gradient is written once (#72). With the block
reaching the sampler, resolve_gradient_fn() has to keep such a model off
the specialized functions, which read the legacy single-term block: crossed
and nested terms now raise their own feature bit (GF_RE_MULTI) alongside
GF_RE_SLOPES, and both fall through to compute_gradient_composite. That
function carried a partial copy of the block -- no Cholesky factor, no LKJ
prior, no re = diag(sigma) L z transform -- so a correlated slope beside a
GP field reached it with the slope scale's gradient at zero against a
finite difference of 0.409, and the runtime gradient check fell back to
numerical. The prior phase and the chain rule are now
re_blocks_prior_grad() / re_blocks_writeback()
(src/hmc_re_blocks_grad.h), called by compute_gradient_analytical() and
by the composite alike, with re_blocks_eta() / re_blocks_scatter() for
the observation loop. Three new finite-difference fixtures cover what could
not be built before: gp_slopes_corr, gp_crossed, and gp_slopes at the
composite it now dispatches to.
A GP spatial main effect carries the blocks beside it (#70).
is_gp_spatial() routes the whole model to cpp_hmc_fit_gp_v2, which was
handed gp_params, ms_gp_params, ms_temporal_params, rsr_params,
temporal_params and the ZI arguments and nothing else. A spatiotemporal
interaction, a latent factor, a TVC or an SVC paired with spatial_gp(),
spatial_multiscale() or spatial_hsgp() was dropped from the ModelData
the sampler ran on, while initialize_hmc_params_full() and
hmc_param_layout() above it still allocated and named that block's
parameters: every column from the dropped block onward was read at the wrong
offset and reported under the wrong name. Measured on 6 units x 5 times with
spatial_gp(~ lon + lat) and spatiotemporal(type = "separable"),
phi_st_space came back on (22.7, 43380.7) against a uniform prior on
(0.01, 10) whose density is -Inf outside those bounds. The four bundles
are built once, above the sampler branch, and both entry points assemble
their blocks through one set of builders (src/hmc_model_data_blocks.h)
rather than each writing its own extraction. spatiotemporal() no longer
refuses a GP main effect.
A specialized gradient declares the blocks it writes (#71).
resolve_gradient_fn() selected a hand-coded gradient by a chain of
negations, one per feature per function, and several of them omitted
has_spatiotemporal, so a model carrying an interaction reached a function
that never writes it. Each entry now names the blocks it handles and the
model's own blocks are one mask, so a model carrying anything outside that
set falls through to the composite; adding a block to GradFeature is a
block every entry that does not name it falls through on. Three combinations
the negations could not express fall out of it: a GP margin and a GMRF margin
are different temporal blocks, a multiscale temporal term is not an SVC term,
and no specialized function writes the random-slope block.
The composite gradient is a catch-all for every sampled field. It had no
NNGP arm at all -- neither the GP main effect, the multi-scale GP, nor the
NNGP approximation of an SVC -- so a model routed to it left those blocks at
zero while compute_log_post() carried them into eta. All three are written
from the same ratiod_gp / ratiod_svc entry points the dedicated kernels
use. A collapsed field alongside a second block routes to the numerical
gradient of the density, its marginal being an inner Laplace at a mode that
moves with every other block in eta.
Twelve cpp_gradient_check() fields cover the combinations, each held to
central differences at 1e-8 to 1e-7, and cpp_gradient_dispatch() reports the
function a model resolves to so the mask is asserted directly rather than
inferred from whether a gradient happens to come out right.
test-gp-sampler-blocks.R fits three GP-plus-block models end to end and
reads the draws against the prior support the densities refuse outside of.
The GP spatiotemporal interaction exists (#68, #69). type = "separable"
reached the sampler and segfaulted on the first log-posterior evaluation, in
every gradient mode. SpatiotemporalData's NNGP members -- nn, nn_idx,
nn_order, nn_dist_space, nn_dist_time, coords, time_values,
nonsep_type, cov_space, cov_time -- were never assigned by anything:
st_params did not carry them and the sampler entry did not write them, so
st_gp_nngp_log_lik() read nn_order[0] on an empty vector, and the four
scalars among them were indeterminate rather than unset. Both types are now
built end to end. The interaction carries its own coords, so an areal main
effect pairs with a continuous interaction over the same units; nonsep_type,
cov_space, cov_time and nn are its own too, and type = "nonsep_gp" is
reachable from the front door for the first time.
That interaction now has a gradient (#68). compute_gradient_composite
and compute_gradient_spatiotemporal_handcoded both skipped the GP types
entirely, eta included, so a fit would have moved the block under nothing but
its sum-to-zero penalty and computed every other block's gradient at a linear
predictor missing the interaction. st_gp_nngp_grad() (src/st_prior_grad.h)
writes the field, its precision and both ranges, sharing the conditional
decomposition with the density through st_gp_nngp_scan() so the two cannot
disagree about what the prior is. sigma2 factors out of the neighbour block,
which collapses its gradient to a closed form with no extra solve; the two
ranges take one further pair of triangular solves against the factor the scan
already built. Checked against central differences at 1e-7 to 1e-5 on four
new cpp_gradient_check() fields (stgp, stgp_matern, stgp_gneiting,
stgp_latent), on both the analytic and the templated density.
The interaction's covariance is one definition (#68). The GP branch
carried its own Cholesky, its own two triangular solves and its own
conditional-variance floor, a fifth copy of what ratiod_cov (src/hmc_cov.h)
holds for every other NNGP path; it now reads that one. Each arm is written as
a correlation with r(0, 0) = 1 and scaled by sigma2 once, so the marginal
variance is what the PC prior is a prior on and the neighbour block's diagonal
is the kernel's own value rather than a constant written beside it.
cpp_test_st_gp_log_lik() drives the shipped density against a dense
multivariate normal built from the kernel formulas, at nn = n - 1 where the
NNGP is exact.
Two non-separability names are gone, for different reasons (#68).
Cressie-Huang (1999) had no branch of its own and was evaluated as the product
covariance. The additive kernel C_s(h) + C_t(u) cannot be an interaction
covariance on a complete grid: its rank there is S + T - 1 of S * T, and
every direction it drops is an interaction -- measured on a 6 x 5 grid, rank
10 of 30 and smallest eigenvalue -1.6e-15, against condition numbers of 8.4
and 10.8 for the product and Gneiting arms.
spatiotemporal_gp() is withdrawn (#68). It emitted type = "st_gp", a
string no other R file read and the sampler entry rejected, and its
observation-level index set had no relationship to the S x T sum-to-zero
penalty the density applies. spatiotemporal(type = "nonsep_gp", nonsep_type =, coords =) is the spelling that works, and carries the same
knobs.
A spatial GP indexed an interaction by row (#68). build_st_index() set
s_idx <- seq_len(N) for a GP spatial component while n_spatial counted
unique locations, so st_flat ran past the S x T grid it indexes the moment
two observations shared a location. It reads obs_to_loc.
compute_gradient_latent_handcoded was selected for models it cannot write
(#71). Its guard checked for no spatial and no temporal main effect and not
for an interaction, so a latent factor alongside one returned the latent
factors at -2.34 against a finite difference of 0.114. The remaining
specialized guards omitting the same check are listed in #71.
The spatiotemporal interaction's AR1 time margin exists (#66).
spatiotemporal(temporal = temporal_ar1(...)) allocated a correlation and
put a Beta prior on it, and nothing else read it:
spatiotemporal_log_prior() opened with (void)rho, type_ii_log_prior()
carried a // AR1 would need rho parameter comment and no branch, and
type_iv_log_prior() fell through both stencils, so its Kronecker quadratic
form was exactly zero and its rank_time read rw2_rank. The field had no
temporal prior curvature at all, and the density and its gradient disagreed
on the normalizer. The AR1 precision R(rho) is now the interaction's time
margin on Type II, Type IV and HSGP-ST -- the quadratic form and the
rho-dependent log|R| = log(1 - rho^2), differentiated together -- and
ratiod_ar1::ar1_precision_apply() / ar1_quadratic_form() /
ar1_log_prior() (src/ar1_shared.h) are the one matrix the scalar AR1
density and the Kronecker margin both read. Measured at the gradient check's
own draw, logit_rho_st went from an analytic 0.000000 against a
finite-difference -0.299363 -- the prior alone, since the density had no
other dependence -- to -9.133676 against -9.133675.
One spatiotemporal interaction gradient (#66). The two hand-coded
gradients that reach an interaction each carried their own copy of every
type. They had drifted: compute_gradient_composite() had no Type II RW2
branch at all and normalized that arm at rank T instead of T - 2, no
AR1 branch anywhere, and neither wrote logit_rho_st.
st_interaction_prior_grad() (src/st_prior_grad.h) is the single
implementation both now call, so a term added to the density reaches both
gradients or neither. st_time_rank() is the matching single rank, read by
the density, the non-centered rank correction and both gradients.
Type I and Type III no longer allocate a correlation. Neither reads a time margin -- Type I is iid over the whole grid, Type III's time margin is unstructured -- so the parameter the layout used to allocate for them moved under its own prior and nothing else. A time margin an interaction reads but cannot express (a GP margin under Type IV) is now refused at the door rather than silently run as RW1.
temporal_ar1()'s correlation lives on (-1, 1) (#67). It was mapped onto
(0, 1), so a negative autocorrelation was not reachable from that block at
all, while the TVC arm, the multi-scale short-term arm and the spatiotemporal
interaction all used rho = 2u - 1. ?temporal_ar1 documented
Uniform(-1, 1) and ?priors_default documented Beta(2, 2), and the code did
neither. The mapping is now 2u - 1 with the shared Beta(2, 2) on
u = (rho + 1) / 2, which changes the posterior of every existing
temporal_ar1() fit. ratiod_ar1::rho_from_logit() and drho_dlogit()
are the one map and the one chain rule; the eight hand-written copies of each
are gone.
temporal_ar1(rho_prior = ) is wired (#67). It was stored on the spec
and never read again, as was ratiod_priors(rho_temporal = ). Both now reach
the density and every gradient as the Beta anchors on u = (rho + 1) / 2,
the block's own prior taking precedence; a spatiotemporal() interaction's
AR1 margin reads its own temporal spec's prior the same way. A rho_prior
that is not a prior_beta() is refused at the door -- which includes
temporal_ar1(rho = 0.8), where rho used to partial-match rho_prior and
store a bare number nothing read.
One AR1 correlation prior, read by every density and every gradient
(#58). The TVC rho carried a Uniform(-1, 1) prior in compute_log_post and
a Beta(2, 2) one in compute_log_post_impl, exactly one log u + log(1 - u)
apart, so the autodiff modes sampled a different posterior from the handcoded
one -- 1.408954 in the density and the whole gap in the gradient at
logit_rho_tvc. The prior is Beta(2, 2) on u = (rho + 1) / 2 everywhere it
applies now, which is what priors_default() documents for a temporal
correlation and what the multi-scale short-term arm already used; the
spatiotemporal interaction's rho moves onto it from Uniform(-1, 1) as well, so
a fit carrying one shifts by that term. ratiod_ar1::log_prior_logit_rho()
and its gradient (src/ar1_shared.h) are the single expression all four
correlations read -- the two densities, both handcoded gradient functions and
the multi-scale arm -- with each site naming its own (a, b), so the
Uniform(0, 1) the plain temporal_ar1() block carries is the same helper at
(1, 1) rather than a second spelling. Four unreachable copies of the prior
in hmc_tvc.h, hmc_tvc_autodiff.h and hmc_temporal_multiscale.h are gone.
tvc_ar1 joins the structures in test-gradient-correctness.R, whose
density-equality assertion it now answers.
One floored 1 - rho^2 (#59). It was clamped at 1e-10 in five places,
clamped at 1e-12 in three, added as + 1e-10 in five more, and left
unfloored in eight -- including the AR1 marginal variance, its templated twin,
the temporal GP's two conditional variances and the Laplace path. + 1e-10
and max(., 1e-10) are different functions, so paths that must agree parted
company wherever either bound: at a correlation near the edge of its range the
two densities were 3.14 apart on temporal_ar1, 27.2 on tvc_ar1 and 1.9e-03
on ms_temporal, and one step further out compute_log_post returned -Inf
with a non-finite handcoded gradient while the templated density stayed
finite. ratiod_ar1::one_minus_rho2() is the one clamp every site reads.
cpp_gradient_check(near_unit_rho = TRUE) places each sampled correlation
where the floor binds, which nothing in the sweep reached before.
Two defects found while testing the above are filed rather than folded in:
no gradient function writes logit_rho_st, so a Type IV interaction with an
AR1 time margin moves its rho under a flat gradient and fuses a log posterior
missing the prior (#66); and temporal_ar1()'s rho lives on (0, 1) while its
documentation describes (-1, 1), with rho_prior inert (#67).
PG GP fits carry their spatial field, so they can predict one (#63).
convert_pg_gp_to_ratiod_fit() read the field draws off "w_gp" while the
sampler returns them under "gp", and [[ on a missing name is NULL, so
.internal$w_gp was NULL on every such fit and predict_spatial_gp_pg()
returned nothing without saying so. fitted() was unaffected -- it reads the
linear predictor the sampler stores directly. The (sigma2, phi) draws that
go with the field were also read from the first chain alone while the field
spanned every chain, so a multi-chain fit paired half its draws with nothing;
both are stacked in the same chain order now and the kriging entry refuses a
pair that is not row-aligned.
A GP prediction is made in the units its fit was made in.
spatial_gp(scale_coords = TRUE), the default, centres and scales the
coordinates before the neighbour structure is built, and scale() puts the
centre and scale on an attribute that subsetting the matrix drops. Nothing
recorded them, so coords.0 was read raw against a scaled training set on
every backend: kriging the training locations of a fit back onto themselves
recovered its own field at correlation -0.01. gp_scale_new_coords() is the
conversion, and the same fixture now reads 0.95.
One kriging implementation, and one covariance family behind it.
kriging_predict() and cov_function() were an R triple loop over
(draw, location, neighbour pair) and an R copy of the covariance family that
knew three of the four kernels, resolving "spherical" to an exponential one
-- the same defect fixed in the compiled family at 267f2cb, in the other copy.
cpp_kriging_predict() replaces both: it reads hmc_cov.h, the family every
fitting path reads, and takes the whole draw matrix at once, so the neighbour
sets are found once rather than per draw. The HMC and PG prediction paths both
go through it.
cov_type_code() is the one name-to-code mapping. Three copies of it in R
each ended in a bare 0L, so a covariance name none of them knew ran
exponential; one of the three was missing "spherical" outright. The C++ side
meets it at ratiod_cov::cov_type_from_int() in cov_type_code.h, which the
PG entries now share, and both ends refuse a value they do not know.
The PG binomial GP sampler now pairs observations with locations by
coordinate rather than by row position (#60). spatial_gp() reduces the
data to unique locations and builds the NNGP on those, but
cpp_pg_binomial_gibbs_gp() was never given the map back: it aggregated the
per-location Polya-Gamma likelihood over row position and discarded every row
past n_spatial, and it wrote the field into only the first n_spatial
linear predictors. On the four-observations-per-location design used to
measure it, 75 of 100 observations reached no location at all, the 25 that did
were credited to locations by row order, and the field never entered three
quarters of eta. obs_to_loc is now passed and read on both sides. The
coordinates handed to the sampler were also the observation-order matrix
rather than the unique locations the neighbour structure indexes, so the
neighbour blocks were assembled from the wrong points; unique_coords is what
goes in now. Recovery of the true field over three seeds: correlation 0.95,
0.98, 0.98.
sigma2_gp and phi_gp are drawn from the NNGP density the field was
sampled under (#62). The sigma2 acceptance ratio was built from an
independent N(0, sigma2) quadratic form carrying no log-determinant, so
raising sigma2 only ever reduced the penalty and the ratio pointed one way
at every value; phi's ratio was the log-scale Jacobian alone, so phi was
drawn from its prior rather than its posterior. Both now difference
sum_k log N(w_k; mu_k, d_k) over the sampler's own conditioning order,
evaluated by the same conditional the field sweep draws from. The sigma2
prior's change of variables was also carrying a full log(sigma2) where the
PC prior on sigma = sqrt(sigma2) gives +0.5 log(sigma2) once the walk's
Jacobian is included. Against a true sigma2 of 0.7, posterior medians over
three seeds are 0.36, 0.60 and 0.65, where the issue measured [94, 1746] and
[24.9, 11568].
spatial_multiscale(cov = ) is read by the multiscale PG sampler (#61).
cov_type was accepted by cpp_pg_binomial_gibbs_multiscale_gp() and never
mentioned again; an exponential kernel was written out inline at six sites, so
three of the four documented choices silently ran the fourth. Both scales now
go through the same pg_nngp_conditional() the single-scale entry uses. That
also replaces what stood beside the kernel: the conditional mean was a raw
covariance-weighted sum of the neighbours in place of the NNGP weights
C^-1 c, the conditional precision was 1 / sigma2 in place of the
conditional variance, nn_idx was indexed by location and its values read as
locations where both are positions in the conditioning order, and each
scale's sigma2 was a conjugate draw under an independent N(0, sigma2)
form that read the PC prior's upper bound as an inverse-gamma rate and its
tail probability not at all.
One convention, checked at the entry. PgNngpGraph carries the neighbour
structure the two entries share, pg_check_nngp_graph() refuses a nn_order
that is not a 0-based permutation of the locations, a nn_idx entry that is
not an earlier position in the order, and a coordinate matrix that is not one
row per location, and pg_check_obs_to_loc() does the same for the
observation map. The 1-based nn_order fixed in 267f2cb was a silent
out-of-bounds read; it is now an error naming the argument.
One walk of the sampler's parameter vector now feeds fitted(),
predict(), ratio(), the draws matrix and the per-structure extractors
(#44). There were four copies of that walk. Three of them were thinner than
the fourth: on the default non-centred parameterization they added the raw
z to the linear predictor instead of sigma_re * z, so every fitted value
and every prediction on a model with random effects was off by the posterior
draw of sigma_re; they recognised one intercept-only random-effect term, so
a crossed or slope-carrying model read parameters that were not random
effects at all; they recognised icar and bym2 among the spatial types and
dropped every other field from the predictor without saying so; and they
applied the logistic link to every family, reporting a number in (0, 1) where
a count ratio belongs. fitted(fit, component = "ratio") and ratio(fit)
now return the same draws, and so does predict(fit, newdata = fit$data).
The block order the R side reads is the order the sampler writes (#44).
hmc_param_layout() (R/hmc_unpack.R) mirrors compute_param_layout() in
src/hmc_sampler.cpp block for block, and is the only place an offset is
computed. The previous walk placed SVC before the areal field, the GP before
the temporal block, TVC before zero-inflation and HSGP among the spatial
types, where the sampler lays them out in a different order; a fit combining
two of those structures had its columns named after the wrong parameters from
the first mismatch onward. It also walked a temporal block for a multiscale
temporal model, which the sampler allocates none for. Every fit now checks
that its layout accounts for exactly the columns the sampler returned, and
errors rather than mislabelling them.
svc(), tvc(), temporal() and spatiotemporal_effects() return their
draws instead of erroring on every fit (#46). No backend stored the draws
they read. The converter now stores each structure's own draws, shaped the
way its extractor indexes them, and the svc() / tvc() error text names the
argument that routes them (spatial = spatial_svc(...),
temporal = temporal_tvc(...)) rather than an argument tratio() does not
have. The four \dontrun{} examples that showed the wrong spelling are
fixed. tvc() gained a group argument, since the draws carry one
trajectory per grouping level.
pp_check() plots a fit (#48). It read object$draws$y_num_rep, and
object$draws is a matrix on every backend, so the call raised
$ operator is invalid for atomic vectors on its second executable line and
never reached the guard written to explain the missing draws. The replicates
are now drawn rather than looked up: posterior_predict() simulates a
numerator and a denominator response from each posterior draw's fitted means
and that draw's dispersion, through the same family machinery
prior_predict() uses.
as_draws(), pp_check() and posterior_predict() are methods on the
generics that own them, not new generics (#52). They were bare
UseMethod() calls, so attaching tulpaRatio after posterior or bayesplot
masked their generic and a user's as_draws() on a draws_matrix or
pp_check(y, yrep) on plain vectors dispatched into this package and found
no method. posterior::as_draws, bayesplot::pp_check and
tulpa::posterior_predict are imported and re-exported instead; bayesplot
moves from Suggests to Imports.
spread_draws(), gather_draws() and point_interval() are renamed to
draws_wide(), draws_long() and draws_interval() (#52). tidybayes owns
those three names and takes different arguments under them (variable[index]
specifications rather than bare parameter names), so a method on its generic
would not have honoured the contract. The old names are removed.
The offline gradient check can build a model for every kernel the
dispatcher can select (#45). cpp_gradient_check covered the areal fields,
the random-effect layouts, HSGP and TVC; it could not build a GP, multi-scale
GP, SVC, temporal-GP, multiscale-temporal or latent-factor model, so the
handcoded gradients those run on under the default mode had no
finite-difference case. make_model() now builds gp, gp_collapsed,
gp_temporal, msgp, msgp_temporal, svc, svc_hsgp, temporal_gp,
ms_temporal, latent, tvc_ar1 and temporal_ar1, and an unrecognised
field name is an error rather than a model with no structured field in it,
which would pass while testing nothing. The temporal block's AR1 arm carries
a rho that its RW1 and RW2 siblings have no counterpart for, so the two cases
that reached it before covered none of the terms that rho enters.
The runtime gradient check differences away from the origin (#45). Every structured block starts at exactly zero, where a quadratic-form prior contributes nothing to any gradient: a wrong precision matrix, a sign error on it and an absent sum-to-zero penalty are all invisible there. The check now moves each field block off zero by a small deterministic offset first. Both call sites share one function, so a mismatch raises an R-level warning from the multi-chain path as well as the single-chain one.
The PG binomial GP sampler indexed one past the end of every per-location
array, and could take R down with it. backend_pg.R passed nn_order
1-based while the C++ uses its values directly as indices into w, coords
and the per-location likelihood accumulators, all sized n_spatial, so each
could reach n_spatial itself. backend_hmc.R converts the same vector with
- 1L; the PG path did not. The symptom was STATUS_HEAP_CORRUPTION
(0xC0000374) killing the R process, reached at some parameter states and not
others -- 20 locations crashed at nn = 5 and nn = 6 while 25 locations
survived -- and, where it did not crash, the GP effects were attached to the
wrong locations.
spatial_gp(cov = ) runs the kernel you asked for under mode = "pg"
(#35). The selector tested cov_type 0 and 1 and sent everything else to an
unconditional else computing Matern 2.5, so "gaussian" and "spherical"
-- both documented and match.arg-validated -- ran a kernel the user did not
choose. It now routes through the shared ratiod_cov family introduced above,
and an unrecognised code is an error rather than whichever kernel a
fallthrough reaches. The neighbour block goes through the shared
factorization at the same time, replacing a seventh distinct regularization
(a 1e-10 value substituted into the pivot). Measured on a 12-location
fixture, against an NNGP conditional written out independently in R: all four
kernels now agree to machine precision, and at the probed location the
conditional mean moves by 0.11 for "gaussian" and 0.34 for "spherical"
against the Matern 2.5 both used to return.
tratio(mode = "pg") honours control$seed and control$chains.
fit_pg_binomial() accepts both and threads seed down to a per-chain
set.seed(seed + chain - 1), but tratio() called it without either, so PG
fits ran off whatever RNG state they inherited -- two runs of the same call
differed by 1.23 on the intercept and 3.19 on sigma_re -- and always used
the backend default of 4 chains whatever chains said. Runs at a fixed seed
are now bit-identical.
spatial_gp(parameterization = "noncentered") gets its handcoded gradient
back (#57). Two defects made it disagree with the density it reports, and
resolve_gradient_fn sends every GP model at the default gradient mode to
compute_gradient_gp_handcoded, so this was the ordinary path for such a fit
rather than an opt-in kernel. nngp_nc_backward read a C_mat it never
filled -- the comment beside it said it deliberately did not rebuild one --
which made every dC/dphi entry zero and so dropped both the -dC * alpha
term from dalpha and the alpha' dC alpha term from dd_dphi; and the
caller subtracted z after nngp_nc_backward had already seeded grad_z
with -z, applying the N(0, 1) prior twice. Measured on the harness: worst
relative deviation 4.610e-01 at log_phi_gp, flat across difference steps
from 1e-3 to 1e-7, and the residual on every one of the 25 z entries equal
to -z to 1.1e-07 with a fitted slope of 0.99999997. Both now measure
0.000e+00 with a V-shaped sweep.
The fits themselves were slow rather than wrong: the runtime gradient check
introduced in this release catches the log_phi mismatch at its probe point
and falls back to numerical gradients, which is what a non-centred GP fit had
been running on. That fallback no longer fires. Before this release the check
differenced at the origin, where z = 0 makes w = 0 and both defects
contribute exactly nothing, so it saw neither.
One covariance family and one neighbour-block factorization behind every
NNGP path (#42, #55, #56). The GP path, the SVC path and their two
templated twins each carried a copy of the four kernels and of the small
dense Cholesky over a neighbour set, and the copies had drifted apart on four
separate numbers: the Gaussian kernel was written with and without the 0.5 in
its exponent, Matern 3/2 carried sqrt(3) exactly in two of them and rounded
to 1.732050808 in the third, one had no spherical case and fell through to
exponential, and the neighbour block was regularized four different ways
(1e-8 always, 1e-8 only once a pivot had already gone non-positive, a
1e-6 floor substituted into the pivot, 1e-4 always) with the conditional
variance floored three more (1e-10, 1e-6, and a 1e-4 blend keeping 1%
of the gradient). The analytic gradient paths factorized the block a third
time again, with no ridge and a pivot replaced by 1e-5. So the value a fit
reported and the gradient it moved on came from different copies. Measured
before: the two densities were 30 nats apart on the NNGP SVC and 0.086 apart
on a Matern GP, with the deviation flat across difference steps from 1e-3 to
1e-7. Everything now routes through src/hmc_cov.h, templated once for the
double evaluation and the ad::Var / fwd::Dual / arena::Var gradients
alike; across gp, gp_matern, gp_gaussian, gp_spherical, gp_temporal,
msgp, msgp_temporal, svc and svc_hsgp the two densities now agree to
the printed digit and every gradient matches central differences of both to
2e-06 or better.
cov = "gaussian" is exp(-0.5 (d / phi)^2) on every path. The analytic
path evaluated exp(-(d / phi)^2) while taking its phi derivative from the
other form, so it was running on half the derivative of the density it
reported; the templated path already used the 0.5. The 0.5 is what makes
phi a lengthscale here as it is for the other three kernels -- it is the
standard squared-exponential and the nu -> infinity limit of Matern at a
fixed phi -- so the shared kernel keeps it. A phi estimated by a previous
cov = "gaussian" fit under a handcoded gradient mode corresponds to
phi / sqrt(2) now; a fit under A, A_r or A_t is unchanged.
cov = "spherical" is differentiated rather than silently run as
exponential (#42). Neither analytic derivative had a spherical case and the
templated dispatcher had no spherical kernel, so both fell through to the
exponential arm. dk/dphi = 1.5 sigma2 (r - r^3) / phi is written out, which
is why dcov_dphi now takes sigma2 -- it is the one case not expressible
through the covariance value. An unrecognised covariance name is an error at
both string parsers, where the GP one used to map anything it did not
recognise to spherical and the SVC one to exponential.
A TVC AR1 field's log tau gradient had the wrong sign on its stationary
normalizer (#43). The AR1 prior parameterizes the stationary variance as
1 / (tau (1 - rho^2)), so -0.5 log(2 pi var) differentiates to +0.5 in
log tau; ar1_grad_log_tau seeded the sum at -0.5, leaving the reported
gradient short by exactly 1 at every value of every parameter (4.46185
against a finite difference of 5.46185 on the harness fixture). The three
AR1 gradient functions, rw1_grad_w and rw2_grad_w also read past their
field on a block too short to carry an increment or a second difference;
those priors are flat there, and they now return that instead.
The fused log posterior of the multi-scale GP and latent-factor gradients
carries the binomial coefficient. compute_obs_ll() dropped it while
compute_log_post() adds it back, so those two kernels reported a value
7156 nats away from the density their gradient describes on a binomial fit.
NUTS consumes that value as log_prob.
Prediction maps a new data set onto the fit's own grouping levels through the
same construction the formula parser used (as.factor(), or interaction()
for a nested term). It previously matched against unique() order, which
differs from the factor's level order whenever the levels are not in order of
first appearance, and attached the wrong group's effect.
Zero-inflation and hurdle terms on count responses now reach the templated
log posterior, so mode = "ess" and the A/A_r/A_t gradient modes stop
silently fitting a plain Poisson/NegBin (#34).
compute_log_post_impl read data.zi_type only inside its
ModelType::BINOMIAL arm. The three count-response arms --
NEGBIN_NEGBIN, POISSON_GAMMA and NEGBIN_GAMMA -- called the plain
log_lik_negbin/log_lik_poisson whatever zi_type said, while still
allocating and sampling a ZI coefficient against a likelihood that did not
contain it. That density is the one ESS evaluates and the one the three
autodiff gradient modes differentiate, so
tratio(..., family = ratiod_zinegbin(), mode = "ess") (and
ratiod_zipois(), ratiod_hurdle_negbin(), ratiod_hurdle_pois(), or an
explicit zi = zi_negbin() and friends) completed with no error and
returned a zero-inflation logit drawn from its prior. Measured on the
finite-difference harness, the templated density returned the same number
for zi_poisson, hurdle_poisson and no zero-inflation at all, and
disagreed with the analytic density by between 1.7 and 46 log units across
the six family/structure combinations; both now agree to 1e-11, and the
analytic ZI gradient -- which nothing had ever checked on a count response
-- matches central differences to machine precision in the handcoded and
arena modes alike. hmc_zi.h's count-response densities are templated on
the scalar type rather than duplicated, so compute_log_post and
compute_log_post_impl<T> evaluate one implementation. make_model() in
test_gradient_check.cpp gained a zi argument and the negbin_negbin and
negbin_gamma families, since it previously hardcoded ZIType::NONE for
every call site and no test combined mode = "ess" with any ZI family.
temporal_tvc(structure = ...): "iid" unblocked, "gp" now errors
instead of silently downgrading to RW1, and the Gibbs backend now rejects
anything but "rw1" instead of silently ignoring the argument (#32).
TemporalType::IID had a complete log-prior and gradient in both the H and
autodiff paths (tvc_term_log_prior, tvc_prior_gradients_ws, and the
templated density via hmc_tvc_autodiff.h's using declarations), and the
C++ string parser already recognized "iid" -- but temporal_tvc()'s
match.arg only offered "rw1"/"rw2"/"ar1"/"gp", so no R call could ever
reach it. structure = "iid" is now accepted, verified against
test_gradient_check.cpp's finite-difference harness in both H and A_r
modes (0 deviation), and recovers a simulated intercept and slope in
test-recovery-constrained-fields.R. Separately, "gp" was accepted and
documented as a first-class option, but TVCData has no GP fields at all
(no lengthscale, no covariance branch) -- the C++ parser silently mapped
any unrecognized string, "gp" included, to RW1, so
temporal_tvc(..., structure = "gp") ran to completion and silently fit an
RW1 varying coefficient. temporal_tvc() now errors immediately if
structure = "gp" is requested, rather than reaching the sampler at all.
Independently, the Gibbs backend's TVC update
(# ---- 5. Update TVC coefficients via univariate MH with RW1 conditional proposal ----) is hardcoded to RW1 and never read the tvc_structure field
R was already sending it, so mode = "gibbs" silently ran RW1 regardless of
what structure asked for; the Gibbs backend now rejects any TVC structure
but "rw1" with a clear error pointing at mode = "hmc". The C++ side's
own duplicated string-to-enum parser in hmc_sampler.cpp was replaced with
a call to the existing ratiod_tvc::parse_tvc_structure(), which had no
callers before this.
spatial_car(..., proper = TRUE) now fits a proper CAR instead of
silently downgrading to ICAR (#31). The R object already computed rho
bounds and printed "Proper CAR"; a complete C++ implementation
(hmc_car_proper.h) already existed too. Neither was wired to the other --
every backend collapsed type = "car_proper" to "icar" before it reached
C++ (# Handle all ICAR-type models (car, car_proper both use ICAR parameterization)), so rho stayed fixed at 1 and the roxygen example's
claim that summary(fit_car) shows a rho_spatial parameter was false. The
HMC backend now threads rho through as its own parameter (Q(rho) = D - rho*W, Beta(1,1) prior by default, logit-transformed), added to
SpatialType, ParamLayout, and the shared spatial_gmrf_prior_grad/
compute_log_post machinery that compute_gradient_composite and every H
gradient mode already read from. phi's gradient is exact; rho's
log-determinant term has no cheap closed form (it needs trace(Q^-1 W)), so
it is central-differenced on that one scalar and chain-ruled to logit scale,
verified against test_gradient_check.cpp's finite-difference harness (0
deviation across all 29 parameters on the car_proper field). Proper CAR is
full rank for rho < 1, so unlike ICAR/BYM2 it is not hard-centered.
Autodiff gradient modes (A, A_r, A_t) refuse at the front door with a
clear error rather than differentiate a density missing the field's prior
entirely, since the log-determinant is not yet expressed in the templated
density; mode = "hmc" with the default gradient_mode = "auto" (or "H")
is unaffected. Gibbs already rejected car_proper (its type whitelist never
included it); Laplace, PG, VI, SGHMC, and ESS now reject it explicitly too,
instead of silently fitting an ICAR model or -- for the three backends that
call compute_log_post_impl<double> directly -- silently dropping the
spatial prior's contribution to the posterior altogether. Fixing this also
surfaced a second, independent bug in predict()'s areal-prediction path,
which read phi_spatial starting one column early for any type it treated
as "ICAR-shaped" once car_proper genuinely had an extra rho parameter
in front of it.
A spatiotemporal Type IV interaction has a non-centered parameterization,
which fixes chains freezing on data with little or no true interaction
(#24). spatiotemporal(..., type = "IV") reached rhat 1.14 at
iter=1000/warmup=500 on data simulated with no space-time signal, and a 4x
longer run moved that to 2.05 while every chain's own posterior sd shrank.
The interaction's precision tau has no likelihood signal pulling it away
from the boundary in that case, so its posterior drifts to very large
values (tau ~ 4.3e8 was observed at warmup's end); under the centered
parameterization the field's conditional scale shrinks as 1/sqrt(tau),
so exploring both regions needs step sizes that differ by orders of
magnitude, and NUTS answered by pinning the trajectory length at
max_treedepth on essentially every iteration and collapsing the step size
toward its floor -- the signature of a funnel between a variance component
and the field it scales, not of a chain that merely needs more draws. A
finite-difference check of compute_gradient_spatiotemporal_handcoded's
Type IV (Kronecker) block, absent until now, confirmed the analytic
gradient was not the cause. spatiotemporal(..., parameterization = "noncentered") samples a tau-free z and reconstructs the interaction as
z / sqrt(tau), decoupling the two: rhat on the same repro drops to 1.00
at iter=1000 and stays at 1.00 at iter=4000, with stable per-chain sd and
a step size about 10x larger. The non-centered path already existed in the
gradient code but was never reachable from R and was itself unchecked
against finite differences; the warmup-end precision-informed mass-matrix
warm-start specific to Type IV (which reads the raw parameter as the
centered field directly) is now skipped under the non-centered
parameterization, which it does not describe. parameterization still
defaults to "centered", since data that do carry real space-time
structure move tau away from the boundary and may not need it.
The Gibbs BYM2 sampler no longer leaks the spatial field's level into the
intercept unchecked (#15). The per-site update held the structured field
(phi) to sum zero by moving its mean into the intercept every sweep,
unconditionally: no acceptance test weighed that move against the
intercept's own prior. dev_notes/gibbs_prior_recovery.R (now
tests/testthat/test-gibbs-prior-recovery.R) fits with zero binomial trials,
where the likelihood is flat and the posterior must equal the prior exactly;
under a N(0, 0.5^2) intercept prior the sampled sd came out at 1049 instead
of 0.5, a ~2000x inflation. The move now competes for acceptance like every
other update that touches the intercept, and gained a dedicated block-level
move (mirroring the one the unstructured component already had) so the
now-resisted level direction still mixes: BYM2 intercept ESS is 749 at
iter=2000/warmup=1000 where it used to need the same bar ICAR clears at
1000/500. The Gibbs priors already matched compute_log_post's (Half-Cauchy
sigma_total, Uniform(0,1) rho, Gamma tau, threaded through fit_gibbs()
via priors) from prior work on this issue; this closes out the remaining
checklist items, including removing the unused pc_lambda left over from an
earlier tau update.
A per-fit cores budget can now reach within-chain work (#17). For a
single-chain HMC fit, cores used to be capped at chains before the user
ever saw it (min(chains, cpp_get_max_threads())), so a one-chain fit always
got exactly one thread for its gradient and likelihood loops regardless of
machine size, and the regions that do carry a num_threads() clause
(compute_gradient_analytical's fallback path, the linalg_fast.h kernels)
had no way to reach more than one. cores is now a genuine budget split two
ways: concurrent_chains = min(cores, chains) chains run at once, each with
cores %/% concurrent_chains threads of its own. Reaching that width when
chains also run concurrently needed OpenMP nesting enabled
(omp_set_max_active_levels(2)): a region inside a chain body is a nested
one, and OpenMP pins nested teams to one thread by default. Most common model
configurations (plain Poisson-gamma/binomial/negbin with simple random
effects) take a vectorized Eigen gradient path that never reaches these OMP
regions at all; the scalar fallback that does carry them serves zero-
inflation, correlated random slopes, and GP/multiscale spatial models.
The two log posteriors are now the same function (#18). compute_log_post
and the templated compute_log_post_impl<double> that the A_r / A / A_t
gradient modes differentiate returned different numbers on every model and
omitted different structures. For those modes the runtime gradient check
differences that same templated density, so a structure missing from it was
invisible to both the gradient and its own numerical reference: they agreed on
the same wrong value. Differencing the two densities against each other
instead, over every structure the parameter layout allocates for, is what
found the following.
| structure | before | after |
| --- | --- | --- |
| HSGP spatial field | absent: all 11 parameters gradient exactly 0, worst relative deviation 1.00 | 1.3e-7 |
| random slopes, uncorrelated | slope variance absent, 1.00 | 1.2e-9 |
| random slopes, correlated | Cholesky factor and variance absent, 1.08 | 1.0e-8 |
| crossed random intercepts | second term absent: 5 effects and its variance 0, 1.33 | 1.2e-9 |
| TVC precision | a different prior on each side, 0.376 | 4.2e-10 |
| collapsed ICAR / BYM2 | read params[-1], access violation | refused |
The residual constant that remained on every model was the binomial
coefficient: the eta-form likelihood drops it as constant in eta, which is
right for a Laplace inner solve and wrong for a log posterior that is reported
and compared. lp - lp_impl was a flat -7156.27 across all fields and is now
0 to floating-point noise.
A collapsed spatial field refuses an autodiff gradient mode instead of
crashing. The collapsed parameterizations marginalize the field out by
locating its mode with Newton and adding a Laplace correction, which is not a
closed-form function of the parameters, so the templated density cannot
express them and the layout allocates no slot for the field. Reading that
absent slot was an access violation. gradient_mode in "A_r", "A", "A_t"
now errors at the front door naming the parameterization; "H" and "N" carry
the analytic density and are unaffected. The same declaration covers the HSGP
multi-scale GP and the non-centred GP, which are expressible but not written
yet (#26).
The TVC precision takes the prior it is configured with. tvc_tau_shape
and tvc_tau_rate reach the sampler from priors and defaulted to
Gamma(2, 0.5), but two of the three implementations ignored them for a
hard-coded PC prior with P(sigma > 1) = 0.01. The density, the specialized
gradient and the composite gradient now all read the configured
hyperparameters, through one derivative helper.
The random-effect prior is one implementation. hmc_re.h computes the
prior and the effective effects for every layout the sampler allocates -- one
or more crossed terms, intercept-only or with slopes, correlated or not,
centred or non-centred -- and both densities call it, so a term cannot be added
to one and not the other. The TVC quadratic forms, per-term prior and eta
assembly move the same way: templated in hmc_tvc.h and re-exported to
ratiod_tvc_ad, replacing copies that had already drifted on the cyclic
wrap-around edge. The RW1/RW2/AR1 templates take a pointer, so a caller holding
a strided slice passes an offset rather than copying into a vector.
cpp_gradient_check() reaches the structures it could not build before
(hsgp, tvc, re, re_crossed, re_slopes, re_slopes_corr,
icar_collapsed, bym2_collapsed), differences either density via
reference =, and reports log_post_impl and impl_gap alongside the
gradient. test-gradient-correctness.R asserts across every field that the two
densities return the same number, that no block has an identically-zero
gradient under any autodiff mode, and that a collapsed model refuses one.
Its tolerance now allows an absolute 1e-6 before the relative comparison: a
central difference of a log posterior of magnitude 1e3-1e4 resolves a gradient
entry to about 1e-7, and the HSGP basis coefficients run to 6e-4, where 1e-4
relative asked for more than the difference can deliver.
Every intrinsic field is now identified at the engine's constant (#12).
s2z_precision(n) returns 1 / (kappa * n)^2 from sd(sum phi) = kappa * n,
so callers pass a PRECISION. The spatiotemporal and TVC penalties took that
constant as an argument and their call sites filled it with a bare 0.001 --
which is kappa itself, weaker by ~2.5e6 at T = 20, pinning the field's sum
at sd 31.6 instead of 0.02. Each penalty now derives the constant from the
length of the sum it pins, so no call site can supply one. The interaction
margins take their own constants rather than sharing one: a space margin sums
S terms and a time margin T.
SVC keeps its own constant, deliberately. Its penalty is rewritten on the
sum rather than the mean -- the two forms are the same penalty, related by a
factor of n_obs -- and both it and its gradient now read one
svc_sum_to_zero_ridge(). That constant is NOT s2z_precision(n_obs): an
NNGP field's prior is proper, so its mean is identified and what is wanted is
a ridge, not the pin a rank-deficient field needs. Substituting the intrinsic
constant was tried and reverted -- at N = 80 it is 156.25 against the
ridge's 0.0125, and at that stiffness 4 chains x 1000 iterations return in 9
seconds with per-chain posterior SD 0, Rhat 15.1 and the slope at 0.095
against a truth of 0.3, where the ridge samples for minutes. Deriving the
right constant for a proper field on its own terms is #25.
The multiscale trend and seasonal arms are pinned (#12). Neither carried a sum-to-zero term on any path. Both are intrinsic and both enter the same linear predictor, so each had a constant null direction aliased with the intercept and with the other -- a two-dimensional ridge the fitted series hides, because only the level moves. Both are now pinned in the log posterior and in the analytic gradient; the proper short-term arm (AR1/IID) identifies its own level and is left alone. On a 500-observation binomial fit against a known intercept of 0.5, over 4 chains:
| trend + seasonal | intercept | posterior SD | Rhat | slope (truth 0.3) | |---|---|---|---|---| | RW1 + cyclic seasonal, before | 2.6101 | 7.1840 | 1.589 | 0.2886 | | RW1 + cyclic seasonal, after | 0.5054 | 0.0170 | 1.012 | 0.2878 | | RW2, before | 3.7200 | 1.1098 | 2.163 | 0.2926 | | RW2, after | 0.4868 | 0.0167 | 1.000 | 0.2907 |
The slope is right either way, which is why nothing downstream flagged it.
The spatiotemporal interaction's margins are identified (#12). Same fit with a Knorr-Held Type IV interaction over 9 units x 5 times: the intercept posterior SD moves from 0.1106 to 0.0109 and Rhat from 3.43 to 1.14. It does not reach 1.05, and a 4x longer run makes it worse rather than better; that residual is tracked separately as #24.
The cyclic RW rank is no longer overcounted (#12). RW1 and RW2 ranks were
computed inline as cyclic ? T : T-1 / T-2 at 16 sites. A cycle-graph
Laplacian still annihilates the constant, so a cyclic RW1 has rank T-1, and
a cyclic RW2 T-1 as well (a linear ramp is not periodic). Every site now
calls tulpa::rw1_rank() / rw2_rank() from the engine, which cannot drift
from the engine's own normalizers.
The spatiotemporal analytic gradient matches the penalty it differentiates.
Three further 0.001 literals sat in the ST gradient paths, including the one
feeding the fused log-posterior reconstruction. With the old constant the
resulting mismatch was too small to see; at the correct precision
test-gradient-correctness.R fails on icar_st and bym2_st until the
gradient carries the same per-margin constants. hmc_sampler.h's
build_and_factorize() keeps its 0.001, renamed diag_ridge: it is a
Tikhonov ridge lambda*I that makes a Kronecker of two intrinsic precisions
factorize, not the rank-1-per-margin sum-to-zero penalty.
The multiscale, TVC and SVC sum-to-zero terms are written once. The
double and autodiff twins of the multiscale kernels, tvc_sum_to_zero_penalty
and svc_sum_to_zero_penalty are collapsed into one body templated over the
scalar type. The missing multiscale pin survived a previous migration pass
precisely because it had to be found in six places; there is now one.
src/soft_sum_to_zero.h is gone in favour of the engine's
tulpa/soft_sum_to_zero.h.
Recovery tests cover the fields that were unconstrained.
test-recovery-constrained-fields.R gains multiscale (RW1 trend + cyclic
seasonal, and RW2 trend), TVC and spatiotemporal cases. Shape and class
assertions pass just as happily when a field's level is free, which is how
all three items reached a release. The TVC and spatiotemporal cases assert
recovery but not Rhat: neither model mixes to 1.05, which reproduces on the
released build for TVC (#23) and is a separate residual for the interaction
(#24). Both comments say so and name the issue to restore them under.
ModelData assigned the CSR adjacency without the
component partition the sum-to-zero augmentation iterates over, so the field's
constant direction carried no precision and its level random-walked:
mean(phi_spatial) reached +200.9, -342.8 and +423.1 at seeds 42, 43 and 44
where the legacy backend holds it within 0.02 of zero, and tau_spatial came
out at 16.89 against the legacy 11.90. The field's shape and the fixed effects
were unaffected, which is why it read as a healthy fit. The bridge now calls
ModelData::set_spatial_adjacency(), which derives the partition and the
component count from the adjacency (tulpa 0.0.108). The two B1d spec-vs-legacy
parity tests run again, having been skipped since the defect was filed.tratio() replaces ratiod() as the package's front door, joining the
tulpa* family of fitting verbs alongside tulpa::tulpa() and
tulpaObs::tobs(). The ratiod_* family constructors, the ratiod_fit
class, and every other ratiod_-prefixed name are unchanged.
tratio() carries statistical arguments in its signature and every perf,
numerical, and tuning knob in a single control = list(), matching the
convention the engine and tulpaObs already follow. Moved into control:
chains, iter, warmup, thin, cores, seed, verbose,
adapt_delta, max_treedepth, metric, riemannian, gradient_mode,
re_param, vi_variant, and the stochastic-gradient knobs (batch_size,
epsilon, alpha, schedule_*, use_schedule) that previously arrived
through .... An unrecognised knob is now an error rather than a silent
no-op, and a warmup that would leave no post-warmup draws is rejected.
refresh is removed. It was a documented argument that no backend ever
read, so passing it changed nothing.
Spatiotemporal models fitted under an autodiff gradient mode now sample the
density they are meant to (#14). The log posterior is written twice, once as
compute_log_post and once as the templated compute_log_post_impl that the
arena, forward, and tape modes differentiate, and the spatiotemporal
interaction was present only in the first. The interaction parameters
therefore carried no prior and contributed nothing to either linear
predictor, so their gradient was exactly zero and every other parameter was
differentiated against a linear predictor short one effect. The
interaction now enters the templated log posterior as well, covering the
Knorr-Held types I to IV, the separable and non-separable GP forms, the HSGP
basis form, and the non-centred Type IV parameterization. The interaction
math in src/hmc_spatiotemporal.h is templated in place rather than copied
into a second autodiff header, so both paths call one implementation.
gradient_mode = "handcoded" was correct throughout and is unchanged, to the
last bit.
The samplers no longer hold their per-thread scratch in the shape that
corrupts the heap on Windows (#16). Every workspace in hmc_sampler.cpp was
a static thread_local object owning heap buffers. Under the mingw toolchain
such an object gets emutls-backed storage and a destructor registered through
__cxa_thread_atexit, and the two are released by thread-exit hooks whose
relative order is unspecified, so the destructor can free through storage that
is already gone and the process dies at some later unrelated free. A worker
thread exits only when libgomp narrows a team, which is why the fault appeared
when a fit asked for fewer chains than the one before it. The 42 declarations
now go through RATIOD_TLS_WORKSPACE (src/tls_workspace.h), which keeps the
workspace behind a constant-initialized pointer: no initialization guard, no
destructor registration, nothing for a dying worker to free. No object in the
package registers a thread-exit destructor any more.
tratio() no longer fails against a tulpa version the requirement allows. It
calls tulpa::tulpa_check_control(), exported from tulpa 0.0.99, while
DESCRIPTION asked only for tulpa (>= 0.0.95).
diagnostics() replaces mcmc_diagnostics(), following the engine rename in
tulpa 0.0.95. It delegates to tulpa::diagnostics(), which selects chain or
approximation diagnostics from the fit's draws provenance rather than from the
function's name. mcmc_diagnostics() is deprecated and still returns the same
value.A fit no longer decides how many threads later fits get (#16). The Laplace
and Polya-Gamma backends take a per-fit cores, which they applied by moving
the process-wide OpenMP thread count and leaving it moved, and it defaults to
one thread. Every region in the session that read that value afterwards was
therefore sized by whichever of those fits ran last, in fits it knew nothing
about. Each backend now restores the previous value when it returns, on the
interrupt path as well as the normal one.
Successive fits in one session no longer corrupt the heap on Windows (#8).
The OpenMP team for chain-parallel work was sized per fit from cores, so a
fit asking for fewer cores than the one before it shrank the team; libgomp
then destroyed the surplus workers and ran their thread_local destructors
while later work was still in flight, and the damage surfaced as
STATUS_HEAP_CORRUPTION at an unrelated free(). Four chains followed by
two was enough to trigger it. The team now only ever grows, sized by the
widest core budget any fit has asked for, and cores bounds how many chains
are in flight instead, so it keeps its meaning. Growing a team was always
safe; only shrinking faulted. Its size never comes from
omp_get_max_threads(): the Laplace and Polya-Gamma backends move that value
through omp_set_num_threads() and leave it moved, so reading it let a
Laplace fit ahead of the first chain fit pin the team at one thread and
serialize the chains of every fit after it.
The Gibbs spatial backend runs its chains in parallel under OpenMP (#7). It
looped over chains in R, so a 4-chain fit cost close to four times a 1-chain
fit; on a 40-site binomial ICAR model that ratio drops from 3.7x to 0.94x.
cores now reaches this backend, bounded by the chain count and the machine.
Chain seeds are still derived in R, so a fit with a given seed returns the
same draws as before, and the parallel and serial paths agree bitwise.
The handcoded gradient no longer disagrees with its own log posterior when a
temporal effect is not shared between numerator and denominator. The
vectorized and fused kernels added the temporal effect to the denominator's
linear predictor and fed the denominator residual back into the temporal
gradient, both unconditionally, while compute_log_post() applied the effect
to the numerator alone. NUTS therefore sampled a density whose gradient
pointed elsewhere, silently, for every family with a real denominator
(poisson_gamma, negbin_negbin, negbin_gamma, gamma_gamma,
lognormal). Binomial models were unaffected. The gradient checks now cover
a continuous-denominator family with shared = FALSE, which is the only
configuration that can see this.
ModelData's scalar members are initialised at declaration. The constructor
set only unique_id, so roughly thirty members held indeterminate values
until assigned; the GP entry point never assigned re_parameterization and
read stack residue to choose between the centred and non-centred
parameterisation.
The scalar gradient fallback reduces its per-thread partial sums in thread-index order rather than in a critical section, whose summation order varies from run to run.
mcmc_diagnostics() now reports per-parameter Rhat and ESS (#4). It handed a
multi-parameter draws array to the single-variable posterior::rhat() /
ess_bulk() / ess_tail(), collapsing every parameter to one scalar that was
then recycled across the returned rows: a well-mixed four-chain fit reported
rhat = 2.124, ess_bulk = 1. Diagnostics now delegate to
tulpa::mcmc_diagnostics(), the engine's per-parameter implementation, which
removes the local re-derivation (compute_diagnostics_basic(),
compute_split_rhat(), compute_ess_basic()). summary(), plot_rhat(),
plot_ess(), diagnostic_summary() and check_diagnostics() all inherit the
correction, and a converged fit no longer prints a false non-convergence
warning.
print() on a ratiod_diagnostic_summary no longer errors when a fit has a
parameter with Rhat > 1.01 or ESS < 400. The worst-Rhat and worst-ESS tables
are column extracts of the diagnostics table and dispatched on
print.ratiod_diagnostics(), which rounded columns the extract does not
carry.
B1a PoC: plain binomial fits can route through tulpa's LikelihoodSpec
path (tulpa::get_nuts_fn()) behind a feature flag
(options(tulpaRatio.use_specs = TRUE) or TULPARATIO_USE_SPECS=1).
Off by default; legacy backend remains the production path until full
parity ships in B5. Posterior means agree with the legacy backend within
Monte Carlo noise.
B1b: the same feature flag now covers all 7 ratio families (binomial, poisson_gamma, negbin_gamma, negbin_negbin, gamma_gamma, lognormal, beta_binomial) for the simplest non-ZI / no-spatial / no-RE / single-chain configuration. Autodiff only — H-kernel port is deferred to B2.
B1c: zero-inflation, hurdle, and one-inflation variants now route through
the LikelihoodSpec path. Each family's templated likelihood dispatches on
data.zi_type and calls a shared mixture helper in
src/lik_specs/lik_helpers.h; per-family allowlist lives in
src/lik_dispatch.cpp and R/backend_hmc.R::SPEC_ZI_COMPAT. Posterior
parity with the legacy backend matches within Monte Carlo noise for every
binomial × {ZI, hurdle, OI, ZOIB} combination; for count families
(poisson_gamma, negbin_*) the legacy A_r path was already silently
ignoring ZI, so the spec path here gives the mathematically correct
posterior — recorded as a behaviour change, not a parity mismatch (see
dev_notes/B1c_zi_surface.md). Autodiff only; ZI variants for the
hand-coded H-kernel remain deferred.
The tulpa ABI test no longer pins a hardcoded version number. It compared
cpp_tulpa_abi_version() against a literal 32L, a second copy of a
constant tulpa already owns, so every tulpa ABI bump broke the test with a
stale-literal failure that said nothing about compatibility. Two unguarded
accessors (cpp_tulpa_compiled_abi_version(), the value baked in from the
LinkingTo headers, and cpp_tulpa_runtime_abi_version(), the value in the
loaded tulpa DLL) let the test assert the invariant that matters — that the
two agree — and report both numbers when they do not. tulpa remains the
single source of truth for the version.
Documentation is regenerated under roxygen2 8.0.0, matching the rest of the
packages that link against tulpa. RoxygenNote migrates to
Config/roxygen2/version. Two external links now resolve to a topic alias
rather than an Rd filename (posterior::as_draws_df(),
loo::stacking_weights()), which clears the "Non-topic package-anchored
link(s)" note from R CMD check; INFERENCE_TIERS drops the \docType{data}
and \format{} entries roxygen2 8.0.0 no longer emits for documented values;
and the package page lists authors alongside the maintainer.
Tiered inference mode system: Inference backends are now organized into three tiers based on epistemic guarantees:
Elliptical Slice Sampling (ESS): New Tier 1 backend for models with
Gaussian priors. Use mode = "ess".
Variational Inference (VI): New Tier 3 backend with mean-field, full-rank,
and low-rank variants. Use mode = "vi".
Stochastic Gradient MCMC: New Tier 3 backends for very large data:
mode = "sghmc"mode = "sgld"Thread-safe autodiff: A_t (tape autodiff) mode now uses RAII-based TapeScope for thread isolation. Parallel chains work correctly with all gradient modes.
Hand-coded gradients: Added O(N) analytical gradients for:
Spatial prediction: New predict() method for ratiod_fit objects
supports prediction at new spatial locations.
fwd::Dual) for thread-safe gradient computationSpatiotemporal interaction effects: New spatiotemporal() function specifies
space-time interactions using Knorr-Held (2000) Type I-IV models:
Non-separable spatiotemporal GP: New spatiotemporal_gp() function for continuous
space-time GP models with Gneiting or Cressie-Huang non-separable covariance.
Spatiotemporal effect extraction: New spatiotemporal_effects() extracts posterior
distributions in array, long, or summary format with visualization support.
Zero-inflated binomial: New ratiod_zibinomial() for proportions with excess zeros.
One-inflated binomial: New ratiod_oibinomial() for proportions with excess ones
(100% success/detection).
Zero-and-one inflated binomial (ZOIB): New ratiod_zoibinomial() for proportions
with excess at both boundaries.
Hurdle binomial: New ratiod_hurdle_binomial() for binomial data with hurdle
at zero.
vignette("spatial-temporal") with examples of space-time modeling.vignette("zero-inflation") covering ZI, OI, ZOIB, and hurdle binomial models.Latent factors for unmeasured confounders: New latent_factor() function allows
specifying latent factors that capture shared unmeasured variation between numerator
and denominator processes. Latent factors are observation-level random effects with
sum-to-zero constraints for identifiability. Use latent = latent_factor(n_factors = 1)
in the ratiod() call.
Latent factor extraction: New latent_factors() function extracts posterior
summaries or full draws for latent factor scores from fitted models.
vignette("random-effects") with section on latent factors.Initial release with:
ratiod_negbin_negbin(), ratiod_binomial(), ratiod_poisson_gamma()