Nothing
R/nested_laplace.R
In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models
Defines functions
.nl_checkpoint_args
.nl_progress_args
.prior_from_spatial_spec
.prior_from_temporal_spec
prior_from_spec
.nl_dispatch_multi
.nl_grid_layout
.nl_grid_crossing
.nl_check_grid_cap
.nl_max_grid_cells
.nl_with_internal_batch
.nl_internal_batch
.nl_dense_grid_warn
.nl_multi_grid_warn
.nl_block_axis_grid
.nl_block_spec_for_cpp
.nl_single_block_width
.nl_role_cols
.nl_block_roles
.is_multi_block_prior
.nl_nngp_field_args
.nl_block_projector
.nl_refuse_projector
.nl_block_is_projected
.nl_attach_grid_hessians
.nl_modal_joint_precision
`%||%`
.nl_attach_evidence
.nl_outer_log_evidence
.nl_normalise_weights_safe
tulpa_normalise_weights_safe
.nl_attach_outer_integration
.nl_screen_log_offset
.nl_dispatch
.mcar_default_logchol_grid
.nl_adj_args
.nl_check_block_fields
.nl_attach_pareto_k
.nl_prune_gate
.nl_check_screen_iters
.nl_check_prune_tol
tulpa_nested_laplace
Documented in
prior_from_spec
tulpa_nested_laplace
tulpa_normalise_weights_safe
#' Nested Laplace approximation for latent Gaussian models
#'
#' @description
#' Generic outer-grid nested Laplace driver. Builds a grid over the
#' hyperparameters of a single latent prior block (spatial or temporal),
#' runs an inner Laplace at each grid point with warm-starting, and
#' integrates over the grid to give proper hyperparameter marginals.
#'
#' Supported priors:
#' * Spatial (areal): `"icar"` (1D grid on tau), `"bym2"`
#' (2D on (sigma, rho)), `"car_proper"` (2D on (tau, rho); rho lives
#' in the eigenvalue interval (1/lambda_min, 1/lambda_max) of
#' `D^{-1} W`).
#' * Spatial (continuous): `"nngp"` (2D on (sigma2, phi_gp)),
#' `"hsgp"` (2D on (sigma2, lengthscale)).
#' * Temporal: `"rw1"`, `"rw2"` (1D grid on tau), `"ar1"` (2D on (tau, rho))
#' * SPDE continuous spatial: see \code{cpp_nested_laplace_spde()} (separate
#' entry, rebuilds Q via SPDE Q-builder).
#'
#' @param y Response vector.
#' @param n_trials Trial sizes (binomial). Pass `1L`-vector otherwise.
#' @param X Fixed-effects design matrix.
#' @param prior A list describing the latent prior block. Required field
#' `type` one of \{"icar", "bym2", "car_proper", "rw1", "rw2", "ar1"\}.
#' Type-specific fields:
#' * icar: `spatial_idx`, `n_spatial_units`, `adj_row_ptr`, `adj_col_idx`,
#' `n_neighbors` (CSR adjacency, 0-based); optional `tau_grid`.
#'
#' Any block accepts an optional `svc_weight`, one weight per observation,
#' which makes it a varying coefficient: observation `i` contributes
#' `svc_weight[i] * f_i` instead of `f_i`, where `f_i` is the block's field
#' at that row (`z[spatial_idx[i]]` for an areal block, `(A u)_i` for an
#' SPDE one). Absent, the field enters unweighted.
#' * bym2: same adjacency; `scale_factor`; optional `sigma_grid`, `rho_grid`.
#' * car_proper: same adjacency; optional `tau_grid`, `rho_grid`,
#' `rho_bounds = c(lower, upper)` (defaults to (0, 1)).
#' * rw1/rw2: `temporal_idx` (1-based), `n_times`; optional `tau_grid`,
#' `cyclic` (default FALSE).
#' * ar1: `temporal_idx`, `n_times`; optional `tau_grid`, `rho_grid`
#' (each `rho_grid` value strictly inside (-1, 1)).
#'
#' An `icar`, `bym2`, `car_proper` or `nngp` block may carry `projector`, an
#' `N x n_spatial_units` (`N x n_spatial` for `nngp`) matrix (dense or a
#' \pkg{Matrix} sparse matrix), in place of `spatial_idx`: observation `i`
#' then receives `(A z)_i`, so restricted spatial regression is
#' `A = S P_perp` with `S` the observation-to-unit incidence. The block is
#' fitted by the multi-block driver, alone or in a list of blocks. An
#' intrinsic field's level is placed from `A 1`: removed with no fold where
#' it is zero, folded into the intercept where it is the intercept column,
#' otherwise left to the augmentation. The `nngp` and `car_proper` priors are
#' proper and carry their level themselves.
#'
#' Grids of a two-axis type (`bym2`, `car_proper`, `ar1`) are PAIRED when
#' both are supplied: the i-th entries of `tau_grid` / `sigma_grid` and
#' `rho_grid` form the i-th outer-grid cell, so the two must have equal
#' length, and a tensor grid is written as `expand.grid()` columns. A single
#' supplied axis is kept as an axis (its distinct values) and crossed with
#' the default of the other.
#'
#' A default grid axis is a starting axis, not a hard ceiling: for `icar`
#' (`tau_grid`) and `bym2` (`sigma_grid`) a posterior mode that rails a
#' boundary node (`pareto_k_regime = "collapsed_edge"`) triggers one
#' mode-Hessian recenter-and-refit, reported through
#' `outer_grid_placement` / `outer_grid_recenter_declined` -- see
#' [tulpa_nested_laplace_joint()]'s return docs. An axis the caller pinned is
#' never moved; mark a grid your own code defaulted with [auto_grid()] to keep
#' the recenter live on it.
#' @param re_idx Optional 1-based RE group index per obs (defaults to no RE).
#' @param n_re_groups RE group count (default 0).
#' @param sigma_re RE standard deviation (default 1).
#' @param family `"binomial"`, `"poisson"`, `"neg_binomial_2"`, etc.
#' @template phi
#' @param offset Optional per-observation offset added to the linear predictor
#' (length `length(y)`), as an `offset()` term in a model formula. Every
#' single- and multi-block kernel carries it, and `fitted_eta` includes it.
#' @param control Optional list of perf/numerical tuning knobs (statistical
#' arguments stay top-level), following the `control` convention of
#' [tulpa()]. Recognised elements (defaults in parentheses):
#' * `max_iter` (`50L`), `tol` (`1e-6`) -- inner Newton iteration budget and
#' tolerance.
#' * `n_threads` (`1L`) -- inner-loop OpenMP threads.
#' * `x_init` (`NULL`) -- warm-start for the first grid point's inner solve.
#' * `keep_grid_hessians` (`TRUE`) -- when `TRUE`, retain per-grid-point
#' fixed-effects marginal Hessian \eqn{H_\beta} and mode \eqn{\hat{\beta}}
#' on the return list as `$grid_hessians` (list of dense \eqn{p\times p}
#' matrices) and `$grid_modes` (list of length-\eqn{p} vectors). Used
#' downstream by simplified-Laplace (SLA) callers to assemble skew-aware
#' marginals -- see the cumulant pooling in [rubins_pool()].
#' * `diagnose_k` (`TRUE`), `k_samples` (`500L`) -- compute the outer
#' Pareto-\eqn{\hat{k}} accuracy diagnostic (`$pareto_k`) by importance
#' sampling the hyperparameter posterior against the Gaussian proposal
#' fitted to the grid, drawing `k_samples` extra inner-marginal evaluations.
#' `k_samples` is a precision knob: the GPD tail size is held at the
#' fraction the default budget implies, so a larger budget sharpens the
#' same k-hat instead of moving it to a deeper quantile of the weight
#' distribution (gcol33/tulpa#631). `k_tail_points` overrides that tail
#' size directly (capped at 20% of the draws).
#' Computed for a single-block, single positive-scale-axis grid; left `NA`
#' (with the grid's quadrature ESS as the fallback diagnostic) for
#' multi-block, multi-axis, or bounded-parameter grids. See [tulpa_psis()].
#' * `diagnose_skew` (`TRUE`), `skew_idx` (`NULL`) -- compute the inner-Laplace
#' skewness diagnostic (`$inner_skew`, gamma_3, Rue Martino & Chopin 2009
#' Sec 3.2.3) at the fitted MAP grid cell: one extra Newton solve, scoring
#' the `p` fixed-effects latent indices by default (pass `skew_idx`, 1-based
#' latent indices, to score additional ones, e.g. specific spatial units --
#' the full latent field is not scored by default since it costs one linear
#' solve per index). This is the complementary layer to `diagnose_k`: the
#' outer diagnostic scores the hyperparameter-grid integration around a
#' FIXED inner Laplace, this scores whether that inner Gaussian
#' approximation is itself a good fit to the latent-field conditional
#' posterior. See [diagnostics()] for the combined whole-fit verdict.
#' * `within_cell` (`"log_quadratic"`) -- the WITHIN-CELL construction the
#' reported per-axis hyperparameter intervals are read with. The outer
#' grid's weights say how much mass each cell holds; they do not say how it
#' is spread inside the cell, and a quantile needs both.
#' `"log_quadratic"` joins each row's node log densities by overlapping
#' quadratics, exact for a Gaussian at any spacing, and continues the end
#' quadratic past the outer node instead of stopping at the box edge; it
#' runs on the rows that resolve their own conditional (`h / sd <= 2`) and
#' reads every other row as boxes, declining to `"box_uniform"` where no
#' row qualifies. `"box_uniform"` puts the cumulative FULL mass at each
#' cell EDGE and interpolates between edges; `"chord"` puts the cumulative
#' MID-mass at each cell coordinate and interpolates between coordinates.
#' The default moved from `"box_uniform"` at gcol33/tulpa#932: a placed
#' axis sits at 1.25 posterior SDs per cell, where the box read is narrow
#' by construction, and against dense references of a two-arm ICAR
#' fixture the log-quadratic read's summed error is 0.038 against the box
#' read's 0.110 on the same grids. `outer_grid_h_over_sd` is how wide a
#' cell is on each axis, and `theta_within_cell` is what each axis was
#' actually read with. Only a `"density"` support admits the box and
#' log-quadratic reads -- a CCD design, a locally refined grid and a
#' posterior sample are not cell partitions that tile -- and an axis they
#' decline on reports the read that ran with a reason rather than erroring.
#' Nothing else moves: point estimates, moments, draws and weights are
#' untouched.
#' * `skew_correct` (`TRUE`) -- consume the inner-Laplace expansion instead of
#' only grading it: report
#' Cornish-Fisher marginal quantiles at each coefficient's own `gamma_3`,
#' about the centre `gamma_1 + gamma_3 / 2` that Rue, Martino & Chopin
#' (2009) eq. (22) implies, from `summary()` / `confint()` wherever the
#' combined inner band says the leading-order expansion is in its regime, and
#' the Gaussian quantiles everywhere else. It is post-processing on the
#' reported quantiles: draws, modes and weights are untouched, so a fit run
#' with it off is bit for bit the fit it was before. A coefficient whose
#' location term could not be formed declines rather than reading it as zero.
#' The band that bounded the relocation itself (`centre_unreliable`) is
#' off (`Inf`): scored over seven fixtures with an exact
#' reference, every finite cutoff declines the coefficients the correction
#' helps most, because a large centre carrying a small `gamma_3` is uniformly
#' weak correlation rather than an expansion out of its regime.
#' `$skew_correction` records the per-coefficient `gamma_3`, `gamma_1` and
#' the centre they form, the band, the inner importance k-hat, the combined
#' band, the eligibility and the reason behind it; the `skew_applied`
#' attribute on `summary()` / `confint()` records what was actually used at
#' the requested level. RMC fit a skew normal here instead; the series
#' correction is the same-order alternative, and unlike a skew normal its
#' skewness does not saturate inside the band it is applied on.
#'
#' MEASURED. Against exact quadrature
#' quantiles of rare-event binomial-logit posteriors it cuts total absolute
#' endpoint error 69.2%, improving both endpoints in every case. Scored over
#' the WHOLE marginal -- paired CRPS against the exact posterior in a
#' 400-replicate prior-predictive experiment -- it reads -0.01643 against
#' the uncorrected Laplace at t = -1.89, essentially all of the -0.01662 the
#' exact posterior itself achieves, and its PIT re-enters the simultaneous
#' SBC band. Applied about the Laplace mode instead of about
#' `gamma_1 + gamma_3 / 2` the same reshaping scored +0.00775 at t = +3.54,
#' a net loss; the location term is what supplies that centre.
#'
#' IT IS ON BY DEFAULT, so `summary()` / `confint()` on a
#' nested-Laplace fit report the corrected quantiles wherever the combined
#' inner band admits the coefficient; `skew_correct = FALSE` restores the
#' uncorrected report exactly. Scored against the mixture read a
#' correction-off fit gives, the flip is t = -1.895 on the rare-event
#' intercept and -3.765 / -3.201 on the small-group Bernoulli design. Across
#' twelve model classes read off one solve per seed, pooled 95% coverage
#' moves 0.9510 -> 0.9542 at a standard error of 0.0070, with every class
#' inside the acceptance the shipped gates use. A fit the correction cannot
#' help -- a coupled one, whose location term is unreachable -- reports what
#' it reported before, to the bit.
#' * `subspace_debias` (`FALSE`) -- correct only the latent directions the
#' inner-layer diagnostics flagged, by exact Metropolis, and leave the rest
#' at their Gaussian conditional. `TRUE` takes every default; a list
#' overrides `band`
#' (the inner-reliability floor a coordinate is selected at, default
#' `"ok"`), `idx` (pin the corrected set explicitly, skipping the selector),
#' `closure` / `closure_max` (grow the set by strongly coupled
#' precision-graph neighbours: a coordinate strongly coupled to a member of
#' the corrected set and left out of it is carried linearly, which is the
#' error the correction removes, so requesting the closure also retains the
#' modal cell's joint precision for the selector to read), the sampler
#' budget `n_iter` / `warmup` / `thin`, and
#' `n_draws`. The selector reads the per-index bands `diagnose_skew` already
#' attached, so it costs no extra solve; the correction itself re-runs the
#' settled grid once with the sampler on, and the fit then reports `$draws`
#' -- the per-cell Metropolis sample for the selected coordinates, the rest
#' from the Gaussian conditional given them -- instead of the
#' Gaussian-mixture moments. An EMPTY selection leaves the fit bit-for-bit
#' identical to the plain path. `$subspace_debias` records what was selected,
#' the bands it was read from, and the per-cell acceptance rate. Requesting
#' it turns `keep_grid_hessians` on, since the recombination reads exactly
#' those per-cell pieces.
#' * `cila` (`FALSE`) -- corrected integrated Laplace, the second inner-layer
#' debias (after Lai, Margossian and Sheldon,
#' arXiv:2605.20345). Where
#' `subspace_debias` selects coordinates and runs exact Metropolis on them,
#' this selects nothing: at every outer cell it draws `n_points` points from
#' the whole inner Gaussian, weights each by the exact joint density it came
#' from, and reports the weighted particles. `TRUE` takes the defaults; a
#' list overrides `n_points` (`1024L`), `variant` (`"qmc"`, a Sobol net;
#' `"is"` for iid draws, `"rqmc"` for the net under `n_shift` random
#' shifts), `n_shift` (`8L`), `n_draws` and `seed`. Below 512 points a
#' cell's particle set is too coarse to be a marginal at all and the request
#' is refused. The corrected per-cell masses become the fit's own
#' `weights` / `log_marginal` and `weights_source` reports `"cila"`; the
#' pre-correction pair is kept as `$cila$laplace`. A cell whose inner solve
#' factorized sparsely draws through the CHOLMOD factor's own triangular and
#' permutation solves; an LDL' factor has no square root to draw with and is
#' declined with `"sparse_factor_not_ll"`.
#' * `auto_recenter` (`TRUE`) -- outer-grid placement policy. `TRUE`
#' re-centres the movable default
#' axes on the posterior mode and refits when the grid either RAILS (an
#' axis's own marginal is maximal at one of its own endpoints) or does not
#' RESOLVE its posterior (an axis's node spacing exceeds 2 posterior SDs in
#' its own coordinate). Both tests read the weights the fit already stored,
#' so a grid that brackets and resolves its mode costs nothing beyond them;
#' when the pass does fire it is a second full grid solve plus a
#' finite-difference mode/Hessian stencil. A recentred axis is
#' `mode +/- 2.5 sd` over 5 nodes, a cell width of 1.25 posterior SDs by
#' construction, against a census median of 3.9 on the fixed spans.
#'
#' Measured over 200 fixed-truth seeds on each of six configurations (icar
#' chain / icar lattice / rw1 / bym2 / iid / nngp), the default moves mean
#' |coverage - nominal| from 0.043 to 0.030 at the 95% level, 0.171 to
#' 0.084 at 80% and 0.243 to 0.129 at 50% against the rail-only policy, at
#' 0.63 times the 95% interval width and 0.76 times the median bias, for
#' 1.71 times the wall clock.
#'
#' Three other values. `"rail"` is the rail test alone, which is what
#' `TRUE` meant before the sizing measurement settled the default. `FALSE`
#' integrates over the grid exactly as given, whatever it is, and records
#' `outer_grid_recenter_declined = "auto_recenter_disabled"`. `"always"`
#' re-centres every movable default axis whatever the fit did, at 2.04
#' times the wall clock; it agrees with the default seed for seed on five
#' of the six measured configurations and differs on the one whose default
#' axes already resolve their posterior, where it takes 50% coverage to
#' 0.135 against a nominal 0.5.
#'
#' Which families carry movable axes is `.NL_REGISTRY_AXIS_FIELD`: icar,
#' rw1, rw2, iid, bym2, nngp, hsgp and spde. car_proper, ar1 and hsgp_mo
#' each carry a correlation axis with no guessable coordinate, and mcar /
#' miid / tgmrf hold their axes in one matrix field; a fit of any of them
#' records which through `outer_grid_recenter_declined` rather than passing
#' in silence. The per-axis policy names are the standalone registry path
#' only -- [tulpa_nested_laplace_joint()] and [fit_st_nested()] refuse them
#' rather than accept them and ignore them.
#' * `max_grid_cells` (`NULL`, no ceiling) -- a ceiling on a multi-block
#' outer grid's cell count, one inner Newton solve per cell; a grid past it
#' is refused before any solve, naming the axes that produced the count.
#' Unset, a grid of any size is solved, and one past 2048 cells is announced
#' with a warning first.
#' * `checkpoint` (`list(path =, resume =)`) -- grid-cell checkpoint /
#' resume. Each solved outer cell is appended to `path`, keyed by its
#' hyperparameter coordinate; `resume = TRUE` loads the finished cells and
#' solves only the rest, and a fingerprint mismatch (different data,
#' settings or grid) errors rather than resuming onto a stale result. See
#' `?tulpa_nested_laplace` "Checkpoint / resume".
#' * `progress`, `progress.every`, `progress.file`, `progress.throttle` --
#' the outer-grid progress reporter: whether to print, how often (in
#' cells), where to, and the minimum seconds between lines.
#' * `prune` (`TRUE`), `prune_tol` (`1e-3`), `screen_iters`
#' (`.NL_SCREEN$iters`, 2 -- the depth gcol33/tulpa#640 measured the screen
#' down to) -- cheap-pass screening of the outer grid. The driver first
#' sweeps the lattice running a `screen_iters`-step inner Newton per cell,
#' each warm-started from the previous screened cell's quasi-mode, computes
#' a screening Laplace log-marginal, softmax-normalises it, and skips the
#' full inner Newton on every cell whose screened weight is below
#' `prune_tol`. The neighbour warm-start keeps each cheap mode near its
#' cell's own mode, so the cheap ranking tracks the full-solve ranking even
#' where the latent mode moves across the grid. Pruned cells get
#' `log_marginal = -Inf`, `n_iter = 0` and inherit the pilot mode; the pilot
#' cell is never pruned, and at least `prune_min_keep` cells (5) are solved
#' in full whatever the tolerance says -- the highest-ranked dropped cells
#' are restored up to that floor, because the outer grid placement pass
#' reads a finite-difference curvature off the cells that were SOLVED and a
#' kept set of one carries none. After the full pass the posterior mass the
#' dropped cells may carry is bounded by lifting each by the worst screening
#' error measured on a kept cell (`prune_dropped_mass_bound`); past 1% the
#' driver solves the dropped cells that contribute most to it until it is
#' under 1% (`prune_n_repaired`), and a safety gate falls back to the full
#' grid (with a warning) if the bound still fails. The bound assumes a
#' dropped cell is screened no worse than the worst kept one. On by default
#' (gcol33/tulpa#925); `FALSE` solves every cell.
#' `prune_tol` must be a finite numeric in `[0, 1)` and `screen_iters` a
#' single integer `>= 1`; both are validated whatever `prune` says.
#' `prune` / `prune_tol` reach both the single-block and the multi-block
#' dispatch; the screening DEPTH is declared by the single-block kernels
#' only, so a multi-block prior refuses a pinned `screen_iters` rather than
#' accepting it and ignoring it.
#' * `prune_log_gap` (unset) -- the screening cut stated in nats instead of as
#' a normalised weight: keep every cell within this many nats of the best
#' screened cell. `prune_tol` is a weight, so what it cuts at is the gap
#' `-(log(prune_tol) + log(Z))`, and the whole range a caller would type
#' (`1e-3` to `1e-12`) is 7 to 28 nats of it -- on an outer surface whose
#' log-marginal spans thousands of nats every setting returns the same kept
#' set. Stating the gap directly keeps the knob's resolution wherever the
#' surface is steep. It reaches the kernel as `prune_tol = exp(-gap)`, so
#' the realised cut is the requested gap less `log(Z)`, at most `log(n_grid)`
#' nats narrower; the fit reports the realised value. Setting both this and
#' `prune_tol` is an error -- they state one cut in two units.
#' * `fitted_var` (`TRUE`) -- fill the per-row predictive variance
#' `fitted_eta_var`, the companion of `fitted_eta` a caller marginalises the
#' per-row linear predictor with. It costs one back-solve per distinct
#' loading vector per cell, which on a design with few repeated rows is the
#' dominant cost of a cell, so a caller reading only `fitted_eta` or the
#' coefficient summaries sets `FALSE` and the fit then carries no
#' `fitted_eta_var`. Every single-block field reports it, the NNGP, HSGP,
#' SPDE and spatiotemporal fields included, each read off its cell's own
#' precision at the mode; a multi-block prior refuses `FALSE` rather than
#' accepting it and ignoring it.
#' * `var_of_means_consistency` (`TRUE`) -- the consistency pass of
#' [tulpa_nested_laplace_joint()]: a field SD (`tau` of an `icar`, `rw1` or
#' `rw2` block, `sigma` of an `iid`, `bym2` or `spde` block) whose marginal
#' has collapsed onto too few nodes to carry a spread has levels added in
#' every row of the other axes that holds the posterior until it resolves, and the new cells are
#' solved and carried like the rest (`var_of_means_consistency_info`
#' records the pass). Other axes keep the nodes they were laid on.
#' * `adaptive_grid` (`FALSE`), `adaptive_grid_edge_thresh` (`0.02`),
#' `adaptive_grid_max_passes` (`1L`) -- the opt-in boundary / interior pass,
#' on the same axes.
#' * `axis_refine` (`NULL`) -- per-axis override of how far those axes may be
#' moved, named by the grid's axes (`"tau"`; `"b1.tau"` on a multi-block
#' prior): `"none"`, `"densify"` (inside the declared span) or `"extend"`.
#' Nodes the caller wrote down densify by default and an engine-placed axis
#' extends, as in [tulpa_nested_laplace_joint()].
#'
#' @return A list with:
#' * `theta_grid`: matrix or vector of grid hyperparameter values.
#' * `log_marginal`: log p(y, mode | theta_k) at each grid point.
#' * `weights`: integration weights normalising to sum 1.
#' * `theta_mean`, `theta_sd`: posterior moments per hyperparameter.
#' * `n_iter`: inner Newton iterations per grid point.
#' * `modes`: matrix `[n_grid x n_x]` of inner modes, when stored.
#' * `pareto_k`, `pareto_k_is_ess`: outer Pareto-\eqn{\hat{k}} and its
#' importance-sampling ESS (`NA` when not computed for the grid; see
#' `control$diagnose_k`).
#' * `pareto_k_proposal_source`, `pareto_k_first_pass`: which proposal family
#' the reported k-hat came from (the outer k-hat is scored against several
#' candidates and the best is kept), and the k-hat of the FIRST pass -- the
#' proposal exactly as this backend placed it, before any candidate refined
#' or replaced it. A large gap between the two says the nodes are badly
#' scaled around the hyperparameter posterior even though the verdict is
#' fine; no gap says the placement was already right.
#' * `inner_skew`, `inner_skew_idx`, `inner_skew_dropped`: the inner-Laplace
#' skewness diagnostic (gamma_3) at each scored latent index and its
#' 1-based index, plus a count of (index, observation) contributions
#' dropped for a non-finite third derivative (see `control$diagnose_skew`).
#' * `timing`: named numeric of wall-clock seconds (`total`, `setup`,
#' `grid`, `postproc`, `diagnostics`); the `grid` phase is the inner
#' Laplace pass that scales with grid size. Surfaced one-line in `print`.
#' * `prune_cheap_log_marginal`, `prune_mask`, `prune_n_pruned`, `prune_tol`:
#' present when the grid was screened and the safety gate did not
#' trip -- the cheap-screen log-marginal per cell (the truncated solve's
#' value plus half its Newton decrement, the second-order estimate of the
#' converged value the cells are ranked by; the decrement itself is
#' `prune_screen_decrement`), the logical mask of
#' pruned cells, the pruned-cell count, and the threshold actually applied.
#' * `prune_log_gap_cut`, `prune_cheap_lm_spread`, `prune_min_keep`,
#' `prune_n_floor_restored`: the same screen read on the scale it operates
#' on -- how many nats below the best screened cell the tolerance cut at,
#' how many nats the screened surface spans, the kept-cell floor applied,
#' and how many cells that floor put back. A cut that is a sliver of the
#' spread is a tolerance with no resolution on this grid.
#' * `prune_dropped_mass_bound`, `prune_n_repaired`, `prune_gate_mass`: the
#' posterior mass the dropped cells may carry, with each lifted by the worst
#' screening error measured on a kept cell; how many dropped cells the
#' driver solved in full to bring that under the gate; and the gate.
#' * `prune_fallback_triggered`, `prune_fallback_reason`: present only when
#' the gate did trip; the rest of the fit is then the full-grid result and
#' the reason string records which check fired.
#' * `prior`: echoed input.
#'
#' @param spec Optional `tulpa_temporal` or `tulpa_spatial` spec object
#' (output of [temporal_rw1()], [temporal_rw2()], [temporal_ar1()],
#' [spatial_car()], [spatial_bym2()], etc.). When supplied alongside
#' `data`, the `prior` list is built automatically via
#' [prior_from_spec()] -- pass either `prior` or `spec`, not both.
#' @param data Data frame used to validate `spec` and resolve
#' time/group/site indices. Required when `spec` is supplied.
#' @param likelihood Optional model-supplied likelihood, replacing the built-in
#' `family`. Pass an external pointer to a `tulpa::NestedLikelihood` (built in
#' a model package's own C++ from a `LikelihoodSpec`); the inner Laplace solve
#' then reads the per-observation score, Fisher weight, and log-likelihood from
#' that spec instead of `family`, so `family`/`phi` are ignored. Used by model
#' packages to fit a custom response without adding a family to tulpa -- for
#' example tulpaObs threads its marginalized single-season occupancy
#' likelihood (a scaled Bernoulli, with the latent occupancy state integrated
#' out) through this. Multi-block `prior` only. Default `NULL` (use `family`).
#' @template hyperprior
#'
#' @references
#' Rue, Martino & Chopin (2009). Approximate Bayesian inference for latent
#' Gaussian models by using integrated nested Laplace approximations.
#' \emph{JRSS-B} 71(2):319-392.
#' @examples
#' \donttest{
#' set.seed(1)
#' S <- 30L # spatial units arranged in a chain
#' nb <- lapply(seq_len(S), function(s) setdiff(c(s - 1L, s + 1L), c(0L, S + 1L)))
#' nn <- lengths(nb)
#' field <- as.numeric(scale(cumsum(rnorm(S, 0, 0.4)))) # smooth spatial field
#' idx <- rep(seq_len(S), each = 6L); n <- length(idx); x <- rnorm(n)
#' y <- rbinom(n, 1L, plogis(-0.2 + 0.6 * x + field[idx]))
#' prior <- list(type = "icar", n_spatial_units = S, spatial_idx = idx,
#' adj_row_ptr = c(0L, cumsum(nn)), adj_col_idx = unlist(nb) - 1L,
#' n_neighbors = nn, tau_grid = c(0.5, 1, 2, 4, 8))
#' fit <- tulpa_nested_laplace(y, rep(1L, n), cbind(1, x), prior = prior,
#' family = "binomial")
#' fit$theta_mean # marginalized ICAR precision
#' }
#' @export
tulpa_nested_laplace <- function(y, n_trials, X, prior = NULL,
spec = NULL, data = NULL,
re_idx = NULL, n_re_groups = 0L, sigma_re = 1.0,
family = "binomial", phi = 1.0, offset = NULL,
likelihood = NULL,
hyperprior = c("proper", "flat"),
control = list()) {
tulpa_check_control(control, .CONTROL_KEYS$nested_laplace, "tulpa_nested_laplace")
# A typo is named here, against the family registry, rather than by the C++
# family lookup, whose message names a header (gcol33/tulpa#886). A consumer
# `likelihood` supplies its own density and bypasses the builtin family.
if (is.null(likelihood)) .family_or_stop(.canonical_family(family))
hyperprior <- .hp_choice(match.arg(hyperprior))
tm <- .tulpa_timer()
# Perf/numerical knobs live in `control = list()` (matching tulpa()); the
# top-level signature carries only statistical arguments.
max_iter <- control$max_iter %||% 50L
tol <- control$tol %||% 1e-6
n_threads <- control$n_threads %||% 1L
x_init <- control$x_init
# Cheap-pass outer-grid screening. The kernels read a positive tolerance as
# "screen the lattice and prune", so the user-facing `prune` toggle gates it:
# FALSE always solves every cell, TRUE screens at `prune_tol`. Both the
# tolerance and the screening depth are validated whatever `prune` says, so a
# misspecified value is reported where the user set it rather than at the next
# fit that happens to switch pruning on.
prune <- isTRUE(control$prune %||% .nl_screen("prune")[["registry"]])
# `prune_log_gap` states the same cut in nats -- the units the screened
# surface is measured in -- and replaces the tolerance when set.
prune_tol <- .nl_prune_tol_from_control(
control,
.nl_check_prune_tol(control$prune_tol %||% .nl_screen("prune_tol")))
prune_tol_eff <- if (prune) prune_tol else 0
screen_iters <- .nl_check_screen_iters(
control$screen_iters %||% .nl_screen("iters"))
fitted_var <- isTRUE(control$fitted_var %||% TRUE)
# Subspace debias. The correction reports per-cell DRAWS
# instead of the Gaussian-mixture moments, and those are built from the same
# per-cell fixed-effect mode + precision `.nested_fixed_moments()` reads, so
# requesting it turns that retention on: without it there is nothing to
# recombine the sampled coordinates with.
sd_cfg <- .subspace_debias_config(control$subspace_debias)
# Corrected integrated Laplace, the second inner debias.
# Same shape as the subspace one: a second pass over the settled grid whose
# per-cell particles are pooled into the reported draws.
cila_cfg <- .cila_config(control$cila)
keep_grid_hessians <- isTRUE(control$keep_grid_hessians %||% TRUE) || !is.null(sd_cfg)
# Retaining the modal cell's joint precision is what makes the debias closure
# available on a grid fit; nothing else reads it, so it rides the request.
keep_joint_prec <- !is.null(sd_cfg) && !identical(sd_cfg$closure, FALSE)
diagnose_k <- isTRUE(control$diagnose_k %||% TRUE)
k_samples <- as.integer(control$k_samples %||% .nl_diag("k_samples"))
k_tail_points <- control$k_tail_points
diagnose_skew <- isTRUE(control$diagnose_skew %||% TRUE)
skew_idx <- control$skew_idx
skew_correct <- isTRUE(control$skew_correct %||% .nl_diag("skew_correct"))
# Which WITHIN-CELL construction the reported per-axis hyperparameter
# intervals are read off. A numerical reconstruction of the
# integration design the fit already ran, so it sits beside `integration` /
# `local_ccd` / `skew_correct` in `control` rather than in the signature,
# which carries only statistical arguments. The default is
# `.nl_diag("within_cell")`, and `"chord"` restores the pre-0.0.188 report
# byte-identically.
within_cell <- .nl_within_cell_mode(control$within_cell)
# Outer-grid refinement: the passes the joint drivers run
# (`.joint_refine_outer_grid()`), over the field SD axes this door can write a
# cell onto. The consistency pass is on by default, the boundary / interior
# pass opt in, as there.
var_of_means_consistency <- isTRUE(control$var_of_means_consistency %||% TRUE)
adaptive_grid <- isTRUE(control$adaptive_grid %||% FALSE)
adaptive_grid_edge_thresh <- control$adaptive_grid_edge_thresh %||% 0.02
adaptive_grid_max_passes <- control$adaptive_grid_max_passes %||% 1L
# Grid-cell checkpoint/resume. `control$checkpoint =
# list(path =, resume =)` makes every grid cell append to `path`; a resume
# loads the finished cells and solves only the rest. A fresh (resume = FALSE)
# run removes any stale file once here, before the first kernel call, so the
# several within-fit kernel calls (initial grid, adaptive refinement) all
# append rather than truncate each other. The path threads into `cargs` and
# reaches both the single-block (.nl_dispatch) and multi-block
# (.nl_dispatch_multi) kernels by name; the k-hat diagnostic re-evaluations
# are run with it stripped so they do not pollute the file.
.ckpt <- .nl_checkpoint_args(control)
if (nzchar(.ckpt$path) && !isTRUE(.ckpt$resume) && file.exists(.ckpt$path)) {
file.remove(.ckpt$path)
}
# The caller's outer-grid ceiling, published for the duration of this fit so
# the initial dispatch and every re-dispatch (the k-hat re-evaluations, the
# subspace-debias refit) read the same value.
.op_grid_cap <- options(tulpa.nl_max_grid_cells = .nl_max_grid_cells(control))
on.exit(options(.op_grid_cap), add = TRUE)
if (!is.null(spec)) {
if (!is.null(prior)) {
stop("Pass either `spec` or `prior`, not both.", call. = FALSE)
}
if (is.null(data)) {
stop("`data` is required when `spec` is supplied.", call. = FALSE)
}
prior <- prior_from_spec(spec, data)
}
# Allow passing a single tgmrf S3 object directly: wrap it in a length-1
# block list so the multi-block dispatch picks it up. This keeps the
# ergonomic API close to the formula sketch (one block,
# one call) while still routing through the validated multi-block path.
if (inherits(prior, "tulpa_latent_block")) {
prior <- list(prior)
}
if (!is.list(prior)) {
stop("`prior` must be a list (single block) or list-of-lists (multi-block).",
call. = FALSE)
}
# This door's outer grid is spanned by the latent blocks' hyperparameters
# alone, so a prior with no block leaves nothing to integrate. The joint door
# integrates per-arm dispersion axes as well and does take an empty prior.
if (!length(prior)) {
stop("`prior` holds no latent block, so tulpa_nested_laplace() has no ",
"hyperparameter to integrate. Fit a fixed-effects model with ",
"tulpa(mode = \"laplace\"), or integrate a dispersion axis with ",
"tulpa_nested_laplace_joint(prior = list(), phi_grid = ).",
call. = FALSE)
}
# A projected areal field is read by the multi-block driver alone, so a
# single projected block routes there as a length-1 block list.
if (.nl_block_is_projected(prior)) prior <- list(prior)
# Outer-axis provenance: record which grid axes the caller
# declared as defaults with `auto_grid()` -- the registry rescue below reads
# it -- and strip the markers so every downstream consumer sees plain grids.
.prov <- .nl_grid_provenance(prior)
prior <- .prov$prior
# And refuse an axis this path cannot read: a grid field
# some OTHER driver spells the same family with is neither integrated nor
# reported here, so a pinned one is an error and a defaulted one is recorded.
.op_axis <- .nl_publish_axis_dropped(
.nl_check_axis_fields(prior, "registry", auto = .prov$auto))
on.exit(options(.op_axis), add = TRUE)
# A model-supplied `likelihood` is only honoured by the spec-driven multi-block
# driver (run_multi_block_nested_laplace); the single-kernel path (.nl_dispatch)
# has no spec hook. Route a single LatentBlock-type prior through the
# multi-block driver by wrapping it as a length-1 block list -- the per-cell
# math is identical (every single-block kernel already calls the same driver).
if (!is.null(likelihood) && !is.null(prior$type)) {
prior <- list(prior)
}
vd <- .validate_glm_design(y, X, n_trials, "tulpa_nested_laplace")
N <- vd$N
n_trials <- vd$n_trials
if (is.null(re_idx)) re_idx <- rep(0L, N)
re_idx <- .validate_re_idx(re_idx, n_re_groups, N, "tulpa_nested_laplace")
if (!is.null(offset)) {
offset <- as.numeric(offset)
if (length(offset) != length(y) || any(!is.finite(offset))) {
stop("`offset` must be a finite numeric vector of length length(y) (",
length(y), ").", call. = FALSE)
}
}
cargs <- list(
y = as.numeric(y),
n = as.integer(n_trials),
offset_nullable = offset,
X = X,
re_idx = as.numeric(re_idx),
n_re_groups = as.integer(n_re_groups),
sigma_re = as.numeric(sigma_re),
family = family,
phi = .phi_to_kernel(family, as.numeric(phi)),
max_iter = as.integer(max_iter),
tol = as.numeric(tol),
n_threads = as.integer(n_threads),
x_init_nullable = x_init,
store_Q = isTRUE(keep_grid_hessians),
checkpoint_path = .ckpt$path,
prune_tol = prune_tol_eff,
screen_iters = screen_iters,
compute_fitted_var = fitted_var,
# Read by the dispatchers, which fold the outer hyperprior in R; never
# handed to a kernel.
hyperprior = hyperprior
)
# cargs without the checkpoint, for the k-hat diagnostic re-evaluations, and
# without the prune: a screened cell carries `log_marginal = -Inf`, which is
# not the inner marginal the importance ratio and the placement stencil read.
# They evaluate the grid they are handed, so every cell of it has to be solved.
cargs_no_ckpt <- utils::modifyList(
cargs, list(checkpoint_path = "", prune_tol = 0))
p_fixed <- ncol(X)
if (.is_multi_block_prior(prior)) {
# The multi-block entry takes the prune tolerance but declares neither the
# screening depth nor the per-row variance switch, so a pinned one is
# refused here rather than accepted and ignored.
if (!is.null(control$screen_iters)) {
stop("`control$screen_iters` applies to a single-block `prior`; ",
"this fit takes the multi-block dispatch.", call. = FALSE)
}
if (!is.null(control$fitted_var) && !isTRUE(control$fitted_var)) {
stop("`control$fitted_var` applies to a single-block `prior`; ",
"this fit takes the multi-block dispatch.", call. = FALSE)
}
tm$mark("setup")
# Every full-grid solve of this fit -- the first one and the placement
# pass's refit -- is screened at the same tolerance and gated the same way,
# so no reported grid is a pruned one whose screen ranking was unreliable.
solve_grid_multi <- function(prior_i) {
out <- .nl_prune_gate(
.nl_dispatch_multi(cargs, prior_i, likelihood = likelihood,
progress = .nl_progress_args(control),
within_cell = within_cell),
prune_tol_eff,
function() .nl_dispatch_multi(
utils::modifyList(cargs, list(prune_tol = 0)), prior_i,
likelihood = likelihood, progress = .nl_progress_args(control),
within_cell = within_cell))
refined <- .nl_refine_registry(
out, function(theta_mat, screen = NULL) .nl_dispatch_multi(
.nl_refine_cargs(cargs_no_ckpt, screen, prune_tol_eff), prior_i,
likelihood = likelihood, theta_grid_override = theta_mat),
prior_i, TRUE, .prov$auto,
consistency = var_of_means_consistency,
adaptive_grid = adaptive_grid,
edge_thresh = adaptive_grid_edge_thresh,
max_passes = adaptive_grid_max_passes,
user = control$axis_refine)
if (is.null(refined$refining_axis)) return(refined)
# The dispatcher integrated the base grid; the appended cells are measured
# and read with the rest.
refined <- .nl_attach_outer_integration(refined, "multi-block outer grid")
.nl_posterior_moments_multi(refined, refined$blocks, refined$axis_offsets,
refined$theta_grid, within = within_cell)
}
res <- solve_grid_multi(prior)
tm$mark("grid")
if (isTRUE(keep_grid_hessians)) {
res <- .nl_attach_grid_hessians(res, p_fixed, keep_joint_prec)
}
tm$mark("postproc")
res <- .nl_attach_pareto_k(res, prior, cargs_no_ckpt, "multi", NULL,
likelihood, k_samples, compute = diagnose_k,
tail_points = k_tail_points)
# The same placement pass the single-block path runs, over the block list.
# Every axis is block-prefixed here, so a moved block is re-crossed on its
# own fields and its neighbours keep theirs; a prior whose blocks the
# family table does not cover declines with a reason rather than leaving
# the fit with no placement record at all.
multi_rescue <- .nl_registry_grid_rescue(
res, NULL, prior,
refit = function(prior_i) {
r <- solve_grid_multi(prior_i)
if (isTRUE(keep_grid_hessians)) r <- .nl_attach_grid_hessians(r, p_fixed, keep_joint_prec)
.nl_attach_pareto_k(r, prior_i, cargs_no_ckpt, "multi", NULL,
likelihood, k_samples, compute = diagnose_k,
tail_points = k_tail_points)
},
refit_log_marginal = function(prior_i, theta_mat) {
.nl_dispatch_multi(cargs_no_ckpt, prior_i, likelihood = likelihood,
theta_grid_override = theta_mat)$log_marginal
},
auto = .prov$auto, multi = TRUE,
policy = .nl_recenter_mode(control$auto_recenter))
res <- multi_rescue$res
prior <- multi_rescue$prior
res$outer_grid_placement <- res$outer_grid_placement %||% "fixed"
res <- .nl_inner_skew_at_theta(res, prior, cargs_no_ckpt, "multi", NULL,
likelihood, p_fixed, skew_idx,
compute = diagnose_skew)
res <- .nl_skew_correction_attach(res, p_fixed, skew_correct)
res <- .nl_subspace_debias_attach(
res, sd_cfg,
redispatch = function(req) .nl_dispatch_multi(
cargs_no_ckpt, prior, likelihood = likelihood, debias = req,
theta_grid_override = .nl_refined_theta(res)),
p_fixed = p_fixed, beta_names = colnames(X))
res <- .nl_cila_attach(
res, cila_cfg,
redispatch = function(req) .nl_dispatch_multi(
cargs_no_ckpt, prior, likelihood = likelihood, cila = req,
theta_grid_override = .nl_refined_theta(res)),
p_fixed = p_fixed, beta_names = colnames(X),
remoments = function(r) .nl_posterior_moments_multi(
r, r$blocks, r$axis_offsets, r$theta_grid, within = within_cell))
tm$mark("diagnostics")
res$prior <- prior
res$timing <- tm$timing()
return(.finalize_fit(res, backend = "nested_laplace",
n_fixed = p_fixed, fixed_names = colnames(X),
extra_class = c("tulpa_nested_laplace", "list"),
data = list(y = y, n_trials = n_trials,
model_matrix = X, family = family,
offset = offset, phi = phi)))
}
if (!is.null(likelihood)) {
stop("`likelihood` (model-supplied spec) needs a latent prior block. ",
"Pass `prior` (a block or list of blocks).", call. = FALSE)
}
if (is.null(prior$type)) {
stop("Single-block `prior` must have a `type` field, ",
"or pass a list of blocks for multi-block.", call. = FALSE)
}
type <- tolower(prior$type)
tm$mark("setup")
# Every full-grid solve of this fit -- the first one and the placement pass's
# refit -- is screened at the same tolerance and gated the same way, so no
# reported grid is a pruned one whose screen ranking was unreliable.
solve_grid <- function(prior_i) {
out <- .nl_prune_gate(
.nl_dispatch(type, cargs, prior_i, prior_i),
prune_tol_eff,
function() .nl_dispatch(type, utils::modifyList(cargs, list(prune_tol = 0)),
prior_i, prior_i))
.nl_refine_registry(
out, function(theta_mat, screen = NULL) {
blk2 <- .nl_registry_write_theta(
list(prior_i), theta_mat, colnames(theta_mat))[[1L]]
.nl_dispatch(type, .nl_refine_cargs(cargs_no_ckpt, screen, prune_tol_eff),
blk2, prior_i)
},
prior_i, FALSE, .prov$auto,
consistency = var_of_means_consistency,
adaptive_grid = adaptive_grid,
edge_thresh = adaptive_grid_edge_thresh,
max_passes = adaptive_grid_max_passes,
user = control$axis_refine)
}
res <- solve_grid(prior)
tm$mark("grid")
# Integrate exp(log_marginal) over the outer grid. `log_marginal` is the
# inner marginal in the internal (log-scale) hyperparameter parameterization
# the grid is laid out on, and each cell carries the prior mass of the node it
# sits at: its cell width on that coordinate times the axis's declared density
# (R/hyper_quadrature.R). On an evenly spaced grid the widths are equal and
# this reduces to the plain marginal softmax, so where the nodes fall changes
# the quadrature error and not the measure being integrated.
res <- .nl_attach_outer_integration(res, "outer grid")
res <- .nl_posterior_moments(res, type, within = within_cell)
if (isTRUE(keep_grid_hessians)) {
res <- .nl_attach_grid_hessians(
res, p_fixed, keep_joint_prec,
var_cols = .nl_role_cols(prior, NULL, res$Q_csc_n, "temporal"))
}
tm$mark("postproc")
res <- .nl_attach_pareto_k(res, prior, cargs_no_ckpt, "single", type, NULL,
k_samples, compute = diagnose_k,
tail_points = k_tail_points)
# Auto-recenter a railed default outer axis (the registry
# generalization of the joint-path fix, extended over the family table):
# a single mode-Hessian recenter-and-refit when the fixed default axis
# does not contain its own mode. A no-op (zero extra fit) unless it actually
# did; a PINNED axis always wins (`.nl_axis_is_pinned()`: a grid the caller
# named and neither marked with `auto_grid()` nor equal to the engine's own
# default axis). Which families are covered, and why each of the rest
# declines, is `.NL_REGISTRY_AXIS_FIELD`.
registry_rescue <- .nl_registry_grid_rescue(
res, type, prior,
refit = function(prior_i) {
r <- solve_grid(prior_i)
r <- .nl_attach_outer_integration(r, "outer grid")
r <- .nl_posterior_moments(r, type, within = within_cell)
if (isTRUE(keep_grid_hessians)) {
r <- .nl_attach_grid_hessians(
r, p_fixed, keep_joint_prec,
var_cols = .nl_role_cols(prior_i, NULL, r$Q_csc_n, "temporal"))
}
.nl_attach_pareto_k(r, prior_i, cargs_no_ckpt, "single", type, NULL,
k_samples, compute = diagnose_k,
tail_points = k_tail_points)
},
refit_log_marginal = function(prior_i, theta_mat) {
blk2 <- .nl_registry_write_theta(
list(prior_i), theta_mat, colnames(theta_mat))[[1L]]
.nl_dispatch(type, cargs_no_ckpt, blk2, prior_i)$log_marginal
},
auto = .prov$auto,
policy = .nl_recenter_mode(control$auto_recenter))
res <- registry_rescue$res
prior <- registry_rescue$prior
res$outer_grid_placement <- res$outer_grid_placement %||% "fixed"
res <- .nl_inner_skew_at_theta(res, prior, cargs_no_ckpt, "single", type,
NULL, p_fixed, skew_idx,
compute = diagnose_skew)
res <- .nl_skew_correction_attach(res, p_fixed, skew_correct)
res <- .nl_subspace_debias_attach(
res, sd_cfg,
redispatch = function(req) .nl_dispatch(
type, utils::modifyList(cargs_no_ckpt, list(debias = req)),
.nl_refined_block(res, prior), prior),
p_fixed = p_fixed, beta_names = colnames(X))
res <- .nl_cila_attach(
res, cila_cfg,
redispatch = function(req) .nl_dispatch(
type, utils::modifyList(cargs_no_ckpt, list(cila = req)),
.nl_refined_block(res, prior), prior),
p_fixed = p_fixed, beta_names = colnames(X),
remoments = function(r) .nl_posterior_moments(r, type, within = within_cell))
tm$mark("diagnostics")
res$prior <- prior
res$timing <- tm$timing()
.finalize_fit(res, backend = "nested_laplace",
n_fixed = p_fixed, fixed_names = colnames(X),
extra_class = c("tulpa_nested_laplace", "list"),
data = list(y = y, n_trials = n_trials, model_matrix = X,
family = family, offset = offset, phi = phi))
}
# --- cheap-pass screening knobs ---------------------------------------------
#
# The kernels read a tolerance in [0, 1) and a screening depth of at least one
# Newton step. A value outside either range is a caller error rather than
# something to clamp: a clamped tolerance reads back off the fit as a setting
# that was applied.
.nl_check_prune_tol <- function(x) {
if (length(x) != 1L || !is.numeric(x) || !is.finite(x) || x < 0 || x >= 1) {
stop("`prune_tol` must be a finite numeric in [0, 1).", call. = FALSE)
}
as.numeric(x)
}
.nl_check_screen_iters <- function(x) {
if (length(x) != 1L || !is.numeric(x) || !is.finite(x) || x < 1 ||
x > .Machine$integer.max || x != round(x)) {
stop("`screen_iters` must be a single integer >= 1.", call. = FALSE)
}
as.integer(x)
}
# Gate one screened grid solve. At a zero tolerance nothing was screened and
# the kernel result is returned untouched, so a fit with pruning off takes the
# same path it always did. Otherwise the shared cheap-pass gate decides whether
# the kept set can be trusted, re-running `resolve_full` (the same solve at
# `prune_tol = 0`) when it cannot.
.nl_prune_gate <- function(res, prune_tol, resolve_full) {
if (prune_tol <= 0) return(res)
.joint_prune_safety_gate(res, resolve_full = resolve_full,
fn = "tulpa_nested_laplace()")
}
# Positive-scale hyperparameter grid fields: one log-transform unconstrains
# them, so the outer Pareto-k-hat (`.nested_grid_pareto_k`) is well-defined and
# needs no bounded-parameter Jacobian. A block carrying a bounded axis (e.g. a
# correlation `rho_grid`) is absent from this set and is declined to the
# quadrature-ESS fallback rather than transformed with a guessed support.
.NL_POS_GRID <- c("sigma_grid", "sigma2_grid", "tau_grid", "phi_gp_grid",
"phi_grid", "range_grid", "scale_grid", "sigma_spatial_grid",
"lambda_grid", "kappa_grid")
# Attach res$pareto_k for a single-block, single positive-scale-axis nested fit.
# Leaves res$pareto_k = NA (the diagnostic layer then reports the grid's
# quadrature ESS) for multi-block, multi-axis, or bounded-parameter grids.
# `refit` re-evaluates the inner marginal at a substituted grid through the SAME
# driver `res` came from; run with the RNG restored so the fit is unperturbed.
.nl_attach_pareto_k <- function(res, prior, cargs, dispatch_kind, type,
likelihood, n_samples, compute = TRUE,
tail_points = NULL) {
res$pareto_k <- NA_real_
res$pareto_k_is_ess <- NA_real_
res$pareto_k_scope <- "outer (hyperparameter) Gaussian proposal"
# The outer-integration regime is read off the stored grid
# weights, so it costs nothing and is attached on every path -- including the
# ones that DECLINE the k-hat below, where knowing the grid collapsed (and
# whether against a boundary) is the only outer signal left.
res <- .joint_attach_pareto_k_regime(res)
# Every decline below says which one it was.
decline <- function(reason, detail = NULL) {
.k_attach_declined(res, .k_decline(reason, detail))
}
if (!compute) return(decline("not_requested")) # field present, not computed
blocks <- if (is.list(prior) && is.null(prior$type)) prior else list(prior)
if (length(blocks) != 1L) return(decline("not_applicable", "multi-block grid"))
blk <- blocks[[1L]]
# The multi-block dispatch does not carry a `type`, so resolve it from the block
# itself for the length-1 case (a single latent block wrapped in a list, or the
# model-supplied `likelihood` path). Without this the registry lookup below is
# `.NL_REGISTRY[[NULL]]`, an error rather than NULL.
type <- type %||% tolower(blk$type)
# A fit that relied on the default grid carries no `*_grid` field on the input
# prior, but the realised res$theta_grid is a valid single positive-scale axis
# regardless. Fill the type's defaults (the same ones .nl_dispatch applied) so
# eligibility keys on the actual grid, not on whether the user named it. A
# multi-axis type then exposes >1 grid field and still declines below.
spec <- .NL_REGISTRY[[type]]
if (!is.null(spec) && is.function(spec$defaults)) {
blk <- tryCatch(spec$defaults(blk, cargs), error = function(e) blk)
}
gfs <- grep("_grid$", names(blk), value = TRUE)
gfs <- gfs[vapply(gfs, function(f) is.numeric(blk[[f]]) && length(blk[[f]]) >= 2L,
logical(1))]
if (length(gfs) != 1L) {
return(decline("not_applicable",
if (length(gfs)) "multi-axis grid" else "no resolvable grid axis"))
}
if (!(gfs %in% .NL_POS_GRID)) {
# A bounded axis (a correlation `rho_grid`, say): declining beats
# transforming it with a guessed support -- and it will decline for every
# fit of this family, which is what the reason says.
return(decline("unguessable_axis", gfs))
}
tg <- as.numeric(res$theta_grid)
if (length(tg) != length(res$weights)) {
return(decline("internal_inconsistency", "weights do not match theta_grid"))
}
if (any(tg <= 0)) {
return(decline("internal_inconsistency",
paste0("non-positive node on the log axis ", gfs)))
}
# The batch is the diagnostic's own importance sample, not a grid the caller
# chose, so it runs with the internal-batch flag set: the cell-count advice
# would name a number `control$k_samples` sets and offer a remedy that does
# not reach it (gcol33/tulpa#614).
# Both dispatches read the integrated axes off the fit's own block, the
# multi-block one with the draws as its grid, so the draws are scored against
# the hyperprior the grid carries.
refit <- function(theta_mat) {
blk2 <- blk; blk2[[gfs]] <- as.numeric(theta_mat[, 1L])
lm <- tryCatch(
.nl_with_internal_batch(
if (dispatch_kind == "multi")
.nl_dispatch_multi(cargs, list(blk), likelihood = likelihood,
theta_grid_override = theta_mat)$log_marginal
else .nl_dispatch(type, cargs, blk2, blk)$log_marginal),
error = function(e) rep(-Inf, nrow(theta_mat)))
if (length(lm) != nrow(theta_mat)) rep(-Inf, nrow(theta_mat)) else lm
}
kd <- .with_preserved_seed(tryCatch(
.nested_grid_pareto_k(matrix(log(tg), ncol = 1L), res$weights, refit,
n_samples, tail_points = tail_points),
error = function(e) NULL))
if (is.null(kd)) return(decline("degenerate_proposal", "the scorer errored"))
res$pareto_k <- kd$pareto_k
res$pareto_k_is_ess <- kd$is_ess
res$pareto_k_proposal_source <- kd$proposal_source %||% NA_character_
res$pareto_k_first_pass <- kd$first_pass_k %||% NA_real_
# The GPD tail the shape was fitted on. A raised or capped tail was
# invisible on this path (gcol33/tulpa#692).
res$pareto_k_tail_points <- kd$tail_points %||% NA_integer_
.k_attach_declined(res, kd)
}
# --- Per-prior registry & dispatcher ---
#
# Each entry describes one supported prior `type`. To add a new prior:
# 1. Add a row here with cpp_fn, defaults, pack, theta.
# 2. Make sure cpp_nested_laplace_<variant> exists in RcppExports.
# No changes needed to .nl_dispatch() or the public driver.
# --- Required block fields, per dispatch path -------------------------------
#
# Each registry entry declares which fields its converter indexes, keyed by the
# path that consumes them. The C++ kernels read those vectors positionally with
# no bounds check, so a block naming a field wrongly -- a typo, a stale name
# after a rename, a block copied from a different family -- arrives with a
# zero-length vector and takes the session down. Checking
# presence at the R boundary turns that into a message naming the fields.
#
# Path keys:
# axis read by defaults() / theta() while the outer grid is built, so
# checked on every path.
# single extra fields pack() indexes for the single-block kernel.
# multi extra fields .nl_block_spec_for_cpp() indexes.
# joint extra fields .joint_block_spec_for_cpp() indexes.
# joint_single
# extra fields .joint_call_kernel_via_multi() indexes on the
# single-block joint path, where the per-observation index lives on
# each arm rather than on the block.
# projected
# the fields a projected block carries on a multi-block path: an areal
# block's graph or an NNGP block's neighbour structure, with
# `projector` in place of `spatial_idx`.
# A key absent from an entry means the type is not dispatchable on that path;
# the path's own converter raises its "only supported inside ..." message.
.NL_REQ_AREAL_GRAPH <- c("n_spatial_units", "adj_row_ptr", "adj_col_idx",
"n_neighbors")
.NL_REQ_AREAL <- c("spatial_idx", .NL_REQ_AREAL_GRAPH)
.NL_REQ_TEMPORAL <- c("temporal_idx", "n_times")
.NL_REQ_NNGP_FIELD <- c("coords", "nn_idx", "nn_dist", "n_spatial", "nn")
.NL_REQ_HSGP_ARM <- c("m_total", "phi", "n_obs_per_arm", "eigenvalues")
.NL_REQ_SPDE_MESH <- c("n_mesh", "A_x", "A_i", "A_p",
"C0_diag", "G1_x", "G1_i", "G1_p", "nu")
# `paths` is one or more keys above; their union is checked. A field that is
# present but empty counts as missing -- a zero-length vector is exactly what
# the kernel indexes out of bounds.
.nl_check_block_fields <- function(p, paths, block_index = NULL) {
type <- tolower(p$type %||% "")
spec <- .NL_REGISTRY[[type]]
if (is.null(spec) || is.null(spec$required)) return(invisible(p))
req <- unique(unlist(spec$required[paths], use.names = FALSE))
if (!length(req)) return(invisible(p))
bad <- req[vapply(req, function(f) length(p[[f]]) == 0L, logical(1))]
if (!length(bad)) return(invisible(p))
stop("prior block ", if (is.null(block_index)) "" else paste0(block_index, " "),
"'", type, "' is missing required field(s): ",
paste(bad, collapse = ", "), ".", call. = FALSE)
}
# Areal CSR adjacency block -- shared by icar/bym2/car_proper.
.nl_adj_args <- function(p) list(
spatial_idx = as.integer(p$spatial_idx),
n_spatial_units = as.integer(p$n_spatial_units),
adj_row_ptr = as.integer(p$adj_row_ptr),
adj_col_idx = as.integer(p$adj_col_idx),
n_neighbors = as.integer(p$n_neighbors)
)
# Default outer grid for a p-field separable MCAR, returned in column-major
# lower-triangle log-Cholesky coordinates of Sigma = L L' (the columns the
# integrator reads). For the common p = 2 case the nodes are placed in the
# interpretable (sigma_1, sigma_2, rho) space -- with nodes AT sensible
# correlation values -- and converted to log-Cholesky, so the rho axis is
# covered evenly (raw log-Cholesky off-diagonal nodes map to awkward, uneven
# correlations and can miss the mode under sharp likelihoods). For p > 2 a raw
# log-Cholesky tensor is used (coarser; the off-diagonal correlations are not
# separately interpretable).
.mcar_default_logchol_grid <- function(p) {
if (p == 2L) {
sig <- .nl_grid_axis("mcar_sd")
rho <- .nl_grid_axis("mcar_rho")
g <- expand.grid(s1 = sig, s2 = sig, rho = rho,
KEEP.OUT.ATTRS = FALSE)
# Sigma = [[s1^2, rho s1 s2], [rho s1 s2, s2^2]] -> lower Cholesky L:
# L11 = s1, L21 = rho s2, L22 = s2 sqrt(1 - rho^2).
out <- cbind(L11 = log(g$s1),
L21 = g$rho * g$s2,
L22 = log(g$s2 * sqrt(1 - g$rho^2)))
colnames(out) <- c("L11", "L21", "L22")
return(out)
}
diag_nodes <- log(.nl_grid_axis("mcar_logchol_diag"))
off_nodes <- .nl_grid_axis("mcar_logchol_off")
m <- p * (p + 1L) / 2L
axes <- vector("list", m)
nm <- character(m)
t <- 1L
for (j in seq_len(p)) for (i in j:p) {
axes[[t]] <- if (i == j) diag_nodes else off_nodes
nm[t] <- sprintf("L%d%d", i, j)
t <- t + 1L
}
g <- as.matrix(do.call(expand.grid, c(axes, list(KEEP.OUT.ATTRS = FALSE))))
colnames(g) <- nm
g
}
.NL_REGISTRY <- list(
icar = list(
cpp_fn = "cpp_nested_laplace_icar",
required = list(single = .NL_REQ_AREAL, multi = .NL_REQ_AREAL,
joint = .NL_REQ_AREAL,
projected = .NL_REQ_AREAL_GRAPH,
joint_single = .NL_REQ_AREAL_GRAPH),
defaults = function(p, a) .nl_fill_family_axes(p, "icar"),
pack = function(p) c(.nl_adj_args(p), list(
tau_grid = as.numeric(p$tau_grid)
)),
theta = function(p) list(grid = as.numeric(p$tau_grid), names = "tau")
),
bym2 = list(
cpp_fn = "cpp_nested_laplace_bym2",
required = list(single = .NL_REQ_AREAL, multi = .NL_REQ_AREAL,
joint = .NL_REQ_AREAL,
projected = .NL_REQ_AREAL_GRAPH,
joint_single = .NL_REQ_AREAL_GRAPH),
defaults = function(p, a) .nl_fill_family_axes(p, "bym2"),
pack = function(p) c(.nl_adj_args(p), list(
scale_factor = as.numeric(p$scale_factor %||% 1.0),
sigma_spatial_grid = as.numeric(p$sigma_grid),
rho_grid = as.numeric(p$rho_grid),
node_prec = if (!is.null(p$node_prec)) as.numeric(p$node_prec)
)),
theta = function(p) list(
grid = cbind(sigma = p$sigma_grid, rho = p$rho_grid),
names = c("sigma", "rho")
)
),
car_proper = list(
cpp_fn = "cpp_nested_laplace_car_proper",
required = list(single = .NL_REQ_AREAL, multi = .NL_REQ_AREAL,
joint = .NL_REQ_AREAL,
projected = .NL_REQ_AREAL_GRAPH,
joint_single = .NL_REQ_AREAL_GRAPH),
defaults = function(p, a) {
# `rho_car_grid` is the joint-API spelling of the correlation axis; accept
# it as an alias so a grid passed under either name is honoured.
p$rho_grid <- p$rho_grid %||% p$rho_car_grid
# Both axes supplied -> treat as the pre-paired integration grid (the
# theta builder pairs them with cbind), matching bym2 / ar1.
if (!is.null(p$tau_grid) && !is.null(p$rho_grid)) {
return(.nl_check_paired_axes(p, c("tau_grid", "rho_grid")))
}
# Otherwise default each missing axis independently and cross the two, so
# supplying only one axis keeps it instead of discarding it.
rb <- p$rho_bounds %||% .nl_grid_par("car_rho", "bounds")
# Margin in from the eigenvalue endpoints -- Q goes singular at the
# boundary, so anchoring the grid in (lower + eps, upper - eps) avoids
# NaN log-determinants from the per-grid-point Cholesky.
eps <- .nl_grid_par("car_rho", "margin") * (rb[2] - rb[1])
g_tau <- if (is.null(p$tau_grid))
.nl_grid_axis("car_tau") else sort(unique(p$tau_grid))
g_rho <- if (is.null(p$rho_grid))
seq(rb[1] + eps, rb[2] - eps,
length.out = .nl_grid_par("car_rho", "n")) else sort(unique(p$rho_grid))
gr <- expand.grid(tau = g_tau, rho = g_rho)
p$tau_grid <- gr$tau; p$rho_grid <- gr$rho
p
},
pack = function(p) c(.nl_adj_args(p), list(
tau_grid = as.numeric(p$tau_grid),
rho_grid = as.numeric(p$rho_grid)
)),
theta = function(p) list(
grid = cbind(tau = p$tau_grid, rho = p$rho_grid),
names = c("tau", "rho")
)
),
mcar = list(
# Separable multivariate CAR (p coupled areal fields sharing Sigma (x) Q^-1).
# Multi-block-only (like iid): the coupled inner solve lives in the joint
# driver. The outer axes are the p(p+1)/2 log-Cholesky coordinates of Sigma.
cpp_fn = NULL,
required = list(axis = "n_fields",
joint = c(.NL_REQ_AREAL, "n_fields", "field_weight")),
defaults = function(p, a) {
if (is.null(p$logchol_grid)) {
p$logchol_grid <- .mcar_default_logchol_grid(as.integer(p$n_fields))
}
p
},
pack = function(p) stop(
"MCAR is only supported inside a multi-block prior (a coupled areal ",
"field declared by spatial(graph, ~ ... | cell)).", call. = FALSE
),
theta = function(p) {
g <- p$logchol_grid
list(grid = g, names = colnames(g))
}
),
miid = list(
# Multivariate IID: p coupled per-group coefficient
# fields sharing a free Sigma -- the non-spatial sibling of mcar (Q = I, so
# no graph). Multi-block-only (like mcar / iid): the coupled inner solve
# lives in the joint driver. The outer axes are the same p(p+1)/2
# log-Cholesky coordinates of Sigma; only the per-field precision (I vs
# D - W) differs, which does not change the outer grid, so the default grid
# is reused from mcar.
cpp_fn = NULL,
required = list(axis = "n_fields",
joint = c("obs_idx", "n_groups", "n_fields", "field_weight")),
defaults = function(p, a) {
if (is.null(p$logchol_grid)) {
p$logchol_grid <- .mcar_default_logchol_grid(as.integer(p$n_fields))
}
p
},
pack = function(p) stop(
"miid is only supported inside a multi-block prior (a correlated ",
"random-slope term, e.g. (1 + x | group)).", call. = FALSE
),
theta = function(p) {
g <- p$logchol_grid
list(grid = g, names = colnames(g))
}
),
nngp = list(
cpp_fn = "cpp_nested_laplace_nngp",
# spatial_idx / nn_order / cov_type are defaulted on the single-block path,
# so they are not required there; the multi-block path reads spatial_idx.
required = list(single = .NL_REQ_NNGP_FIELD,
multi = c(.NL_REQ_NNGP_FIELD, "spatial_idx"),
joint = c(.NL_REQ_NNGP_FIELD, "spatial_idx"),
projected = .NL_REQ_NNGP_FIELD),
defaults = function(p, a) .nl_fill_family_axes(p, "nngp"),
pack = function(p) {
# spatial_idx (1-based, length N) maps obs -> spatial unit. Default to
# 1..n_spatial when N == n_spatial (the legacy one-obs-per-location case).
spatial_idx <- if (!is.null(p$spatial_idx)) {
as.integer(p$spatial_idx)
} else {
seq_len(as.integer(p$n_spatial))
}
c(list(spatial_idx = spatial_idx),
.nl_nngp_field_args(p),
list(sigma2_grid = as.numeric(p$sigma2_grid),
phi_gp_grid = as.numeric(p$phi_gp_grid)))
},
theta = function(p) list(
grid = cbind(sigma2 = p$sigma2_grid, phi_gp = p$phi_gp_grid),
names = c("sigma2", "phi_gp")
)
),
spde = list(
# Continuous Matern SPDE field on a shared FEM mesh. INDEXED_MULTI: each
# observation projects onto ~3 mesh nodes through the FEM matrix A, so the
# block is consumed via obs_indices (the A CSC slots), not a one-node
# spatial_idx. Same SpdeQBuilder and make_spde_block factory the
# single-Laplace occupancy SPDE path (cpp_nested_laplace_spde) uses.
cpp_fn = NULL, # multi-block only; reached via .nl_dispatch_multi
# The single-block converter reads a scalar `n_obs`; the joint one reads a
# per-arm `n_obs_per_arm` alongside per-arm A CSC triples.
required = list(multi = c(.NL_REQ_SPDE_MESH, "n_obs"),
joint = c(.NL_REQ_SPDE_MESH, "n_obs_per_arm")),
defaults = function(p, a) {
if (is.null(p$range_grid) || is.null(p$sigma_grid)) {
# Geometric grid centred on the PC-prior mode (prior_range[1] /
# prior_sigma[1]). The span is kept tight (mode / SPAN .. mode * SPAN)
# rather than the wide [mode*0.3, mode*3] window: a rectangular
# Cartesian grid (no CCD mode-find here) that runs to a very small
# range / large sigma corner lets the binary-occupancy field over-fit
# there, and the corner cell -- being the highest inner likelihood --
# dominates the weighted field even though the PC prior disfavours it.
# Anchoring the grid near the regularised prior mode keeps the
# integrated field in the well-identified region the single-Laplace
# SPDE path fits at the fixed mode, while still integrating modest
# hyperparameter uncertainty.
n_g <- as.integer(p$n_grid %||% .nl_grid_par("spde_registry", "n"))
span <- as.numeric(p$grid_span %||% .nl_grid_par("spde_registry", "span"))
rng_mode <- (.hp_range_anchor(p)$anchor %||%
.nl_grid_par("spde_registry", "centre_range"))[1]
sig_mode <- (p$prior_sigma %||% .nl_scale_anchor())[1]
rg <- exp(seq(log(rng_mode / span), log(rng_mode * span), length.out = n_g))
sg <- exp(seq(log(sig_mode / span), log(sig_mode * span), length.out = n_g))
gr <- expand.grid(range = rg, sigma = sg)
p$range_grid <- gr$range
p$sigma_grid <- gr$sigma
}
p
},
pack = function(p) stop(
"spde is only supported inside a multi-block latent prior on the ",
"nested-Laplace path. Pass it as a one-element list to ",
"tulpa_nested_laplace(prior = list(blk), ...).",
call. = FALSE
),
theta = function(p) list(
grid = cbind(range = as.numeric(p$range_grid),
sigma = as.numeric(p$sigma_grid)),
names = c("range", "sigma")
)
),
iid = list(
# Unstructured (iid) Gaussian latent, integrated over a single sigma axis.
# Multi-block-only: there is no single-arm iid kernel (a lone iid field is a
# plain random intercept). Reached via .nl_dispatch_multi -- the C++
# build_blocks_from_spec builds the block with d_fac(k) = sigma_k and an
# N(0, 1) prior. A one-point sigma_grid conditions the field at a FIXED
# sigma (an EM-estimated random-effect term on a nested block carries its
# own grid this way; see .fit_block_via_nested_laplace).
cpp_fn = NULL,
required = list(multi = c("obs_idx", "n_units"),
joint = c("obs_idx", "n_units")),
defaults = function(p, a) .nl_fill_family_axes(p, "iid"),
pack = function(p) stop(
"iid is only supported inside a multi-block latent prior on the ",
"nested-Laplace path. Pass it as one element of a block list to ",
"tulpa_nested_laplace(prior = list(...), ...).",
call. = FALSE
),
theta = function(p) list(grid = as.numeric(p$sigma_grid), names = "sigma")
),
hsgp = list(
cpp_fn = "cpp_nested_laplace_hsgp",
# The single-block kernel takes the basis as one matrix (phi_basis) with its
# eigenvalues (lambda_eig); the joint one takes a per-arm list (phi) with
# `eigenvalues` and the basis / row counts it needs to slice them.
required = list(single = c("phi_basis", "lambda_eig"),
multi = c("phi_basis", "lambda_eig"),
joint = .NL_REQ_HSGP_ARM),
defaults = function(p, a) .nl_fill_family_axes(p, "hsgp"),
pack = function(p) list(
phi_basis = as.matrix(p$phi_basis),
lambda_eig = as.numeric(p$lambda_eig),
sigma2_grid = as.numeric(p$sigma2_grid),
lengthscale_grid = as.numeric(p$lengthscale_grid)
),
theta = function(p) list(
grid = cbind(sigma2 = p$sigma2_grid, lengthscale = p$lengthscale_grid),
names = c("sigma2", "lengthscale")
)
),
hsgp_mo = list(
# Multi-output (co-regionalization) HSGP block. K = n_arms
# correlated latent fields share basis + eigenvalues, cross-output
# coefficients live in Sigma. First ship: K = 2 with raw axes
# (sigma_1, sigma_2, rho, ell) and the natural Cartesian product over
# paired per-axis grids. Multi-block-only -- there is no single-arm
# version of multi-output HSGP.
cpp_fn = NULL,
required = list(joint = .NL_REQ_HSGP_ARM),
defaults = function(p, a) {
# All four grids treated as paired axes (each row of theta_grid is
# one (sigma_1, sigma_2, rho, ell) tuple). When any grid is missing
# we fill in a small Cartesian default; when the user supplies grids
# they must already be aligned (use expand.grid() on the R side to
# form a Cartesian, then pass the four aligned columns).
.nl_fill_family_axes(p, "hsgp_mo")
},
pack = function(p) stop(
"hsgp_mo is only supported inside a multi-block joint prior. ",
"Wrap the block in list() and pass to tulpa_nested_laplace_joint(prior = list(blk), ...).",
call. = FALSE
),
theta = function(p) list(
grid = cbind(sigma_1 = as.numeric(p$sigma_1_grid),
sigma_2 = as.numeric(p$sigma_2_grid),
rho = as.numeric(p$rho_grid),
ell = as.numeric(p$lengthscale_grid)),
names = c("sigma_1", "sigma_2", "rho", "ell")
)
),
rw1 = list(
cpp_fn = "cpp_nested_laplace_temporal",
required = list(single = .NL_REQ_TEMPORAL, multi = .NL_REQ_TEMPORAL,
joint = .NL_REQ_TEMPORAL),
defaults = function(p, a) .nl_fill_family_axes(p, "rw1"),
pack = function(p) list(
temporal_idx = as.integer(p$temporal_idx),
n_times = as.integer(p$n_times),
temporal_type = "rw1",
tau_grid = as.numeric(p$tau_grid),
rho_grid = numeric(0),
cyclic = isTRUE(p$cyclic),
n_groups = as.integer(p$n_groups %||% 1L)
),
theta = function(p) list(grid = as.numeric(p$tau_grid), names = "tau")
),
rw2 = list(
cpp_fn = "cpp_nested_laplace_temporal",
required = list(single = .NL_REQ_TEMPORAL, multi = .NL_REQ_TEMPORAL,
joint = .NL_REQ_TEMPORAL),
defaults = function(p, a) .nl_fill_family_axes(p, "rw2"),
pack = function(p) list(
temporal_idx = as.integer(p$temporal_idx),
n_times = as.integer(p$n_times),
temporal_type = "rw2",
tau_grid = as.numeric(p$tau_grid),
rho_grid = numeric(0),
cyclic = FALSE,
n_groups = as.integer(p$n_groups %||% 1L)
),
theta = function(p) list(grid = as.numeric(p$tau_grid), names = "tau")
),
ar1 = list(
cpp_fn = "cpp_nested_laplace_temporal",
required = list(single = .NL_REQ_TEMPORAL, multi = .NL_REQ_TEMPORAL,
joint = .NL_REQ_TEMPORAL),
defaults = function(p, a) {
# The temporal kernel takes rho as given, so an AR1 correlation outside
# (-1, 1) -- a non-stationary process -- is refused here, by name.
.nl_check_axis_open_interval(p, "rho_grid", -1, 1,
"an AR1 correlation")
.nl_fill_family_axes(p, "ar1")
},
pack = function(p) list(
temporal_idx = as.integer(p$temporal_idx),
n_times = as.integer(p$n_times),
temporal_type = "ar1",
tau_grid = as.numeric(p$tau_grid),
rho_grid = as.numeric(p$rho_grid),
cyclic = FALSE,
n_groups = as.integer(p$n_groups %||% 1L)
),
theta = function(p) list(
grid = cbind(tau = p$tau_grid, rho = p$rho_grid),
names = c("tau", "rho")
)
),
lf = list(
# Latent factor block. F = 1 only. No outer-grid axes:
# identifiability is handled by tight Gaussian anchors on (u_1,
# lambda_1) inside the C++ factory; everything else is jointly
# optimized by the inner Newton. Returns a 1 x 0 grid so the joint
# Cartesian product treats this block as a "no axis" contributor.
cpp_fn = NULL,
required = list(joint = c("obs_idx", "n_latent")),
defaults = function(p, a) p,
pack = function(p) stop(
"lf is only supported inside a multi-block prior. ",
"Wrap the block in list(): tulpa_nested_laplace_joint(prior = list(blk), ...).",
call. = FALSE
),
theta = function(p) list(
grid = matrix(numeric(0), nrow = 1L, ncol = 0L),
names = character(0)
)
),
tgmrf = list(
# User-defined GMRF block (see ?tgmrf). Only available inside a multi-
# block prior. The outer theta-grid is the Cartesian product of
# per-axis grids; per-axis defaults are 5 equispaced points between
# block$bounds$lower[j] and block$bounds$upper[j], or block$init[j] +/-
# 2 if no bounds are supplied. Precomputed Q_k / logdet / log p(theta_k)
# at every grid row are packed into the C++ block spec by
# `.nl_block_spec_for_cpp`.
cpp_fn = NULL,
# The multi-block converter builds the per-grid Q / logdet / log-prior from
# the block's own closures; the joint one takes them already materialised.
required = list(
axis = c("theta_dim", "theta_names"),
multi = c("n_latent", "obs_idx"),
joint = c("n_latent", "obs_idx", "Q_csc_p_per_grid", "Q_csc_i_per_grid",
"Q_csc_x_per_grid", "logdet_Q_per_grid",
"log_prior_theta_per_grid")),
defaults = function(p, a) {
if (is.null(p$theta_grid_built)) {
d <- p$theta_dim
axes <- vector("list", d)
for (j in seq_len(d)) {
if (!is.null(p$bounds)) {
lo <- p$bounds$lower[j]
hi <- p$bounds$upper[j]
} else {
hw <- .nl_grid_par("tgmrf_axis", "half_width")
lo <- p$init[j] - hw
hi <- p$init[j] + hw
}
axes[[j]] <- seq(lo, hi,
length.out = .nl_grid_par("tgmrf_axis", "n"))
}
names(axes) <- p$theta_names
gr <- do.call(expand.grid, axes)
p$theta_grid_built <- as.matrix(gr)
}
p
},
pack = function(p) stop(
"tgmrf is only supported inside a multi-block prior. ",
"Wrap the block in list(): tulpa_nested_laplace(prior = list(blk), ...).",
call. = FALSE
),
theta = function(p) list(
grid = p$theta_grid_built,
names = p$theta_names
)
)
)
.nl_dispatch <- function(type, a, p, declared) {
spec <- .NL_REGISTRY[[type]]
if (is.null(spec)) {
stop("Unknown prior type: ", type,
". Supported: ", paste(names(.NL_REGISTRY), collapse = ", "), ".",
call. = FALSE)
}
# A multi-block-only type has no single-arm kernel; let its own pack() raise
# the "only supported inside a multi-block prior" message rather than a
# field-name list for fields the single path would never read.
if (!is.null(spec$cpp_fn)) .nl_check_block_fields(p, c("axis", "single"))
hyperprior <- .hp_choice(a$hyperprior)
a$hyperprior <- NULL
# A refinement round's screen record (`.hyper_refine_screen()`): the cells
# are a batch joining a solved grid. Read here, never handed to a kernel.
screen <- a$screen_refine
a$screen_refine <- NULL
p <- spec$defaults(p, a)
th <- spec$theta(p)
tg <- .nl_theta_matrix(list(theta_grid = th$grid, theta_names = th$names))
# `declared` is the block whose grid the fit integrates: `p` itself for a grid
# solve, the fit's block when `p` carries a stencil or a probe row written onto
# its grid fields. The integrated columns are read off that grid, since the
# rows written onto `p` vary a column the fit holds fixed or hold one it
# integrates at a single value (gcol33/tulpa#760).
dth <- spec$theta(spec$defaults(declared, a))
hp_declared <- .hp_declare(
.nl_theta_matrix(list(theta_grid = dth$grid, theta_names = dth$names)))
# The block's hyperprior, folded where every caller of the kernel reads it:
# the grid solve, the placement refit and its stencil, the k-hat refit. A
# screened solve ranks its cells with it too.
hp <- list(.nl_block_log_hyperprior(p, tg, hyperprior, declared = hp_declared))
if ((a$prune_tol %||% 0) > 0) {
a$screen_log_offset <- .nl_screen_log_offset(tg, hp, screen = screen)
}
out <- do.call(spec$cpp_fn, c(spec$pack(p), a))
out$theta_grid <- th$grid
out$theta_names <- th$names
.nl_fold_hyperprior(out, hp)
}
# Per-cell log weight the cheap screen ranks with: the hyperprior R folds into
# `log_marginal` after the kernel returns, plus the cell measure the grid's
# weights are built with (`.nl_attach_outer_integration()`), in the grid's cell
# order. Both are fixed by the grid before any cell is solved, and together with
# the kernel's log-marginal they are the posterior weight a cell ends up with,
# so the screen and the full pass rank the same quantity.
#
# `screen` is the record a refinement round hands its kernel
# (`.hyper_refine_screen()`): the cells are a batch joining a solved grid, so
# their measure is read on the merged grid (`log_measure`) rather than on the
# batch, and the grid's own mass (`log_ref`) travels as the offset's `log_ref`
# attribute, which the kernel's screen normalises against.
.nl_screen_log_offset <- function(tg, hp_parts, specs = NULL, screen = NULL) {
tg <- as.matrix(tg)
n <- nrow(tg)
lp <- numeric(n)
axes <- character(0)
for (hp in hp_parts) {
lp <- lp + hp$lp
axes <- union(axes, hp$axes)
}
lq <- if (!is.null(screen)) {
if (length(screen$log_measure) != n) {
stop("A refinement screen measures ", length(screen$log_measure),
" cells; the kernel call holds ", n, ".", call. = FALSE)
}
screen$log_measure
} else {
if (is.null(specs) && !is.null(colnames(tg))) {
specs <- .joint_axis_specs_from_grid(tg, folded_axes = axes)
}
.hyper_log_quad_weights(tg, specs)
}
off <- lp + (if (length(lq) == n) lq else 0)
off[is.na(off)] <- -Inf
if (!is.null(screen)) attr(off, "log_ref") <- as.numeric(screen$log_ref)
off
}
# The outer grid's cell measure, its weights and its log evidence, set together
# from the grid and the log-marginals the fit carries. `log_marginal` already
# holds the hyperprior density the dispatch folded in (`R/hyperprior_default.R`),
# so a folded axis is measured by its cell widths alone and its support is not
# clipped to the span the nodes were declared over.
.nl_attach_outer_integration <- function(res, what = "outer grid") {
tg <- .nl_theta_matrix(res)
# A node a refinement pass appended was not declared, so the specs are read
# off the base cells and the cells' tags carry the rest.
refining <- res$refining_axis
base <- !nzchar(.hyper_slice_home(refining, nrow(tg)))
specs <- .joint_axis_specs_from_grid(
tg[base, , drop = FALSE], folded_axes = res$log_hyperprior_axes)
res$log_quad <- .hyper_log_quad_weights(tg, specs, refining = refining)
res$axis_support <- .hyper_grid_supports(tg, specs, refining = refining)
res$weights <- .nl_normalise_weights_safe(res$log_marginal, what,
log_quad = res$log_quad)
res <- .nl_prune_merged_record(res)
.nl_attach_evidence(res, tg, specs)
}
# Normalise log-marginals to integration weights summing to 1. The single
# weight normaliser for every outer hyperparameter grid: it drops non-finite
# nodes (an inner Newton diverging in a grid corner returns +Inf / NaN) from
# the max-shift, zeroes their weight, and renormalises over the finite cells;
# all-NA (with a warning) when no node carries mass, never NaN. There is no
# unguarded variant -- a single `max(lm)` over a grid with one +Inf / NaN cell
# would collapse every weight to NaN and break the grid sampler. `what` names
# the grid in the degenerate-case warning.
# `log_quad` adds the per-cell log quadrature weight of the outer grid (see
# R/hyper_quadrature.R), so posterior weight is likelihood times prior mass
# rather than likelihood alone. Callers integrating over a hyperparameter grid
# pass it; callers weighting something that is not a quadrature rule leave it
# NULL.
#' Normalise outer-grid log-marginals to posterior cell weights
#'
#' The single weight normaliser for every outer hyperparameter grid: softmax
#' of `lm` (optionally plus a per-cell log quadrature weight) into weights
#' summing to 1. Non-finite nodes (an inner Newton diverging in a grid
#' corner) are dropped from the max-shift and zeroed before renormalising
#' over the finite cells; an all-non-finite grid returns all-`NA` with a
#' warning rather than `NaN`. Intended for reconstructing a fit's outer-grid
#' posterior weights (with [tulpa_theta_matrix()] and
#' [tulpa_grid_log_quad()]) when the fit doesn't already carry `fit$weights`.
#'
#' @param lm Numeric vector of per-cell log-marginal-likelihood values.
#' @param what Label for the grid, used only in the degenerate-case warning.
#' @param log_quad Optional per-cell log quadrature weight (prior mass) of
#' the outer grid, e.g. from [tulpa_grid_log_quad()]; `NULL` normalises the
#' likelihood alone.
#' @return Numeric vector of posterior cell weights summing to 1, the same
#' length as `lm`, or all-`NA` if no cell carries finite mass.
#' @seealso [tulpa_theta_matrix()], [tulpa_grid_log_quad()]
#' @examples
#' tulpa_normalise_weights_safe(c(-1, 0, -2))
#' # With per-cell log quadrature weights (cell prior masses):
#' tulpa_normalise_weights_safe(c(-1, 0, -2), log_quad = log(c(0.2, 0.5, 0.3)))
#' @export
tulpa_normalise_weights_safe <- function(lm, what = "grids / data",
log_quad = NULL) {
.nl_normalise_weights_safe(lm, what = what, log_quad = log_quad)
}
.nl_normalise_weights_safe <- function(lm, what = "grids / data",
log_quad = NULL) {
if (!is.null(log_quad)) {
if (length(log_quad) != length(lm)) {
stop(sprintf("`log_quad` has length %d but `lm` has length %d (%s).",
length(log_quad), length(lm), what), call. = FALSE)
}
lm <- lm + log_quad
lm[is.na(lm)] <- -Inf
}
finite_lm <- lm[is.finite(lm)]
if (length(finite_lm) == 0L) {
warning(sprintf("All grid points returned non-finite log_marginal -- check %s.", what),
call. = FALSE)
return(rep(NA_real_, length(lm)))
}
w <- exp(lm - max(finite_lm))
w[!is.finite(w)] <- 0
if (sum(w) == 0) return(rep(NA_real_, length(lm)))
w / sum(w)
}
# Log evidence of an outer hyperparameter grid.
#
# `log_marginal` carries a normalised hyperprior density on each axis's
# integration coordinate and `log_measure` is each cell's ABSOLUTE log volume on
# that coordinate (its widths, and a declared atom's probability), so
#
# log p(y) = log sum_k exp(log_marginal_k + log_measure_k)
#
# is the evidence under the prior itself, not under a copy of it restricted to
# the grid and renormalised there. Where the grid covers the posterior the
# value does not move with where the nodes were laid or how many there are.
#
# A cell whose inner solve failed (`log_marginal` NaN, or `+Inf` from a non-PD
# Hessian) contributes nothing: its likelihood is unknown and is not guessed.
# A cell screened out by the pruning pass carries `-Inf` and contributes
# nothing either, which is what screening it asserted.
.nl_outer_log_evidence <- function(log_marginal, log_measure) {
lq <- as.numeric(log_measure)
if (length(lq) != length(log_marginal)) {
stop(sprintf(paste0("Outer evidence: `log_marginal` has length %d but the ",
"measure has %d."),
length(log_marginal), length(lq)), call. = FALSE)
}
known <- !is.na(log_marginal) & log_marginal < Inf & !is.na(lq)
if (!any(known)) return(NA_real_)
ev <- .tulpa_logsumexp(log_marginal[known] + lq[known])
if (is.finite(ev)) ev else NA_real_
}
# Record a grid's evidence on `res`, or why it has none. The evidence exists
# only when every integrated axis carries a normalised density; otherwise it is
# `NA` with `log_evidence_declined = "improper_hyperprior"` and
# `log_evidence_declined_axes` naming each such axis and its reason.
.nl_attach_evidence <- function(res, tg, specs,
log_measure = .hyper_log_quad_weights(
tg, specs, absolute = TRUE,
refining = res$refining_axis)) {
uncovered <- .nl_evidence_uncovered(tg, specs, res)
res$log_evidence_declined_axes <- if (length(uncovered)) uncovered else NULL
if (length(uncovered)) {
res$log_evidence <- NA_real_
res$log_evidence_declined <- "improper_hyperprior"
return(res)
}
n <- length(res$log_marginal)
if (!is.null(tg) && nrow(as.matrix(tg)) == n &&
!length(.hp_integrated_axes(as.matrix(tg)))) {
# Every axis is pinned, so each is part of the model and the grid integrates
# nothing: its cells are copies of one point, which together carry measure 1.
log_measure <- rep(-log(n), n)
}
res$log_evidence <- .nl_outer_log_evidence(res$log_marginal, log_measure)
res$log_evidence_declined <- NULL
res
}
`%||%` <- function(x, y) if (is.null(x)) y else x
# ----------------------------------------------------------------------------
# Per-grid-point fixed-effects Hessian extraction.
#
# When tulpa_nested_laplace() is called with keep_grid_hessians = TRUE the
# inner Laplace solves are run with store_Q = TRUE, so each grid point's full
# posterior precision (over the joint latent vector x = (beta, RE, latent))
# is returned in CSC lower-triangle form. Here we
#
# 1. Reconstruct sparse symmetric Q_k from the CSC triple at grid point k.
# 2. Compute Sigma_{beta beta} = (Q_k^{-1})[1:p, 1:p] via a single sparse
# solve with the first p columns of the identity as RHS.
# 3. Invert that small p x p block to get the marginal fixed-effects
# Hessian H_beta = Sigma_{beta beta}^{-1}.
#
# Per-grid memory cost: O(p^2). The dominant compute is the sparse solve,
# already required for the integrated log-marginal -- so adding it here is
# essentially free.
# ----------------------------------------------------------------------------
# `keep_joint` additionally retains the MODAL cell's joint precision as
# `res$H_joint`, which is what the subspace-debias coupling closure reads to
# grow the corrected set over strongly coupled neighbours
# (`.subspace_closure()`). One cell, not `n_grid`: the closure is a SELECTION
# decision, and the selection is made once. Off by default, so a fit that did
# not ask for closure retains nothing extra (gcol33/tulpa#862).
# The MODAL cell's joint posterior precision, assembled from the per-cell CSC
# scratch the inner solves return under `store_Q`. This is what the
# subspace-debias coupling closure reads (`.subspace_closure()`): a coordinate
# strongly coupled to a member of the corrected set and left OUT of it is being
# carried linearly, which is precisely the error the correction removes, so the
# closure needs the precision graph to find it.
#
# ONE cell, not `n_grid`. The closure is a selection decision made once, and
# the modal cell is the one the reported mixture is centred on. Both drivers
# read this -- the single/multi-block grid through
# `.nl_attach_grid_hessians()`, the joint through `.joint_finalize_grid_fixed()`
# -- so the two cannot disagree about which cell the coupling came from
# (gcol33/tulpa#862). NULL when the fit retained no precision, which the caller
# records as a declined closure rather than as "no neighbours found".
.nl_modal_joint_precision <- function(res) {
Q_p <- res$Q_csc_p_per_grid
Q_i <- res$Q_csc_i_per_grid
Q_x <- res$Q_csc_x_per_grid
n_x <- res$Q_csc_n
if (is.null(Q_p) || is.null(Q_i) || is.null(Q_x) || is.null(n_x)) return(NULL)
lm <- res$log_marginal
if (is.null(lm) || !any(is.finite(lm))) return(NULL)
k <- which.max(replace(lm, !is.finite(lm), -Inf))
if (k > length(Q_p) || length(Q_p[[k]]) != n_x + 1L) return(NULL)
L <- Matrix::sparseMatrix(i = Q_i[[k]], p = Q_p[[k]], x = Q_x[[k]],
dims = c(n_x, n_x), index1 = FALSE)
L + Matrix::t(L) - Matrix::Diagonal(n_x, Matrix::diag(L))
}
#
# `var_cols` names latent coordinates whose per-cell MARGINAL variance,
# diag(Q_k^{-1}), is retained alongside (`grid_field_var`, n_grid x
# length(var_cols), with the columns in `grid_field_var_cols`). The same solve
# delivers it: those unit columns join E. It defaults to the temporal field (the
# `temporal =` block, or a latent(temporal_ar*()) block), so
# temporal() reports each cell's Gaussian -- mode AND curvature -- mixed over
# the grid rather than the cell modes alone, which carried only the spread
# between cells and covered the truth at 57% for a nominal 95%
# (gcol33/tulpa#904). Only a field whose read-back needs it is retained: one
# column per coordinate per cell, and a large spatial field is not asked for.
#
# The field's full within-cell covariance (`grid_field_cov`, a list of
# n_var x n_var matrices parallel to the cells, NULL where a cell was pruned)
# is retained too, so temporal()'s draws carry the correlation BETWEEN time
# points a cell's Gaussian implies and not only each point's marginal: a
# derived quantity over time (a trend, a difference between two years, a
# cumulative sum) read off draws with independent coordinates has the wrong
# spread. Capped at `.NL_FIELD_COV_MAX` doubles over the grid; past the cap only
# the marginals are kept and the draws fall back to independent coordinates.
.NL_FIELD_COV_MAX <- 4e6
.nl_attach_grid_hessians <- function(res, p_fixed, keep_joint = FALSE,
var_cols = c(
.nl_role_cols(res$blocks,
res$block_latent_offsets,
res$Q_csc_n, "temporal"),
.nl_temporal_latent_block(res)$cols)) {
Q_p <- res$Q_csc_p_per_grid
Q_i <- res$Q_csc_i_per_grid
Q_x <- res$Q_csc_x_per_grid
n_x <- res$Q_csc_n
if (is.null(Q_p) || is.null(Q_i) || is.null(Q_x) || is.null(n_x)) {
warning("keep_grid_hessians = TRUE but the inner solve did not return Q. ",
"This is a tulpa bug; please report.",
call. = FALSE)
return(res)
}
n_grid <- length(Q_p)
grid_hessians <- vector("list", n_grid)
grid_modes <- vector("list", n_grid)
var_cols <- as.integer(var_cols)
var_cols <- var_cols[var_cols > p_fixed & var_cols <= n_x]
n_var <- length(var_cols)
field_var <- if (n_var) matrix(NA_real_, n_grid, n_var) else NULL
# The field's whole within-cell covariance is the lower-right block of the
# same S, so retaining it costs no further solve -- only the storage, which
# grows as n_var^2 per cell and is capped by `.NL_FIELD_COV_MAX`.
keep_cov <- n_var > 0L && as.numeric(n_grid) * n_var^2 <= .NL_FIELD_COV_MAX
field_cov <- if (keep_cov) vector("list", n_grid)
E <- matrix(0, nrow = n_x, ncol = p_fixed + n_var)
E[cbind(c(seq_len(p_fixed), var_cols), seq_len(p_fixed + n_var))] <- 1
modes_mat <- res$modes # n_grid x n_x
# An intrinsic block (ICAR, BYM2's structured part, RW1 / RW2) stands for a
# field whose constant is a HARD sum-to-zero constraint. The inner solve
# augments that direction so the Hessian inverts, and centres it into the
# intercept, so the stored precision carries a finite 1 / (tau J) variance
# on the level -- and the unconditioned inverse puts it into the intercept's
# standard error (3-4x the exact posterior SD on a 50-node ICAR,
# gcol33/tulpa#901) and into an intrinsic temporal field's marginal variance.
# The driver names each intrinsic block's constraint group
# (`constraint_cols`), and both reads are conditioned on them by the kriging
# correction a joint fit's retained block and posterior draws use
# (cpp_joint_inner_vcov_blocks, Rue & Held 2005):
# Cov(x | A x = 0) = V - W (A W)^-1 W', V = Q^-1 E, W = Q^-1 A'.
# The mode needs no correction: the centring already put it on the
# constraint.
constr <- res$constraint_cols
n_con <- length(constr)
if (n_con) {
E <- cbind(E, vapply(constr, function(cc) {
a <- numeric(n_x)
a[as.integer(cc)] <- 1
a
}, numeric(n_x)))
}
sel <- seq_len(p_fixed + n_var)
for (k in seq_len(n_grid)) {
# A cell the cheap screen pruned is never solved, so it holds no precision
# to retain and its CSC slot is empty. Skip the write and the slot stays
# blank, which is what `.nested_fixed_moments()` already reads: it drops any
# cell whose block is absent, and a pruned cell carries log_marginal = -Inf
# and hence zero weight, so nothing the report reads was in it. Skipping
# rather than assigning NULL is what keeps `grid_modes` / `grid_hessians`
# parallel to `weights` -- assigning NULL removes the element and shortens
# the list, which NAs the whole coefficient table (gcol33/tulpa#345).
if (length(Q_p[[k]]) != n_x + 1L) next
grid_modes[[k]] <- if (!is.null(modes_mat)) {
as.numeric(modes_mat[k, seq_len(p_fixed)])
} else {
rep(NA_real_, p_fixed)
}
L <- Matrix::sparseMatrix(
i = Q_i[[k]], p = Q_p[[k]], x = Q_x[[k]],
dims = c(n_x, n_x), index1 = FALSE
)
diag_L <- Matrix::diag(L)
Qk <- L + Matrix::t(L) - Matrix::Diagonal(n_x, diag_L)
# V = Qk^{-1} E selects the fixed-effect (and read field) columns of the
# joint inverse. A dense GP / NNGP field precision at close locations makes
# the joint precision ill-conditioned, which trips Matrix's
# condition-number guard. Retry with a negligible diagonal jitter (scaled
# to the precision) so the marginal-SE extraction stays available rather
# than aborting the whole fit.
V <- tryCatch(
as.matrix(Matrix::solve(Qk, E)),
error = function(e) {
jit <- 1e-8 * mean(Matrix::diag(Qk))
as.matrix(Matrix::solve(Qk + Matrix::Diagonal(n_x, jit), E))
}
)
rows <- c(seq_len(p_fixed), var_cols)
S <- V[rows, sel, drop = FALSE]
if (n_con) {
W <- V[, p_fixed + n_var + seq_len(n_con), drop = FALSE]
AW <- crossprod(E[, p_fixed + n_var + seq_len(n_con), drop = FALSE], W)
Wr <- W[rows, , drop = FALSE]
S <- S - Wr %*% solve((AW + t(AW)) / 2, t(Wr))
}
Sigma_bb <- S[seq_len(p_fixed), seq_len(p_fixed), drop = FALSE]
grid_hessians[[k]] <- .spd_inverse(Sigma_bb)
if (n_var) {
fi <- p_fixed + seq_len(n_var)
field_var[k, ] <- diag(S)[fi]
if (keep_cov) {
Sf <- S[fi, fi, drop = FALSE]
field_cov[[k]] <- (Sf + t(Sf)) / 2
}
}
}
res$grid_hessians <- grid_hessians
res$grid_modes <- grid_modes
if (n_var) {
res$grid_field_var <- field_var
res$grid_field_var_cols <- var_cols
if (keep_cov) res$grid_field_cov <- field_cov
}
if (isTRUE(keep_joint)) res$H_joint <- .nl_modal_joint_precision(res)
# Strip the verbose CSC scratch fields once Hessians are assembled.
res$Q_csc_p_per_grid <- NULL
res$Q_csc_i_per_grid <- NULL
res$Q_csc_x_per_grid <- NULL
res$Q_csc_n <- NULL
res
}
# --- Multi-block dispatch ---
#
# A multi-block prior is `list(block1, block2, ...)` where each blockN is a
# list with a `type` field (same shape as a single-block prior). At each
# outer-grid point all blocks share a Newton solve; the joint grid is the
# Cartesian product of per-block axes.
# Block types whose field may reach the observations through a projector: arm
# k's predictor receives A_k z in place of the unit gather z[spatial_idx]
# (restricted spatial regression, A_k = P_k S_k). Both multi-block drivers read
# one; an areal block's level is placed from A_k 1, an NNGP block's prior
# carries its own.
.NL_PROJECTABLE_AREAL <- c("icar", "bym2", "car_proper")
.NL_PROJECTABLE_TYPES <- c(.NL_PROJECTABLE_AREAL, "nngp")
.nl_block_is_projected <- function(p) {
is.list(p) && tolower(p$type %||% "") %in% .NL_PROJECTABLE_TYPES &&
!is.null(p$projector)
}
# A block `projector` is read by the two multi-block drivers, the single-arm
# one behind tulpa_nested_laplace() and the joint one. The single-block joint
# path would gather the field through `spatial_idx` and drop the projector, so
# it refuses it.
.nl_refuse_projector <- function(prior, door) {
blocks <- if (is.list(prior) && !is.null(prior$type)) list(prior) else prior
if (!is.list(blocks)) return(invisible(NULL))
has <- vapply(blocks, function(b) is.list(b) && !is.null(b$projector),
logical(1))
if (any(has)) {
stop(door, " does not read a block `projector`; a projected areal field ",
"is fitted by tulpa_nested_laplace_joint() with a list of blocks, ",
"`prior = list(<block>)`.", call. = FALSE)
}
invisible(NULL)
}
# A projected block's per-arm projectors as the column-compressed triple
# the C++ block builders read (0-based `A_i` / `A_p`, the SPDE projector's
# layout). `p$projector` is a list of n_arms matrices, each N_k x n_units
# (dense or a Matrix sparse matrix), or one matrix shared by every arm; a NULL
# entry is an arm the field does not reach. The projector replaces the unit
# gather, so a block carrying `spatial_idx` as well, or a per-row
# `svc_weight` that belongs inside A_k, is refused.
.nl_block_projector <- function(p, n_arms, block_index, arm_n_obs) {
type <- tolower(p$type)
lab <- paste0("Block ", block_index, " (type '", type, "')")
if (!is.null(p$spatial_idx)) {
stop(lab, ": pass either `projector` or `spatial_idx`, not both; a ",
"projector already maps the observations to the units.",
call. = FALSE)
}
if (!is.null(p$svc_weight)) {
stop(lab, ": `projector` and `svc_weight` cannot be combined; scale ",
"the projector's rows by the weight.", call. = FALSE)
}
if (is.null(arm_n_obs)) {
stop(lab, ": a projector needs the arms' observation counts.",
call. = FALSE)
}
n_units <- as.integer(p$n_spatial_units %||% p$n_spatial)
projector <- p$projector
if (!is.list(projector)) projector <- rep(list(projector), n_arms)
if (length(projector) != n_arms) {
stop(lab, ": `projector` must be a list of length n_arms (", n_arms,
"), or one matrix shared by every arm.", call. = FALSE)
}
csc <- lapply(seq_len(n_arms), function(k) {
A <- projector[[k]]
if (is.null(A)) {
return(list(x = numeric(0), i = integer(0), p = integer(n_units + 1L)))
}
if (!(is.matrix(A) || methods::is(A, "Matrix")) ||
nrow(A) != arm_n_obs[k] || ncol(A) != n_units) {
stop(lab, ": `projector[[", k, "]]` must be a ", arm_n_obs[k],
" x ", n_units, " matrix (arm ", k, "'s observations by the ",
"block's units).", call. = FALSE)
}
A <- methods::as(methods::as(Matrix::Matrix(A, sparse = TRUE),
"generalMatrix"), "CsparseMatrix")
if (!all(is.finite(A@x))) {
stop(lab, ": `projector[[", k, "]]` holds a non-finite entry.",
call. = FALSE)
}
list(x = as.numeric(A@x), i = as.integer(A@i), p = as.integer(A@p))
})
list(A_x = lapply(csc, `[[`, "x"),
A_i = lapply(csc, `[[`, "i"),
A_p = lapply(csc, `[[`, "p"))
}
# An NNGP block's neighbour structure as the C++ block builders read it, on
# either the single-block kernel or the multi-block driver. nn_order is passed
# 0-based, the convention cpp_laplace_fit_gp also uses.
.nl_nngp_field_args <- function(p) {
n_spatial <- as.integer(p$n_spatial)
list(
coords = as.matrix(p$coords),
nn_idx = as.matrix(p$nn_idx),
nn_dist = as.matrix(p$nn_dist),
nn_order = as.integer((p$nn_order %||% seq_len(n_spatial)) - 1L),
n_spatial = n_spatial,
nn = as.integer(p$nn),
cov_type = as.integer(p$cov_type %||% 2L) # default Matern-5/2
)
}
.is_multi_block_prior <- function(p) {
is.list(p) && is.null(p$type) &&
all(vapply(p, function(x) is.list(x) && !is.null(x$type), logical(1)))
}
# What each block of a nested prior is FOR: "spatial", "temporal" or "other".
# tulpa() tags its `spatial =` / `temporal =` blocks with `role`, and on a tagged
# prior every untagged block (an s(x) smoother, a latent() block, a `(1 | g)`
# iid block) is "other" -- a type does not say what a block is for, since a
# smoother is an rw1 / rw2 block and a random intercept an iid one. A prior
# built by hand carries no tags, and its blocks are classified by type
# (`.SPATIAL_NL_TYPES` / `.TEMPORAL_NL_TYPES`). `prior` is a single block or a
# list of blocks; the result has one entry per block.
.nl_block_roles <- function(prior) {
if (is.null(prior) || !is.list(prior) || !length(prior)) return(character(0))
blocks <- if (!is.null(prior$type)) list(prior) else prior
tagged <- vapply(blocks, function(b) {
r <- if (is.list(b)) b$role else NULL
if (is.character(r) && length(r) == 1L) r else NA_character_
}, character(1))
if (any(!is.na(tagged))) return(ifelse(is.na(tagged), "other", tagged))
types <- vapply(blocks, function(b) tolower(b$type %||% ""), character(1))
ifelse(types %in% .SPATIAL_NL_TYPES, "spatial",
ifelse(types %in% .TEMPORAL_NL_TYPES, "temporal", "other"))
}
# Columns of a nested fit's latent vector (`modes`, the per-cell precision)
# holding the first block whose role is `role`, or NULL. A multi-block solve
# reports where each block starts (`block_latent_offsets`, from the C++ block
# builders, which are what decide a block's width); a single-block solve lays
# out [fixed, RE, field], so its field is the tail. `blocks` is the fit's block
# list (or single block), `n_x` the latent length.
.nl_role_cols <- function(blocks, offsets, n_x, role) {
roles <- .nl_block_roles(blocks)
k <- which(roles == role)[1L]
if (is.na(k) || is.null(n_x)) return(NULL)
if (length(roles) == 1L && !is.null(blocks$type)) {
w <- .nl_single_block_width(blocks)
if (is.null(w) || w < 1L || w > n_x) return(NULL)
return(as.integer(seq.int(n_x - w + 1L, n_x)))
}
if (length(offsets) != length(roles) + 1L) return(NULL)
lo <- offsets[k]; hi <- offsets[k + 1L]
if (!is.finite(lo) || !is.finite(hi) || hi <= lo || hi > n_x) return(NULL)
as.integer(seq.int(lo + 1L, hi))
}
# Width of a single temporal block's field (the only single-block field whose
# coefficients a read-back slices; the rest report through their own
# accessors). NULL when the block does not say.
.nl_single_block_width <- function(p) {
nt <- p$n_times
if (is.null(nt)) return(NULL)
as.integer(nt) * max(as.integer(p$n_groups %||% 1L), 1L)
}
# Structural spec for one block sent to cpp_nested_laplace_multi. Grid values
# go in theta_grid separately; this function returns only the type-specific
# data needed to build a LatentBlock (idx vector, sizes, adjacency, etc.).
#
# `block_joint_grid` is the joint-grid sub-matrix (n_joint x block$theta_dim)
# already realised from the Cartesian product. tgmrf needs it to precompute
# Q(theta_k) at every joint-grid row; built-in types ignore the argument
# because they read theta values directly from theta_grid in the C++ side.
.nl_block_spec_for_cpp <- function(p, block_joint_grid = NULL, n_obs = NULL,
block_index = NULL) {
type <- tolower(p$type)
projected <- .nl_block_is_projected(p)
.nl_check_block_fields(p, c("axis", if (projected) "projected" else "multi"),
block_index)
# Optional per-observation design weight, attached after the type branch has
# built the spec: it turns any block into a varying coefficient whose eta
# contribution is svc_weight[i] * <the block's field at i>, and the C++ side
# honours it on every contribution kind. Reading it per type would decide
# which fields can vary by where the read happens to be written.
.with_svc <- function(out) {
if (is.null(p$svc_weight)) return(out)
out$svc_weight <- as.numeric(p$svc_weight)
out
}
if (type %in% .NL_PROJECTABLE_AREAL) {
out <- list(
type = type,
spatial_idx = if (!projected) as.integer(p$spatial_idx),
n_spatial_units = as.integer(p$n_spatial_units),
adj_row_ptr = as.integer(p$adj_row_ptr),
adj_col_idx = as.integer(p$adj_col_idx),
n_neighbors = as.integer(p$n_neighbors)
)
if (projected) {
out$projector <- .nl_block_projector(p, 1L, block_index %||% 1L, n_obs)
}
if (type == "bym2") {
out$scale_factor <- as.numeric(p$scale_factor %||% 1.0)
if (!is.null(p$node_prec)) out$node_prec <- as.numeric(p$node_prec)
}
.with_svc(out)
} else if (type == "nngp") {
out <- c(list(type = "nngp",
spatial_idx = if (!projected) as.integer(p$spatial_idx)),
.nl_nngp_field_args(p))
if (projected) {
out$projector <- .nl_block_projector(p, 1L, block_index %||% 1L, n_obs)
}
.with_svc(out)
} else if (type == "hsgp") {
.with_svc(list(
type = "hsgp",
phi_basis = as.matrix(p$phi_basis),
lambda_eig = as.numeric(p$lambda_eig)
))
} else if (type %in% c("rw1", "rw2", "ar1")) {
out <- list(
type = type,
temporal_idx = as.integer(p$temporal_idx),
n_times = as.integer(p$n_times)
)
if (type %in% c("rw1", "rw2")) out$cyclic <- isTRUE(p$cyclic)
.with_svc(out)
} else if (type == "iid") {
.with_svc(list(
type = "iid",
obs_idx = as.integer(p$obs_idx),
n_units = as.integer(p$n_units)
))
} else if (type == "spde") {
# FEM mesh + projection. The C++ side wraps the single A CSC in a length-1
# per-arm list and builds the LatentBlock via make_spde_block. Optional
# rational coefficients (fractional nu) pass through when present.
out <- list(
type = "spde",
n_mesh = as.integer(p$n_mesh),
n_obs = as.integer(p$n_obs),
A_x = as.numeric(p$A_x),
A_i = as.integer(p$A_i),
A_p = as.integer(p$A_p),
C0_diag = as.numeric(p$C0_diag),
G1_x = as.numeric(p$G1_x),
G1_i = as.integer(p$G1_i),
G1_p = as.integer(p$G1_p),
nu = as.numeric(p$nu)
)
if (!is.null(p$rational_poles) && !is.null(p$rational_weights)) {
out$rational_poles <- as.numeric(p$rational_poles)
out$rational_weights <- as.numeric(p$rational_weights)
}
.with_svc(out)
} else if (type == "tgmrf") {
# Precompute Q(theta_k) for every joint-grid row and pack CSC triples +
# logdet + log p(theta_k). The C++ side reads them directly at grid
# point k - no R callback during the inner Newton loop. Iterating over
# block_joint_grid (rather than the block's own grid) is what lets
# tgmrf compose with other blocks in a multi-block prior: when other
# blocks vary across joint rows but tgmrf does not, the duplicated Q
# matrices land at the right joint index.
#
# The Q / log_prior(theta) evaluation either goes through the user's R
# closure (`backend == "r"`, the default) or through the registered C++
# spec (`backend == "cpp"`). Wire format downstream is identical in both
# cases.
if (is.null(block_joint_grid)) {
stop("Internal error: tgmrf block_spec called without block_joint_grid.",
call. = FALSE)
}
tg <- as.matrix(block_joint_grid)
n_cells <- nrow(tg)
n_latent <- p$n_latent
Q_csc_p <- vector("list", n_cells)
Q_csc_i <- vector("list", n_cells)
Q_csc_x <- vector("list", n_cells)
logdet_Q <- numeric(n_cells)
log_pi_th <- numeric(n_cells)
use_cpp <- isTRUE(identical(p$backend, "cpp"))
if (use_cpp) {
if (is.null(p$cpp_id) || !nzchar(p$cpp_id)) {
stop("Internal error: tgmrf block has backend = 'cpp' ",
"but no cpp_id.", call. = FALSE)
}
if (!cpp_tgmrf_registry_has(p$cpp_id)) {
stop("tgmrf C++ spec '", p$cpp_id, "' is not registered in the ",
"current R session. Was the user .cpp file sourceCpp'd? ",
"(Re-load via tgmrf_cpp(cpp_file = ..., id = '", p$cpp_id,
"').)", call. = FALSE)
}
}
for (k in seq_len(n_cells)) {
th <- as.numeric(tg[k, ])
names(th) <- p$theta_names
if (use_cpp) {
evk <- cpp_tgmrf_eval(p$cpp_id, th)
if (as.integer(evk$n) != n_latent) {
stop("Registered tgmrf_cpp spec '", p$cpp_id,
"' returned Q of size ", evk$n, " at grid row ", k,
" but the block was registered with n_latent = ", n_latent,
".", call. = FALSE)
}
Q_csc_p[[k]] <- as.integer(evk$p)
Q_csc_i[[k]] <- as.integer(evk$i)
Q_csc_x[[k]] <- as.numeric(evk$x)
Qk <- Matrix::sparseMatrix(
i = Q_csc_i[[k]], p = Q_csc_p[[k]], x = Q_csc_x[[k]],
dims = c(n_latent, n_latent), index1 = FALSE
)
lp <- as.numeric(evk$log_prior)
} else {
Qk <- p$Q(th)
Qk <- methods::as(methods::as(Qk, "generalMatrix"), "CsparseMatrix")
if (nrow(Qk) != n_latent || ncol(Qk) != n_latent) {
stop("`Q(theta)` returned ", nrow(Qk), " x ", ncol(Qk),
" at grid row ", k,
" but the block was registered with n_latent = ", n_latent, ".",
call. = FALSE)
}
Q_csc_p[[k]] <- as.integer(Qk@p)
Q_csc_i[[k]] <- as.integer(Qk@i)
Q_csc_x[[k]] <- as.numeric(Qk@x)
lp <- as.numeric(p$prior(th))
}
# sparse Cholesky log-determinant; falls back to dense det for tiny
# blocks where Matrix::Cholesky balks.
ld <- tryCatch({
ch <- Matrix::Cholesky(methods::as(Qk, "symmetricMatrix"), LDL = FALSE)
# log|Q| = 2 * sum(log(diag(L))) for the LL' factor.
2 * Matrix::determinant(ch, sqrt = TRUE, logarithm = TRUE)$modulus
}, error = function(e) {
as.numeric(determinant(Qk, logarithm = TRUE)$modulus)
})
logdet_Q[k] <- as.numeric(ld)
if (!is.finite(lp)) {
stop("`prior(theta)` returned a non-finite value at grid row ", k,
" (theta = ", paste(format(th, digits = 4), collapse = ", "),
").", call. = FALSE)
}
log_pi_th[k] <- lp
}
list(
type = "tgmrf",
n_latent = as.integer(n_latent),
obs_idx = as.integer(p$obs_idx),
Q_csc_p_per_grid = Q_csc_p,
Q_csc_i_per_grid = Q_csc_i,
Q_csc_x_per_grid = Q_csc_x,
logdet_Q_per_grid = as.numeric(logdet_Q),
log_prior_theta_per_grid = as.numeric(log_pi_th)
)
} else {
stop("Block type '", type, "' is not supported in multi-block priors.",
call. = FALSE)
}
}
# Per-block axis grid as a matrix `[n_block_cells x n_axes_for_block]`.
# Single-axis types (icar / rw1 / rw2 / iid) get a 1-column matrix; 2-axis
# types (bym2 / car_proper / ar1) get a 2-column matrix.
.nl_block_axis_grid <- function(p) {
spec <- .NL_REGISTRY[[tolower(p$type)]]
if (is.null(spec)) {
stop("Unknown block type: ", p$type, call. = FALSE)
}
.nl_check_block_fields(p, "axis")
p <- spec$defaults(p, list())
th <- spec$theta(p)
if (is.matrix(th$grid)) {
grid <- th$grid
if (is.null(colnames(grid))) colnames(grid) <- th$names
} else {
grid <- matrix(as.numeric(th$grid), ncol = 1)
colnames(grid) <- th$names
}
list(grid = grid, names = colnames(grid), prepared = p)
}
# Wall-clock threshold, in seconds, above which a multi-block outer-grid
# solve gets a "this was slow" warning. Cell count alone does not predict
# solve time -- a family/likelihood pair, N, or a block's own per-cell cost
# (gcol33/tulpa#638) can make the same cell count cheap or expensive, and the
# engine's own default per-block grids already exceed any static cell-count
# threshold that fires on genuinely slow custom grids (gcol33/tulpa#820) --
# so the trigger is measured cost, read off the "grid" bucket both multi-block
# dispatchers already mark with `.tulpa_timer()`.
.NL_MULTI_GRID_WARN_SECONDS <- 2
# Shared by the single-response multi-block dispatch (`tulpa_nested_laplace()`
# in this file) and the joint multi-block dispatch
# (`R/nested_laplace_joint_multi.R`): both solve their outer grid in one call
# and mark its elapsed time as the "grid" timing bucket, so both warn off that
# measurement through this one function rather than keeping two copies of the
# same message. `layout` is `.nl_grid_crossing()`'s account of the axes that
# produced the count, so the advice names which grids to shrink.
.nl_multi_grid_warn <- function(elapsed, n_cells, remedy, layout = "") {
if (isTRUE(elapsed > .NL_MULTI_GRID_WARN_SECONDS) && !.nl_internal_batch()) {
warning(sprintf(
"Multi-block outer grid (%s cells) took %s to solve.%s %s",
format(n_cells, scientific = FALSE, trim = TRUE),
.format_duration(elapsed), layout, remedy
), call. = FALSE)
}
invisible(NULL)
}
# Cell count above which a dense multi-block outer tensor is announced with a
# warning BEFORE its solve starts. Each cell is one inner Laplace solve, so the
# count is the run's size in solves; it does not know the per-solve cost, which
# is why it warns rather than refuses and why the post-solve warning above reads
# measured time.
.NL_DENSE_GRID_WARN_CELLS <- 2048L
# Shared by both multi-block dispatchers: the dense tensor about to be solved
# has `n_cells` inner solves, `layout` (`.nl_grid_crossing()`) names the axes
# that produced them and `remedy` is the dispatcher's advice. Gated like the
# timing warning, so an internal batch the caller did not size stays quiet.
.nl_dense_grid_warn <- function(n_cells, remedy, layout = "") {
if (isTRUE(n_cells > .NL_DENSE_GRID_WARN_CELLS) && !.nl_internal_batch()) {
warning(sprintf(
paste0("Multi-block outer grid is a dense tensor of %s cells, one inner ",
"solve each.%s %s"),
format(n_cells, scientific = FALSE, trim = TRUE), layout, remedy
), call. = FALSE)
}
invisible(NULL)
}
# Whether the dispatch currently running is an INTERNAL batch rather than the
# grid the caller asked for. The soft-cap warning above is advice addressed to
# whoever chose the cell count, and on an internal batch that is not the caller:
# the outer Pareto-k diagnostic re-evaluates `log_marginal` at `k_samples`
# importance draws by substituting them for the block's grid axis, so a 7-node
# fit warned about a cell count the diagnostic's own budget set and advised
# reducing per-block grid sizes, which does not reach the number in the message
# (gcol33/tulpa#614).
#
# Transported on an option rather than an argument because the two warning sites
# sit behind different dispatchers -- `.nl_dispatch_multi()` takes `cargs`, the
# joint multi-block driver does not -- which is the same reason
# `tulpa.nl_progress` and `tulpa.nl_checkpoint` are options.
.nl_internal_batch <- function() {
isTRUE(getOption("tulpa.nl_internal_batch", FALSE))
}
# Run `expr` with that flag set. Nests and restores, so an internal batch that
# itself re-dispatches stays internal.
.nl_with_internal_batch <- function(expr) {
op <- options(tulpa.nl_internal_batch = TRUE)
on.exit(options(op), add = TRUE)
force(expr)
}
# The caller's ceiling on a dense multi-block outer tensor, in inner solves.
# Unset there is none: a cell count does not measure the run's cost, so the
# default is the pre-solve warning above, and a caller who wants a hard stop
# (a batch of fits on a shared node, say) sets `control$max_grid_cells`. The
# front doors publish the resolved value on the scoped
# `tulpa.nl_max_grid_cells` option, so every grid built inside the fit -- the
# multi-block dispatch, the joint multi-block dispatch, the k-hat
# re-evaluations -- reads one ceiling without carrying it in each signature,
# the transport `tulpa.nl_progress` and `tulpa.nl_checkpoint` already use.
.nl_max_grid_cells <- function(control = list()) {
v <- control$max_grid_cells
if (is.null(v)) v <- getOption("tulpa.nl_max_grid_cells", NULL)
if (is.null(v)) return(Inf)
if (!is.numeric(v) || length(v) != 1L || is.na(v) || v < 1) {
stop("`control$max_grid_cells` must be a single number >= 1 ",
"(the cell-count ceiling on a multi-block outer grid).", call. = FALSE)
}
as.numeric(v)
}
# The single enforcement of that ceiling, shared by both multi-block
# dispatchers, before any inner solve. `layout` is `.nl_grid_crossing()`'s
# account of the axes that produced the count and `remedy` the dispatcher's
# advice.
.nl_check_grid_cap <- function(n_cells, max_cells, remedy, layout = "") {
if (n_cells <= max_cells) return(invisible(n_cells))
fmt <- function(x) format(x, scientific = FALSE, trim = TRUE)
stop(sprintf(
paste0("Multi-block outer grid has %s cells, past control$max_grid_cells ",
"= %s.%s %s"),
fmt(n_cells), fmt(max_cells), layout, remedy
), call. = FALSE)
}
# The axes a dense outer tensor crosses, as the sentence the timing warning
# carries: " It crosses b1 (...) x phi_pos 4." `block_grids` are the per-block
# axis grids, so a fit that set no grid learns which default axes produced its
# count (gcol33/tulpa#913); `phi_axes` are the active named per-arm dispersion
# axes crossed on top, each a factor on the number of inner solves.
.nl_grid_crossing <- function(block_grids = NULL, phi_axes = NULL) {
parts <- c(if (length(block_grids)) .nl_grid_layout(block_grids),
if (length(phi_axes))
paste0("phi_", names(phi_axes), " ", lengths(phi_axes)))
if (length(parts)) paste0(" It crosses ", paste(parts, collapse = " x "), ".")
else ""
}
# One block's grid as "b<k> (<rows> rows: <axis> <levels> x ...)", the levels
# being each axis's distinct node count. A block grid is a tensor of its axes
# for the grid-axis families and a paired node list for others (BYM2's
# (sigma, rho) nodes), so the row count is stated alongside the per-axis
# levels rather than derived from them.
.nl_grid_layout <- function(block_grids) {
parts <- vapply(seq_along(block_grids), function(b) {
g <- as.matrix(block_grids[[b]])
nm <- colnames(g) %||% paste0("axis", seq_len(ncol(g)))
lv <- vapply(seq_len(ncol(g)), function(j)
length(unique(as.numeric(g[, j]))), integer(1))
sprintf("b%d (%d rows: %s)", b, nrow(g),
paste(nm, lv, collapse = " x "))
}, character(1))
paste(parts, collapse = " x ")
}
.nl_dispatch_multi <- function(cargs, prior_list, likelihood = NULL,
progress = .nl_progress_args(list(progress = FALSE)),
within_cell = NULL,
# Single-cell re-dispatch for the inner-Laplace
# skewness diagnostic (.nl_inner_skew_at_theta()):
# when supplied, `theta_grid_override` (a 1-row
# matrix at the outer grid's MAP cell) replaces
# the Cartesian-product grid construction below,
# reusing the SAME blocks_spec machinery so the
# re-solve is at exactly the fitted MAP theta.
# The SPDE-prior fold / weight normalisation /
# posterior-moments tail is skipped in that case
# -- the caller only wants inner_skew, not a
# re-normalised single-cell "posterior".
theta_grid_override = NULL,
compute_skew = FALSE, skew_idx = NULL,
debias = NULL,
# Corrected integrated Laplace:
# the kernel-facing request list, off by default.
cila = NULL) {
# Inject a default obs_idx for any tgmrf block that didn't supply one.
# The C++ scatter needs obs_idx[i] -> latent slot for each observation;
# the canonical "one obs per latent slot" case has N == n_latent and the
# identity mapping. Heterogeneous (N != n_latent) layouts must supply
# obs_idx explicitly at tgmrf() construction.
N_obs <- length(cargs$y)
for (b in seq_along(prior_list)) {
blk <- prior_list[[b]]
if (identical(tolower(blk$type %||% ""), "tgmrf") && is.null(blk$obs_idx)) {
if (blk$n_latent != N_obs) {
stop("tgmrf block has n_latent = ", blk$n_latent,
" but data has N = ", N_obs, " observations. Provide an explicit ",
"`obs_idx` argument to tgmrf() mapping observations to latent slots.",
call. = FALSE)
}
prior_list[[b]]$obs_idx <- seq_len(N_obs)
}
}
per_block <- lapply(prior_list, .nl_block_axis_grid)
prepared <- lapply(per_block, function(x) x$prepared)
block_grids <- lapply(per_block, function(x) x$grid)
axis_counts <- vapply(block_grids, ncol, integer(1))
axis_offsets <- as.integer(c(0L, cumsum(axis_counts)))
# Block-prefix the axis names so they're unambiguous (e.g. block1.tau).
axis_names <- unlist(lapply(seq_along(block_grids), function(b) {
paste0("b", b, ".", colnames(block_grids[[b]]))
}))
# Set only on the Cartesian-product branch below, so the post-solve timing
# check (which reads these) is a no-op on a single-cell override re-dispatch
# (skew / k-hat / recenter probes -- already never the caller's own grid).
n_cells <- NA_integer_
grid_warn_remedy <- NULL
grid_layout <- ""
if (!is.null(theta_grid_override)) {
joint_grid <- as.matrix(theta_grid_override)
total_axes <- axis_offsets[length(axis_offsets)]
if (ncol(joint_grid) != total_axes) {
stop("theta_grid_override must have ", total_axes,
" columns (sum of block axes), got ", ncol(joint_grid), ".",
call. = FALSE)
}
colnames(joint_grid) <- axis_names
} else {
# Cartesian product of per-block row indices. This path has no CCD
# alternative (unlike the joint multi-block dispatch), so the only
# remedy for a slow or oversized grid is a smaller one.
row_counts <- vapply(block_grids, nrow, integer(1))
idx <- do.call(expand.grid, lapply(row_counts, seq_len))
n_cells <- nrow(idx)
grid_warn_remedy <- "Reduce per-block grid sizes."
grid_layout <- .nl_grid_crossing(block_grids)
.nl_check_grid_cap(n_cells, .nl_max_grid_cells(), grid_warn_remedy,
grid_layout)
.nl_dense_grid_warn(n_cells, grid_warn_remedy, grid_layout)
# Concatenate per-block axis grids into the joint theta_grid.
joint_grid <- do.call(cbind, lapply(seq_along(block_grids), function(b) {
block_grids[[b]][idx[[b]], , drop = FALSE]
}))
colnames(joint_grid) <- axis_names
}
blocks_spec <- lapply(seq_along(prepared), function(b) {
cols <- (axis_offsets[b] + 1L):axis_offsets[b + 1L]
block_joint <- joint_grid[, cols, drop = FALSE]
.nl_block_spec_for_cpp(prepared[[b]], block_joint,
n_obs = length(cargs$y), block_index = b)
})
# Every block's hyperprior over its own columns, folded on the override path
# too: a placement stencil re-evaluates the marginal through it and has to
# difference the target the grid integrates. A screened solve ranks with it.
# The integrated columns are read off the block's declared grid, never off the
# override, whose rows hold a fixed axis at other values or an integrated one
# at a single value (gcol33/tulpa#760).
hp_parts <- lapply(seq_along(prepared), function(b) {
cols <- (axis_offsets[b] + 1L):axis_offsets[b + 1L]
tg_b <- joint_grid[, cols, drop = FALSE]
colnames(tg_b) <- colnames(block_grids[[b]])
blk <- prepared[[b]]
blk$log_prior_theta_per_grid <- blocks_spec[[b]]$log_prior_theta_per_grid
.hp_prefix(.nl_block_log_hyperprior(blk, tg_b, .hp_choice(cargs$hyperprior),
declared = .hp_declare(block_grids[[b]])),
paste0("b", b, "."))
})
prune_tol <- as.numeric(cargs$prune_tol %||% 0)
grid_solve_start <- proc.time()[["elapsed"]]
out <- cpp_nested_laplace_multi(
y = cargs$y,
n = cargs$n,
X = cargs$X,
re_idx = cargs$re_idx,
n_re_groups = cargs$n_re_groups,
sigma_re = cargs$sigma_re,
blocks_spec = blocks_spec,
theta_grid = joint_grid,
axis_offsets = axis_offsets,
family = cargs$family,
phi = cargs$phi,
max_iter = cargs$max_iter,
tol = cargs$tol,
n_threads = cargs$n_threads,
x_init_nullable = cargs$x_init_nullable,
store_Q = isTRUE(cargs$store_Q),
# Cheap-pass screening tolerance; 0 (the default) solves every cell. A
# single-cell re-dispatch inherits it harmlessly -- the driver screens only
# a grid with more than one cell.
prune_tol = prune_tol,
likelihood = likelihood,
progress = isTRUE(progress$progress),
progress_every = as.integer(progress$progress_every),
progress_throttle = as.numeric(progress$progress_throttle),
progress_file = as.character(progress$progress_file),
checkpoint_path = as.character(cargs$checkpoint_path %||% ""),
compute_skew = isTRUE(compute_skew),
skew_idx = skew_idx,
debias = debias,
cila = cila,
offset_nullable = cargs$offset_nullable,
screen_log_offset = if (prune_tol > 0)
.nl_screen_log_offset(joint_grid, hp_parts,
screen = cargs$screen_refine)
)
if (!is.null(grid_warn_remedy)) {
.nl_multi_grid_warn(proc.time()[["elapsed"]] - grid_solve_start,
n_cells, grid_warn_remedy, grid_layout)
}
out$theta_grid <- joint_grid
out$theta_names <- axis_names
out$axis_offsets <- axis_offsets
out$blocks <- prepared
out <- .nl_fold_hyperprior(out, hp_parts)
if (!is.null(theta_grid_override)) return(out) # skew re-dispatch / the
# auto-recenter FD stencil
out <- .nl_attach_outer_integration(out, "multi-block outer grid")
out <- .nl_posterior_moments_multi(out, prepared, axis_offsets, joint_grid,
within = .nl_within_cell_mode(within_cell))
out
}
#' Build a `prior` list for [tulpa_nested_laplace()] from a tulpa spec object
#'
#' @description
#' Validates a `tulpa_temporal` or `tulpa_spatial` specification against
#' `data`, then converts it to the prior list shape consumed by
#' [tulpa_nested_laplace()]. Mainly an internal helper for callers that already
#' have a fitted spec; users typically pass `spec` + `data` directly to
#' `tulpa_nested_laplace()` instead.
#'
#' Supported spec types:
#' * `tulpa_temporal` with `type in {"rw1", "rw2", "ar1"}`
#' * `tulpa_spatial`, areal `type in {"car", "icar", "car_proper", "bym2"}`.
#' An areal spec returns the graph fields only -- `type`, `spatial_idx`,
#' `n_spatial_units`, `adj_row_ptr`, `adj_col_idx`, `n_neighbors` -- plus
#' `scale_factor` for bym2 and the eigenvalue-derived `rho_bounds` for
#' proper CAR. No grid field is returned for any areal type: the outer axes
#' are built later by the family's registry `defaults()` closure, and proper
#' CAR's `(tau, rho)` grid is built there from those bounds.
#' * `tulpa_gp` / `tulpa_hsgp`, continuous `type in {"gp", "nngp", "hsgp"}`:
#' validated against `data` and routed to the `nngp` / `hsgp` nested kernel
#' through the shared converter.
#'
#' SPDE is the one continuous field not built here -- it carries its own
#' (range, sigma) FEM integrator; call [fit_spde()] with the `spatial_spde()`
#' spec.
#'
#' @param spec A `tulpa_temporal` or `tulpa_spatial` object.
#' @param data Data frame the spec resolves time/group/site indices against.
#' @return A `prior` list ready for [tulpa_nested_laplace()].
#' @examples
#' d <- data.frame(year = rep(1:10, each = 3))
#' pr <- prior_from_spec(temporal_ar1("year"), d)
#' str(pr)
#' @export
prior_from_spec <- function(spec, data) {
if (inherits(spec, "tulpa_temporal")) {
return(.prior_from_temporal_spec(spec, data))
}
if (inherits(spec, "tulpa_spatial")) {
return(.prior_from_spatial_spec(spec, data))
}
stop("`spec` must inherit from 'tulpa_temporal' or 'tulpa_spatial'.",
call. = FALSE)
}
.prior_from_temporal_spec <- function(spec, data) {
spec <- validate_temporal(spec, data)
type <- tolower(spec$type)
if (!type %in% c("rw1", "rw2", "ar1")) {
stop("tulpa_nested_laplace() supports temporal types rw1, rw2, ar1; got '",
type, "'.", call. = FALSE)
}
out <- list(
type = type,
temporal_idx = as.integer(spec$time_index),
n_times = as.integer(spec$n_times)
)
if (type %in% c("rw1", "rw2")) out$cyclic <- isTRUE(spec$cyclic)
out
}
.prior_from_spatial_spec <- function(spec, data) {
type <- tolower(spec$type)
# Continuous coordinate-addressed fields (gp/nngp/hsgp). Validate to populate
# the spec's resolved coordinate structure (NNGP neighbour info / HSGP basis),
# then build the nested prior block through the shared converter
# `.spatial_spec_to_nl_prior` -- the single source of truth for spec -> nested
# prior, also used by the tulpa() front door. gp/nngp drive the `nngp` kernel,
# hsgp the `hsgp` kernel; both are single-block (cpp_fn present in .NL_REGISTRY).
if (type %in% c("gp", "nngp")) {
return(.spatial_spec_to_nl_prior(validate_gp(spec, data)))
}
if (type == "hsgp") {
return(.spatial_spec_to_nl_prior(validate_hsgp(spec, data)))
}
if (type == "spde") {
stop("tulpa_nested_laplace() does not host an SPDE field: SPDE carries its ",
"own (range, sigma) FEM integrator. Call fit_spde() with the ",
"spatial_spde() spec.", call. = FALSE)
}
if (!type %in% .NL_FRONTDOOR_AREAL) {
stop("Unknown spatial type '", type, "'. Areal: ",
paste(.NL_FRONTDOOR_AREAL, collapse = "/"),
"; continuous: gp/nngp/hsgp; SPDE via fit_spde().", call. = FALSE)
}
validate_spatial(spec, data)
# Areal field: resolve the per-observation spatial_idx (1-based row indices
# into the adjacency) onto the spec, then build the block through the shared
# converter. Empty cells (no data) are allowed -- the ICAR / CAR / BYM2 prior
# is well-defined on every graph node; unobserved cells contribute no
# likelihood term. See `.resolve_spatial_idx()` for the resolution rules.
adj <- spec$adjacency
n_spatial_units <- nrow(as.matrix(adj))
if (!is.null(spec$level) && spec$level == "group" &&
!is.null(spec$group_var)) {
if (!(spec$group_var %in% names(data))) {
stop("Spatial group variable '", spec$group_var,
"' not found in data.", call. = FALSE)
}
spec$spatial_idx <- .resolve_spatial_idx(
values = data[[spec$group_var]],
n_spatial_units = n_spatial_units,
adjacency = adj,
group_var = spec$group_var
)
} else {
if (nrow(data) != n_spatial_units) {
stop("Observation-level spatial spec requires nrow(data) == ",
"nrow(adjacency).", call. = FALSE)
}
spec$spatial_idx <- seq_len(n_spatial_units)
}
.spatial_spec_to_nl_prior(spec)
}
# Normalize the outer-grid progress knobs from a `control` list into the four
# scalars the C++ entry points accept. On by default; `progress = FALSE` turns
# off the flushed cell-k/n_grid + ETA reporter (see the GridProgress reporter in
# inst/include/tulpa/nested_progress.h). `progress.file`
# adds a heartbeat file for detached runs where Rcout flushing is unreliable.
# The inner refinement / EM / CCD-probe call sites pass `progress = FALSE`
# explicitly so only the top-level fit ticks.
.nl_progress_args <- function(control) {
# Control keys win; otherwise fall back to the scoped `tulpa.nl_progress`
# option a caller (e.g. tulpaObs) may have set for a whole fit, so progress
# reaches grids run through fixed-control internal paths (the EM per-block
# nested solve) without threading the knobs through every layer.
# `[[` (exact) not `$` on `progress`: it is a strict prefix of the other three
# keys, so a `$` read returns whichever ONE of them is set (and NULL when two
# are, the match being ambiguous). `isTRUE()` of a file path or a cadence is
# FALSE, which turned asking for the heartbeat file into switching the console
# bar off.
has_ctrl <- !is.null(control[["progress"]]) || !is.null(control$progress.every) ||
!is.null(control$progress.throttle) || !is.null(control$progress.file)
if (!has_ctrl) {
opt <- getOption("tulpa.nl_progress", NULL)
if (is.list(opt)) return(opt)
}
list(
progress = isTRUE(control[["progress"]] %||% TRUE),
progress_every = as.integer(control$progress.every %||% 0L),
progress_throttle = as.numeric(control$progress.throttle %||% 2),
progress_file = as.character(control$progress.file %||% "")
)
}
# Resolve a checkpoint spec to a normalized `list(path, resume)`.
# `control$checkpoint = list(path = , resume = TRUE)` enables it; `resume =
# FALSE` starts over (the caller removes any prior file before the first
# kernel call). Like `.nl_progress_args`, an absent control key falls back to
# the scoped `tulpa.nl_checkpoint` option (when `use_option`, the default) so
# the checkpoint reaches grids run through fixed-control internal paths;
# `tulpa_sample_glmm()`'s per-chain checkpoint is a distinct namespace from
# the nested-Laplace grid-cell one this option serves, so it reads
# `use_option = FALSE` and ignores it.
.nl_checkpoint_args <- function(control, use_option = TRUE) {
cp <- control$checkpoint
if (is.null(cp)) {
if (use_option) {
opt <- getOption("tulpa.nl_checkpoint", NULL)
if (is.list(opt)) return(opt)
}
return(list(path = "", resume = TRUE))
}
if (!is.list(cp) || is.null(cp$path) || !is.character(cp$path) ||
length(cp$path) != 1L || is.na(cp$path) || !nzchar(cp$path)) {
stop("`control$checkpoint` must be a list with a single non-empty ",
"character `path` (e.g. ",
"list(path = \"fit.ckpt\", resume = TRUE)).", call. = FALSE)
}
resume <- cp$resume %||% TRUE
if (!is.logical(resume) || length(resume) != 1L || is.na(resume)) {
stop("`control$checkpoint$resume` must be a single TRUE or FALSE.",
call. = FALSE)
}
list(path = path.expand(cp$path), resume = isTRUE(resume))
}
Try the tulpa package in your browser
Any scripts or data that you put into this service are public.
tulpa documentation built on Oct. 19, 2026, 5:09 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .nl_checkpoint_args .nl_progress_args .prior_from_spatial_spec .prior_from_temporal_spec prior_from_spec .nl_dispatch_multi .nl_grid_layout .nl_grid_crossing .nl_check_grid_cap .nl_max_grid_cells .nl_with_internal_batch .nl_internal_batch .nl_dense_grid_warn .nl_multi_grid_warn .nl_block_axis_grid .nl_block_spec_for_cpp .nl_single_block_width .nl_role_cols .nl_block_roles .is_multi_block_prior .nl_nngp_field_args .nl_block_projector .nl_refuse_projector .nl_block_is_projected .nl_attach_grid_hessians .nl_modal_joint_precision `%||%` .nl_attach_evidence .nl_outer_log_evidence .nl_normalise_weights_safe tulpa_normalise_weights_safe .nl_attach_outer_integration .nl_screen_log_offset .nl_dispatch .mcar_default_logchol_grid .nl_adj_args .nl_check_block_fields .nl_attach_pareto_k .nl_prune_gate .nl_check_screen_iters .nl_check_prune_tol tulpa_nested_laplace
Documented in prior_from_spec tulpa_nested_laplace tulpa_normalise_weights_safe
#' Nested Laplace approximation for latent Gaussian models
#'
#' @description
#' Generic outer-grid nested Laplace driver. Builds a grid over the
#' hyperparameters of a single latent prior block (spatial or temporal),
#' runs an inner Laplace at each grid point with warm-starting, and
#' integrates over the grid to give proper hyperparameter marginals.
#'
#' Supported priors:
#' * Spatial (areal): `"icar"` (1D grid on tau), `"bym2"`
#' (2D on (sigma, rho)), `"car_proper"` (2D on (tau, rho); rho lives
#' in the eigenvalue interval (1/lambda_min, 1/lambda_max) of
#' `D^{-1} W`).
#' * Spatial (continuous): `"nngp"` (2D on (sigma2, phi_gp)),
#' `"hsgp"` (2D on (sigma2, lengthscale)).
#' * Temporal: `"rw1"`, `"rw2"` (1D grid on tau), `"ar1"` (2D on (tau, rho))
#' * SPDE continuous spatial: see \code{cpp_nested_laplace_spde()} (separate
#' entry, rebuilds Q via SPDE Q-builder).
#'
#' @param y Response vector.
#' @param n_trials Trial sizes (binomial). Pass `1L`-vector otherwise.
#' @param X Fixed-effects design matrix.
#' @param prior A list describing the latent prior block. Required field
#' `type` one of \{"icar", "bym2", "car_proper", "rw1", "rw2", "ar1"\}.
#' Type-specific fields:
#' * icar: `spatial_idx`, `n_spatial_units`, `adj_row_ptr`, `adj_col_idx`,
#' `n_neighbors` (CSR adjacency, 0-based); optional `tau_grid`.
#'
#' Any block accepts an optional `svc_weight`, one weight per observation,
#' which makes it a varying coefficient: observation `i` contributes
#' `svc_weight[i] * f_i` instead of `f_i`, where `f_i` is the block's field
#' at that row (`z[spatial_idx[i]]` for an areal block, `(A u)_i` for an
#' SPDE one). Absent, the field enters unweighted.
#' * bym2: same adjacency; `scale_factor`; optional `sigma_grid`, `rho_grid`.
#' * car_proper: same adjacency; optional `tau_grid`, `rho_grid`,
#' `rho_bounds = c(lower, upper)` (defaults to (0, 1)).
#' * rw1/rw2: `temporal_idx` (1-based), `n_times`; optional `tau_grid`,
#' `cyclic` (default FALSE).
#' * ar1: `temporal_idx`, `n_times`; optional `tau_grid`, `rho_grid`
#' (each `rho_grid` value strictly inside (-1, 1)).
#'
#' An `icar`, `bym2`, `car_proper` or `nngp` block may carry `projector`, an
#' `N x n_spatial_units` (`N x n_spatial` for `nngp`) matrix (dense or a
#' \pkg{Matrix} sparse matrix), in place of `spatial_idx`: observation `i`
#' then receives `(A z)_i`, so restricted spatial regression is
#' `A = S P_perp` with `S` the observation-to-unit incidence. The block is
#' fitted by the multi-block driver, alone or in a list of blocks. An
#' intrinsic field's level is placed from `A 1`: removed with no fold where
#' it is zero, folded into the intercept where it is the intercept column,
#' otherwise left to the augmentation. The `nngp` and `car_proper` priors are
#' proper and carry their level themselves.
#'
#' Grids of a two-axis type (`bym2`, `car_proper`, `ar1`) are PAIRED when
#' both are supplied: the i-th entries of `tau_grid` / `sigma_grid` and
#' `rho_grid` form the i-th outer-grid cell, so the two must have equal
#' length, and a tensor grid is written as `expand.grid()` columns. A single
#' supplied axis is kept as an axis (its distinct values) and crossed with
#' the default of the other.
#'
#' A default grid axis is a starting axis, not a hard ceiling: for `icar`
#' (`tau_grid`) and `bym2` (`sigma_grid`) a posterior mode that rails a
#' boundary node (`pareto_k_regime = "collapsed_edge"`) triggers one
#' mode-Hessian recenter-and-refit, reported through
#' `outer_grid_placement` / `outer_grid_recenter_declined` -- see
#' [tulpa_nested_laplace_joint()]'s return docs. An axis the caller pinned is
#' never moved; mark a grid your own code defaulted with [auto_grid()] to keep
#' the recenter live on it.
#' @param re_idx Optional 1-based RE group index per obs (defaults to no RE).
#' @param n_re_groups RE group count (default 0).
#' @param sigma_re RE standard deviation (default 1).
#' @param family `"binomial"`, `"poisson"`, `"neg_binomial_2"`, etc.
#' @template phi
#' @param offset Optional per-observation offset added to the linear predictor
#' (length `length(y)`), as an `offset()` term in a model formula. Every
#' single- and multi-block kernel carries it, and `fitted_eta` includes it.
#' @param control Optional list of perf/numerical tuning knobs (statistical
#' arguments stay top-level), following the `control` convention of
#' [tulpa()]. Recognised elements (defaults in parentheses):
#' * `max_iter` (`50L`), `tol` (`1e-6`) -- inner Newton iteration budget and
#' tolerance.
#' * `n_threads` (`1L`) -- inner-loop OpenMP threads.
#' * `x_init` (`NULL`) -- warm-start for the first grid point's inner solve.
#' * `keep_grid_hessians` (`TRUE`) -- when `TRUE`, retain per-grid-point
#' fixed-effects marginal Hessian \eqn{H_\beta} and mode \eqn{\hat{\beta}}
#' on the return list as `$grid_hessians` (list of dense \eqn{p\times p}
#' matrices) and `$grid_modes` (list of length-\eqn{p} vectors). Used
#' downstream by simplified-Laplace (SLA) callers to assemble skew-aware
#' marginals -- see the cumulant pooling in [rubins_pool()].
#' * `diagnose_k` (`TRUE`), `k_samples` (`500L`) -- compute the outer
#' Pareto-\eqn{\hat{k}} accuracy diagnostic (`$pareto_k`) by importance
#' sampling the hyperparameter posterior against the Gaussian proposal
#' fitted to the grid, drawing `k_samples` extra inner-marginal evaluations.
#' `k_samples` is a precision knob: the GPD tail size is held at the
#' fraction the default budget implies, so a larger budget sharpens the
#' same k-hat instead of moving it to a deeper quantile of the weight
#' distribution (gcol33/tulpa#631). `k_tail_points` overrides that tail
#' size directly (capped at 20% of the draws).
#' Computed for a single-block, single positive-scale-axis grid; left `NA`
#' (with the grid's quadrature ESS as the fallback diagnostic) for
#' multi-block, multi-axis, or bounded-parameter grids. See [tulpa_psis()].
#' * `diagnose_skew` (`TRUE`), `skew_idx` (`NULL`) -- compute the inner-Laplace
#' skewness diagnostic (`$inner_skew`, gamma_3, Rue Martino & Chopin 2009
#' Sec 3.2.3) at the fitted MAP grid cell: one extra Newton solve, scoring
#' the `p` fixed-effects latent indices by default (pass `skew_idx`, 1-based
#' latent indices, to score additional ones, e.g. specific spatial units --
#' the full latent field is not scored by default since it costs one linear
#' solve per index). This is the complementary layer to `diagnose_k`: the
#' outer diagnostic scores the hyperparameter-grid integration around a
#' FIXED inner Laplace, this scores whether that inner Gaussian
#' approximation is itself a good fit to the latent-field conditional
#' posterior. See [diagnostics()] for the combined whole-fit verdict.
#' * `within_cell` (`"log_quadratic"`) -- the WITHIN-CELL construction the
#' reported per-axis hyperparameter intervals are read with. The outer
#' grid's weights say how much mass each cell holds; they do not say how it
#' is spread inside the cell, and a quantile needs both.
#' `"log_quadratic"` joins each row's node log densities by overlapping
#' quadratics, exact for a Gaussian at any spacing, and continues the end
#' quadratic past the outer node instead of stopping at the box edge; it
#' runs on the rows that resolve their own conditional (`h / sd <= 2`) and
#' reads every other row as boxes, declining to `"box_uniform"` where no
#' row qualifies. `"box_uniform"` puts the cumulative FULL mass at each
#' cell EDGE and interpolates between edges; `"chord"` puts the cumulative
#' MID-mass at each cell coordinate and interpolates between coordinates.
#' The default moved from `"box_uniform"` at gcol33/tulpa#932: a placed
#' axis sits at 1.25 posterior SDs per cell, where the box read is narrow
#' by construction, and against dense references of a two-arm ICAR
#' fixture the log-quadratic read's summed error is 0.038 against the box
#' read's 0.110 on the same grids. `outer_grid_h_over_sd` is how wide a
#' cell is on each axis, and `theta_within_cell` is what each axis was
#' actually read with. Only a `"density"` support admits the box and
#' log-quadratic reads -- a CCD design, a locally refined grid and a
#' posterior sample are not cell partitions that tile -- and an axis they
#' decline on reports the read that ran with a reason rather than erroring.
#' Nothing else moves: point estimates, moments, draws and weights are
#' untouched.
#' * `skew_correct` (`TRUE`) -- consume the inner-Laplace expansion instead of
#' only grading it: report
#' Cornish-Fisher marginal quantiles at each coefficient's own `gamma_3`,
#' about the centre `gamma_1 + gamma_3 / 2` that Rue, Martino & Chopin
#' (2009) eq. (22) implies, from `summary()` / `confint()` wherever the
#' combined inner band says the leading-order expansion is in its regime, and
#' the Gaussian quantiles everywhere else. It is post-processing on the
#' reported quantiles: draws, modes and weights are untouched, so a fit run
#' with it off is bit for bit the fit it was before. A coefficient whose
#' location term could not be formed declines rather than reading it as zero.
#' The band that bounded the relocation itself (`centre_unreliable`) is
#' off (`Inf`): scored over seven fixtures with an exact
#' reference, every finite cutoff declines the coefficients the correction
#' helps most, because a large centre carrying a small `gamma_3` is uniformly
#' weak correlation rather than an expansion out of its regime.
#' `$skew_correction` records the per-coefficient `gamma_3`, `gamma_1` and
#' the centre they form, the band, the inner importance k-hat, the combined
#' band, the eligibility and the reason behind it; the `skew_applied`
#' attribute on `summary()` / `confint()` records what was actually used at
#' the requested level. RMC fit a skew normal here instead; the series
#' correction is the same-order alternative, and unlike a skew normal its
#' skewness does not saturate inside the band it is applied on.
#'
#' MEASURED. Against exact quadrature
#' quantiles of rare-event binomial-logit posteriors it cuts total absolute
#' endpoint error 69.2%, improving both endpoints in every case. Scored over
#' the WHOLE marginal -- paired CRPS against the exact posterior in a
#' 400-replicate prior-predictive experiment -- it reads -0.01643 against
#' the uncorrected Laplace at t = -1.89, essentially all of the -0.01662 the
#' exact posterior itself achieves, and its PIT re-enters the simultaneous
#' SBC band. Applied about the Laplace mode instead of about
#' `gamma_1 + gamma_3 / 2` the same reshaping scored +0.00775 at t = +3.54,
#' a net loss; the location term is what supplies that centre.
#'
#' IT IS ON BY DEFAULT, so `summary()` / `confint()` on a
#' nested-Laplace fit report the corrected quantiles wherever the combined
#' inner band admits the coefficient; `skew_correct = FALSE` restores the
#' uncorrected report exactly. Scored against the mixture read a
#' correction-off fit gives, the flip is t = -1.895 on the rare-event
#' intercept and -3.765 / -3.201 on the small-group Bernoulli design. Across
#' twelve model classes read off one solve per seed, pooled 95% coverage
#' moves 0.9510 -> 0.9542 at a standard error of 0.0070, with every class
#' inside the acceptance the shipped gates use. A fit the correction cannot
#' help -- a coupled one, whose location term is unreachable -- reports what
#' it reported before, to the bit.
#' * `subspace_debias` (`FALSE`) -- correct only the latent directions the
#' inner-layer diagnostics flagged, by exact Metropolis, and leave the rest
#' at their Gaussian conditional. `TRUE` takes every default; a list
#' overrides `band`
#' (the inner-reliability floor a coordinate is selected at, default
#' `"ok"`), `idx` (pin the corrected set explicitly, skipping the selector),
#' `closure` / `closure_max` (grow the set by strongly coupled
#' precision-graph neighbours: a coordinate strongly coupled to a member of
#' the corrected set and left out of it is carried linearly, which is the
#' error the correction removes, so requesting the closure also retains the
#' modal cell's joint precision for the selector to read), the sampler
#' budget `n_iter` / `warmup` / `thin`, and
#' `n_draws`. The selector reads the per-index bands `diagnose_skew` already
#' attached, so it costs no extra solve; the correction itself re-runs the
#' settled grid once with the sampler on, and the fit then reports `$draws`
#' -- the per-cell Metropolis sample for the selected coordinates, the rest
#' from the Gaussian conditional given them -- instead of the
#' Gaussian-mixture moments. An EMPTY selection leaves the fit bit-for-bit
#' identical to the plain path. `$subspace_debias` records what was selected,
#' the bands it was read from, and the per-cell acceptance rate. Requesting
#' it turns `keep_grid_hessians` on, since the recombination reads exactly
#' those per-cell pieces.
#' * `cila` (`FALSE`) -- corrected integrated Laplace, the second inner-layer
#' debias (after Lai, Margossian and Sheldon,
#' arXiv:2605.20345). Where
#' `subspace_debias` selects coordinates and runs exact Metropolis on them,
#' this selects nothing: at every outer cell it draws `n_points` points from
#' the whole inner Gaussian, weights each by the exact joint density it came
#' from, and reports the weighted particles. `TRUE` takes the defaults; a
#' list overrides `n_points` (`1024L`), `variant` (`"qmc"`, a Sobol net;
#' `"is"` for iid draws, `"rqmc"` for the net under `n_shift` random
#' shifts), `n_shift` (`8L`), `n_draws` and `seed`. Below 512 points a
#' cell's particle set is too coarse to be a marginal at all and the request
#' is refused. The corrected per-cell masses become the fit's own
#' `weights` / `log_marginal` and `weights_source` reports `"cila"`; the
#' pre-correction pair is kept as `$cila$laplace`. A cell whose inner solve
#' factorized sparsely draws through the CHOLMOD factor's own triangular and
#' permutation solves; an LDL' factor has no square root to draw with and is
#' declined with `"sparse_factor_not_ll"`.
#' * `auto_recenter` (`TRUE`) -- outer-grid placement policy. `TRUE`
#' re-centres the movable default
#' axes on the posterior mode and refits when the grid either RAILS (an
#' axis's own marginal is maximal at one of its own endpoints) or does not
#' RESOLVE its posterior (an axis's node spacing exceeds 2 posterior SDs in
#' its own coordinate). Both tests read the weights the fit already stored,
#' so a grid that brackets and resolves its mode costs nothing beyond them;
#' when the pass does fire it is a second full grid solve plus a
#' finite-difference mode/Hessian stencil. A recentred axis is
#' `mode +/- 2.5 sd` over 5 nodes, a cell width of 1.25 posterior SDs by
#' construction, against a census median of 3.9 on the fixed spans.
#'
#' Measured over 200 fixed-truth seeds on each of six configurations (icar
#' chain / icar lattice / rw1 / bym2 / iid / nngp), the default moves mean
#' |coverage - nominal| from 0.043 to 0.030 at the 95% level, 0.171 to
#' 0.084 at 80% and 0.243 to 0.129 at 50% against the rail-only policy, at
#' 0.63 times the 95% interval width and 0.76 times the median bias, for
#' 1.71 times the wall clock.
#'
#' Three other values. `"rail"` is the rail test alone, which is what
#' `TRUE` meant before the sizing measurement settled the default. `FALSE`
#' integrates over the grid exactly as given, whatever it is, and records
#' `outer_grid_recenter_declined = "auto_recenter_disabled"`. `"always"`
#' re-centres every movable default axis whatever the fit did, at 2.04
#' times the wall clock; it agrees with the default seed for seed on five
#' of the six measured configurations and differs on the one whose default
#' axes already resolve their posterior, where it takes 50% coverage to
#' 0.135 against a nominal 0.5.
#'
#' Which families carry movable axes is `.NL_REGISTRY_AXIS_FIELD`: icar,
#' rw1, rw2, iid, bym2, nngp, hsgp and spde. car_proper, ar1 and hsgp_mo
#' each carry a correlation axis with no guessable coordinate, and mcar /
#' miid / tgmrf hold their axes in one matrix field; a fit of any of them
#' records which through `outer_grid_recenter_declined` rather than passing
#' in silence. The per-axis policy names are the standalone registry path
#' only -- [tulpa_nested_laplace_joint()] and [fit_st_nested()] refuse them
#' rather than accept them and ignore them.
#' * `max_grid_cells` (`NULL`, no ceiling) -- a ceiling on a multi-block
#' outer grid's cell count, one inner Newton solve per cell; a grid past it
#' is refused before any solve, naming the axes that produced the count.
#' Unset, a grid of any size is solved, and one past 2048 cells is announced
#' with a warning first.
#' * `checkpoint` (`list(path =, resume =)`) -- grid-cell checkpoint /
#' resume. Each solved outer cell is appended to `path`, keyed by its
#' hyperparameter coordinate; `resume = TRUE` loads the finished cells and
#' solves only the rest, and a fingerprint mismatch (different data,
#' settings or grid) errors rather than resuming onto a stale result. See
#' `?tulpa_nested_laplace` "Checkpoint / resume".
#' * `progress`, `progress.every`, `progress.file`, `progress.throttle` --
#' the outer-grid progress reporter: whether to print, how often (in
#' cells), where to, and the minimum seconds between lines.
#' * `prune` (`TRUE`), `prune_tol` (`1e-3`), `screen_iters`
#' (`.NL_SCREEN$iters`, 2 -- the depth gcol33/tulpa#640 measured the screen
#' down to) -- cheap-pass screening of the outer grid. The driver first
#' sweeps the lattice running a `screen_iters`-step inner Newton per cell,
#' each warm-started from the previous screened cell's quasi-mode, computes
#' a screening Laplace log-marginal, softmax-normalises it, and skips the
#' full inner Newton on every cell whose screened weight is below
#' `prune_tol`. The neighbour warm-start keeps each cheap mode near its
#' cell's own mode, so the cheap ranking tracks the full-solve ranking even
#' where the latent mode moves across the grid. Pruned cells get
#' `log_marginal = -Inf`, `n_iter = 0` and inherit the pilot mode; the pilot
#' cell is never pruned, and at least `prune_min_keep` cells (5) are solved
#' in full whatever the tolerance says -- the highest-ranked dropped cells
#' are restored up to that floor, because the outer grid placement pass
#' reads a finite-difference curvature off the cells that were SOLVED and a
#' kept set of one carries none. After the full pass the posterior mass the
#' dropped cells may carry is bounded by lifting each by the worst screening
#' error measured on a kept cell (`prune_dropped_mass_bound`); past 1% the
#' driver solves the dropped cells that contribute most to it until it is
#' under 1% (`prune_n_repaired`), and a safety gate falls back to the full
#' grid (with a warning) if the bound still fails. The bound assumes a
#' dropped cell is screened no worse than the worst kept one. On by default
#' (gcol33/tulpa#925); `FALSE` solves every cell.
#' `prune_tol` must be a finite numeric in `[0, 1)` and `screen_iters` a
#' single integer `>= 1`; both are validated whatever `prune` says.
#' `prune` / `prune_tol` reach both the single-block and the multi-block
#' dispatch; the screening DEPTH is declared by the single-block kernels
#' only, so a multi-block prior refuses a pinned `screen_iters` rather than
#' accepting it and ignoring it.
#' * `prune_log_gap` (unset) -- the screening cut stated in nats instead of as
#' a normalised weight: keep every cell within this many nats of the best
#' screened cell. `prune_tol` is a weight, so what it cuts at is the gap
#' `-(log(prune_tol) + log(Z))`, and the whole range a caller would type
#' (`1e-3` to `1e-12`) is 7 to 28 nats of it -- on an outer surface whose
#' log-marginal spans thousands of nats every setting returns the same kept
#' set. Stating the gap directly keeps the knob's resolution wherever the
#' surface is steep. It reaches the kernel as `prune_tol = exp(-gap)`, so
#' the realised cut is the requested gap less `log(Z)`, at most `log(n_grid)`
#' nats narrower; the fit reports the realised value. Setting both this and
#' `prune_tol` is an error -- they state one cut in two units.
#' * `fitted_var` (`TRUE`) -- fill the per-row predictive variance
#' `fitted_eta_var`, the companion of `fitted_eta` a caller marginalises the
#' per-row linear predictor with. It costs one back-solve per distinct
#' loading vector per cell, which on a design with few repeated rows is the
#' dominant cost of a cell, so a caller reading only `fitted_eta` or the
#' coefficient summaries sets `FALSE` and the fit then carries no
#' `fitted_eta_var`. Every single-block field reports it, the NNGP, HSGP,
#' SPDE and spatiotemporal fields included, each read off its cell's own
#' precision at the mode; a multi-block prior refuses `FALSE` rather than
#' accepting it and ignoring it.
#' * `var_of_means_consistency` (`TRUE`) -- the consistency pass of
#' [tulpa_nested_laplace_joint()]: a field SD (`tau` of an `icar`, `rw1` or
#' `rw2` block, `sigma` of an `iid`, `bym2` or `spde` block) whose marginal
#' has collapsed onto too few nodes to carry a spread has levels added in
#' every row of the other axes that holds the posterior until it resolves, and the new cells are
#' solved and carried like the rest (`var_of_means_consistency_info`
#' records the pass). Other axes keep the nodes they were laid on.
#' * `adaptive_grid` (`FALSE`), `adaptive_grid_edge_thresh` (`0.02`),
#' `adaptive_grid_max_passes` (`1L`) -- the opt-in boundary / interior pass,
#' on the same axes.
#' * `axis_refine` (`NULL`) -- per-axis override of how far those axes may be
#' moved, named by the grid's axes (`"tau"`; `"b1.tau"` on a multi-block
#' prior): `"none"`, `"densify"` (inside the declared span) or `"extend"`.
#' Nodes the caller wrote down densify by default and an engine-placed axis
#' extends, as in [tulpa_nested_laplace_joint()].
#'
#' @return A list with:
#' * `theta_grid`: matrix or vector of grid hyperparameter values.
#' * `log_marginal`: log p(y, mode | theta_k) at each grid point.
#' * `weights`: integration weights normalising to sum 1.
#' * `theta_mean`, `theta_sd`: posterior moments per hyperparameter.
#' * `n_iter`: inner Newton iterations per grid point.
#' * `modes`: matrix `[n_grid x n_x]` of inner modes, when stored.
#' * `pareto_k`, `pareto_k_is_ess`: outer Pareto-\eqn{\hat{k}} and its
#' importance-sampling ESS (`NA` when not computed for the grid; see
#' `control$diagnose_k`).
#' * `pareto_k_proposal_source`, `pareto_k_first_pass`: which proposal family
#' the reported k-hat came from (the outer k-hat is scored against several
#' candidates and the best is kept), and the k-hat of the FIRST pass -- the
#' proposal exactly as this backend placed it, before any candidate refined
#' or replaced it. A large gap between the two says the nodes are badly
#' scaled around the hyperparameter posterior even though the verdict is
#' fine; no gap says the placement was already right.
#' * `inner_skew`, `inner_skew_idx`, `inner_skew_dropped`: the inner-Laplace
#' skewness diagnostic (gamma_3) at each scored latent index and its
#' 1-based index, plus a count of (index, observation) contributions
#' dropped for a non-finite third derivative (see `control$diagnose_skew`).
#' * `timing`: named numeric of wall-clock seconds (`total`, `setup`,
#' `grid`, `postproc`, `diagnostics`); the `grid` phase is the inner
#' Laplace pass that scales with grid size. Surfaced one-line in `print`.
#' * `prune_cheap_log_marginal`, `prune_mask`, `prune_n_pruned`, `prune_tol`:
#' present when the grid was screened and the safety gate did not
#' trip -- the cheap-screen log-marginal per cell (the truncated solve's
#' value plus half its Newton decrement, the second-order estimate of the
#' converged value the cells are ranked by; the decrement itself is
#' `prune_screen_decrement`), the logical mask of
#' pruned cells, the pruned-cell count, and the threshold actually applied.
#' * `prune_log_gap_cut`, `prune_cheap_lm_spread`, `prune_min_keep`,
#' `prune_n_floor_restored`: the same screen read on the scale it operates
#' on -- how many nats below the best screened cell the tolerance cut at,
#' how many nats the screened surface spans, the kept-cell floor applied,
#' and how many cells that floor put back. A cut that is a sliver of the
#' spread is a tolerance with no resolution on this grid.
#' * `prune_dropped_mass_bound`, `prune_n_repaired`, `prune_gate_mass`: the
#' posterior mass the dropped cells may carry, with each lifted by the worst
#' screening error measured on a kept cell; how many dropped cells the
#' driver solved in full to bring that under the gate; and the gate.
#' * `prune_fallback_triggered`, `prune_fallback_reason`: present only when
#' the gate did trip; the rest of the fit is then the full-grid result and
#' the reason string records which check fired.
#' * `prior`: echoed input.
#'
#' @param spec Optional `tulpa_temporal` or `tulpa_spatial` spec object
#' (output of [temporal_rw1()], [temporal_rw2()], [temporal_ar1()],
#' [spatial_car()], [spatial_bym2()], etc.). When supplied alongside
#' `data`, the `prior` list is built automatically via
#' [prior_from_spec()] -- pass either `prior` or `spec`, not both.
#' @param data Data frame used to validate `spec` and resolve
#' time/group/site indices. Required when `spec` is supplied.
#' @param likelihood Optional model-supplied likelihood, replacing the built-in
#' `family`. Pass an external pointer to a `tulpa::NestedLikelihood` (built in
#' a model package's own C++ from a `LikelihoodSpec`); the inner Laplace solve
#' then reads the per-observation score, Fisher weight, and log-likelihood from
#' that spec instead of `family`, so `family`/`phi` are ignored. Used by model
#' packages to fit a custom response without adding a family to tulpa -- for
#' example tulpaObs threads its marginalized single-season occupancy
#' likelihood (a scaled Bernoulli, with the latent occupancy state integrated
#' out) through this. Multi-block `prior` only. Default `NULL` (use `family`).
#' @template hyperprior
#'
#' @references
#' Rue, Martino & Chopin (2009). Approximate Bayesian inference for latent
#' Gaussian models by using integrated nested Laplace approximations.
#' \emph{JRSS-B} 71(2):319-392.
#' @examples
#' \donttest{
#' set.seed(1)
#' S <- 30L # spatial units arranged in a chain
#' nb <- lapply(seq_len(S), function(s) setdiff(c(s - 1L, s + 1L), c(0L, S + 1L)))
#' nn <- lengths(nb)
#' field <- as.numeric(scale(cumsum(rnorm(S, 0, 0.4)))) # smooth spatial field
#' idx <- rep(seq_len(S), each = 6L); n <- length(idx); x <- rnorm(n)
#' y <- rbinom(n, 1L, plogis(-0.2 + 0.6 * x + field[idx]))
#' prior <- list(type = "icar", n_spatial_units = S, spatial_idx = idx,
#' adj_row_ptr = c(0L, cumsum(nn)), adj_col_idx = unlist(nb) - 1L,
#' n_neighbors = nn, tau_grid = c(0.5, 1, 2, 4, 8))
#' fit <- tulpa_nested_laplace(y, rep(1L, n), cbind(1, x), prior = prior,
#' family = "binomial")
#' fit$theta_mean # marginalized ICAR precision
#' }
#' @export
tulpa_nested_laplace <- function(y, n_trials, X, prior = NULL,
spec = NULL, data = NULL,
re_idx = NULL, n_re_groups = 0L, sigma_re = 1.0,
family = "binomial", phi = 1.0, offset = NULL,
likelihood = NULL,
hyperprior = c("proper", "flat"),
control = list()) {
tulpa_check_control(control, .CONTROL_KEYS$nested_laplace, "tulpa_nested_laplace")
# A typo is named here, against the family registry, rather than by the C++
# family lookup, whose message names a header (gcol33/tulpa#886). A consumer
# `likelihood` supplies its own density and bypasses the builtin family.
if (is.null(likelihood)) .family_or_stop(.canonical_family(family))
hyperprior <- .hp_choice(match.arg(hyperprior))
tm <- .tulpa_timer()
# Perf/numerical knobs live in `control = list()` (matching tulpa()); the
# top-level signature carries only statistical arguments.
max_iter <- control$max_iter %||% 50L
tol <- control$tol %||% 1e-6
n_threads <- control$n_threads %||% 1L
x_init <- control$x_init
# Cheap-pass outer-grid screening. The kernels read a positive tolerance as
# "screen the lattice and prune", so the user-facing `prune` toggle gates it:
# FALSE always solves every cell, TRUE screens at `prune_tol`. Both the
# tolerance and the screening depth are validated whatever `prune` says, so a
# misspecified value is reported where the user set it rather than at the next
# fit that happens to switch pruning on.
prune <- isTRUE(control$prune %||% .nl_screen("prune")[["registry"]])
# `prune_log_gap` states the same cut in nats -- the units the screened
# surface is measured in -- and replaces the tolerance when set.
prune_tol <- .nl_prune_tol_from_control(
control,
.nl_check_prune_tol(control$prune_tol %||% .nl_screen("prune_tol")))
prune_tol_eff <- if (prune) prune_tol else 0
screen_iters <- .nl_check_screen_iters(
control$screen_iters %||% .nl_screen("iters"))
fitted_var <- isTRUE(control$fitted_var %||% TRUE)
# Subspace debias. The correction reports per-cell DRAWS
# instead of the Gaussian-mixture moments, and those are built from the same
# per-cell fixed-effect mode + precision `.nested_fixed_moments()` reads, so
# requesting it turns that retention on: without it there is nothing to
# recombine the sampled coordinates with.
sd_cfg <- .subspace_debias_config(control$subspace_debias)
# Corrected integrated Laplace, the second inner debias.
# Same shape as the subspace one: a second pass over the settled grid whose
# per-cell particles are pooled into the reported draws.
cila_cfg <- .cila_config(control$cila)
keep_grid_hessians <- isTRUE(control$keep_grid_hessians %||% TRUE) || !is.null(sd_cfg)
# Retaining the modal cell's joint precision is what makes the debias closure
# available on a grid fit; nothing else reads it, so it rides the request.
keep_joint_prec <- !is.null(sd_cfg) && !identical(sd_cfg$closure, FALSE)
diagnose_k <- isTRUE(control$diagnose_k %||% TRUE)
k_samples <- as.integer(control$k_samples %||% .nl_diag("k_samples"))
k_tail_points <- control$k_tail_points
diagnose_skew <- isTRUE(control$diagnose_skew %||% TRUE)
skew_idx <- control$skew_idx
skew_correct <- isTRUE(control$skew_correct %||% .nl_diag("skew_correct"))
# Which WITHIN-CELL construction the reported per-axis hyperparameter
# intervals are read off. A numerical reconstruction of the
# integration design the fit already ran, so it sits beside `integration` /
# `local_ccd` / `skew_correct` in `control` rather than in the signature,
# which carries only statistical arguments. The default is
# `.nl_diag("within_cell")`, and `"chord"` restores the pre-0.0.188 report
# byte-identically.
within_cell <- .nl_within_cell_mode(control$within_cell)
# Outer-grid refinement: the passes the joint drivers run
# (`.joint_refine_outer_grid()`), over the field SD axes this door can write a
# cell onto. The consistency pass is on by default, the boundary / interior
# pass opt in, as there.
var_of_means_consistency <- isTRUE(control$var_of_means_consistency %||% TRUE)
adaptive_grid <- isTRUE(control$adaptive_grid %||% FALSE)
adaptive_grid_edge_thresh <- control$adaptive_grid_edge_thresh %||% 0.02
adaptive_grid_max_passes <- control$adaptive_grid_max_passes %||% 1L
# Grid-cell checkpoint/resume. `control$checkpoint =
# list(path =, resume =)` makes every grid cell append to `path`; a resume
# loads the finished cells and solves only the rest. A fresh (resume = FALSE)
# run removes any stale file once here, before the first kernel call, so the
# several within-fit kernel calls (initial grid, adaptive refinement) all
# append rather than truncate each other. The path threads into `cargs` and
# reaches both the single-block (.nl_dispatch) and multi-block
# (.nl_dispatch_multi) kernels by name; the k-hat diagnostic re-evaluations
# are run with it stripped so they do not pollute the file.
.ckpt <- .nl_checkpoint_args(control)
if (nzchar(.ckpt$path) && !isTRUE(.ckpt$resume) && file.exists(.ckpt$path)) {
file.remove(.ckpt$path)
}
# The caller's outer-grid ceiling, published for the duration of this fit so
# the initial dispatch and every re-dispatch (the k-hat re-evaluations, the
# subspace-debias refit) read the same value.
.op_grid_cap <- options(tulpa.nl_max_grid_cells = .nl_max_grid_cells(control))
on.exit(options(.op_grid_cap), add = TRUE)
if (!is.null(spec)) {
if (!is.null(prior)) {
stop("Pass either `spec` or `prior`, not both.", call. = FALSE)
}
if (is.null(data)) {
stop("`data` is required when `spec` is supplied.", call. = FALSE)
}
prior <- prior_from_spec(spec, data)
}
# Allow passing a single tgmrf S3 object directly: wrap it in a length-1
# block list so the multi-block dispatch picks it up. This keeps the
# ergonomic API close to the formula sketch (one block,
# one call) while still routing through the validated multi-block path.
if (inherits(prior, "tulpa_latent_block")) {
prior <- list(prior)
}
if (!is.list(prior)) {
stop("`prior` must be a list (single block) or list-of-lists (multi-block).",
call. = FALSE)
}
# This door's outer grid is spanned by the latent blocks' hyperparameters
# alone, so a prior with no block leaves nothing to integrate. The joint door
# integrates per-arm dispersion axes as well and does take an empty prior.
if (!length(prior)) {
stop("`prior` holds no latent block, so tulpa_nested_laplace() has no ",
"hyperparameter to integrate. Fit a fixed-effects model with ",
"tulpa(mode = \"laplace\"), or integrate a dispersion axis with ",
"tulpa_nested_laplace_joint(prior = list(), phi_grid = ).",
call. = FALSE)
}
# A projected areal field is read by the multi-block driver alone, so a
# single projected block routes there as a length-1 block list.
if (.nl_block_is_projected(prior)) prior <- list(prior)
# Outer-axis provenance: record which grid axes the caller
# declared as defaults with `auto_grid()` -- the registry rescue below reads
# it -- and strip the markers so every downstream consumer sees plain grids.
.prov <- .nl_grid_provenance(prior)
prior <- .prov$prior
# And refuse an axis this path cannot read: a grid field
# some OTHER driver spells the same family with is neither integrated nor
# reported here, so a pinned one is an error and a defaulted one is recorded.
.op_axis <- .nl_publish_axis_dropped(
.nl_check_axis_fields(prior, "registry", auto = .prov$auto))
on.exit(options(.op_axis), add = TRUE)
# A model-supplied `likelihood` is only honoured by the spec-driven multi-block
# driver (run_multi_block_nested_laplace); the single-kernel path (.nl_dispatch)
# has no spec hook. Route a single LatentBlock-type prior through the
# multi-block driver by wrapping it as a length-1 block list -- the per-cell
# math is identical (every single-block kernel already calls the same driver).
if (!is.null(likelihood) && !is.null(prior$type)) {
prior <- list(prior)
}
vd <- .validate_glm_design(y, X, n_trials, "tulpa_nested_laplace")
N <- vd$N
n_trials <- vd$n_trials
if (is.null(re_idx)) re_idx <- rep(0L, N)
re_idx <- .validate_re_idx(re_idx, n_re_groups, N, "tulpa_nested_laplace")
if (!is.null(offset)) {
offset <- as.numeric(offset)
if (length(offset) != length(y) || any(!is.finite(offset))) {
stop("`offset` must be a finite numeric vector of length length(y) (",
length(y), ").", call. = FALSE)
}
}
cargs <- list(
y = as.numeric(y),
n = as.integer(n_trials),
offset_nullable = offset,
X = X,
re_idx = as.numeric(re_idx),
n_re_groups = as.integer(n_re_groups),
sigma_re = as.numeric(sigma_re),
family = family,
phi = .phi_to_kernel(family, as.numeric(phi)),
max_iter = as.integer(max_iter),
tol = as.numeric(tol),
n_threads = as.integer(n_threads),
x_init_nullable = x_init,
store_Q = isTRUE(keep_grid_hessians),
checkpoint_path = .ckpt$path,
prune_tol = prune_tol_eff,
screen_iters = screen_iters,
compute_fitted_var = fitted_var,
# Read by the dispatchers, which fold the outer hyperprior in R; never
# handed to a kernel.
hyperprior = hyperprior
)
# cargs without the checkpoint, for the k-hat diagnostic re-evaluations, and
# without the prune: a screened cell carries `log_marginal = -Inf`, which is
# not the inner marginal the importance ratio and the placement stencil read.
# They evaluate the grid they are handed, so every cell of it has to be solved.
cargs_no_ckpt <- utils::modifyList(
cargs, list(checkpoint_path = "", prune_tol = 0))
p_fixed <- ncol(X)
if (.is_multi_block_prior(prior)) {
# The multi-block entry takes the prune tolerance but declares neither the
# screening depth nor the per-row variance switch, so a pinned one is
# refused here rather than accepted and ignored.
if (!is.null(control$screen_iters)) {
stop("`control$screen_iters` applies to a single-block `prior`; ",
"this fit takes the multi-block dispatch.", call. = FALSE)
}
if (!is.null(control$fitted_var) && !isTRUE(control$fitted_var)) {
stop("`control$fitted_var` applies to a single-block `prior`; ",
"this fit takes the multi-block dispatch.", call. = FALSE)
}
tm$mark("setup")
# Every full-grid solve of this fit -- the first one and the placement
# pass's refit -- is screened at the same tolerance and gated the same way,
# so no reported grid is a pruned one whose screen ranking was unreliable.
solve_grid_multi <- function(prior_i) {
out <- .nl_prune_gate(
.nl_dispatch_multi(cargs, prior_i, likelihood = likelihood,
progress = .nl_progress_args(control),
within_cell = within_cell),
prune_tol_eff,
function() .nl_dispatch_multi(
utils::modifyList(cargs, list(prune_tol = 0)), prior_i,
likelihood = likelihood, progress = .nl_progress_args(control),
within_cell = within_cell))
refined <- .nl_refine_registry(
out, function(theta_mat, screen = NULL) .nl_dispatch_multi(
.nl_refine_cargs(cargs_no_ckpt, screen, prune_tol_eff), prior_i,
likelihood = likelihood, theta_grid_override = theta_mat),
prior_i, TRUE, .prov$auto,
consistency = var_of_means_consistency,
adaptive_grid = adaptive_grid,
edge_thresh = adaptive_grid_edge_thresh,
max_passes = adaptive_grid_max_passes,
user = control$axis_refine)
if (is.null(refined$refining_axis)) return(refined)
# The dispatcher integrated the base grid; the appended cells are measured
# and read with the rest.
refined <- .nl_attach_outer_integration(refined, "multi-block outer grid")
.nl_posterior_moments_multi(refined, refined$blocks, refined$axis_offsets,
refined$theta_grid, within = within_cell)
}
res <- solve_grid_multi(prior)
tm$mark("grid")
if (isTRUE(keep_grid_hessians)) {
res <- .nl_attach_grid_hessians(res, p_fixed, keep_joint_prec)
}
tm$mark("postproc")
res <- .nl_attach_pareto_k(res, prior, cargs_no_ckpt, "multi", NULL,
likelihood, k_samples, compute = diagnose_k,
tail_points = k_tail_points)
# The same placement pass the single-block path runs, over the block list.
# Every axis is block-prefixed here, so a moved block is re-crossed on its
# own fields and its neighbours keep theirs; a prior whose blocks the
# family table does not cover declines with a reason rather than leaving
# the fit with no placement record at all.
multi_rescue <- .nl_registry_grid_rescue(
res, NULL, prior,
refit = function(prior_i) {
r <- solve_grid_multi(prior_i)
if (isTRUE(keep_grid_hessians)) r <- .nl_attach_grid_hessians(r, p_fixed, keep_joint_prec)
.nl_attach_pareto_k(r, prior_i, cargs_no_ckpt, "multi", NULL,
likelihood, k_samples, compute = diagnose_k,
tail_points = k_tail_points)
},
refit_log_marginal = function(prior_i, theta_mat) {
.nl_dispatch_multi(cargs_no_ckpt, prior_i, likelihood = likelihood,
theta_grid_override = theta_mat)$log_marginal
},
auto = .prov$auto, multi = TRUE,
policy = .nl_recenter_mode(control$auto_recenter))
res <- multi_rescue$res
prior <- multi_rescue$prior
res$outer_grid_placement <- res$outer_grid_placement %||% "fixed"
res <- .nl_inner_skew_at_theta(res, prior, cargs_no_ckpt, "multi", NULL,
likelihood, p_fixed, skew_idx,
compute = diagnose_skew)
res <- .nl_skew_correction_attach(res, p_fixed, skew_correct)
res <- .nl_subspace_debias_attach(
res, sd_cfg,
redispatch = function(req) .nl_dispatch_multi(
cargs_no_ckpt, prior, likelihood = likelihood, debias = req,
theta_grid_override = .nl_refined_theta(res)),
p_fixed = p_fixed, beta_names = colnames(X))
res <- .nl_cila_attach(
res, cila_cfg,
redispatch = function(req) .nl_dispatch_multi(
cargs_no_ckpt, prior, likelihood = likelihood, cila = req,
theta_grid_override = .nl_refined_theta(res)),
p_fixed = p_fixed, beta_names = colnames(X),
remoments = function(r) .nl_posterior_moments_multi(
r, r$blocks, r$axis_offsets, r$theta_grid, within = within_cell))
tm$mark("diagnostics")
res$prior <- prior
res$timing <- tm$timing()
return(.finalize_fit(res, backend = "nested_laplace",
n_fixed = p_fixed, fixed_names = colnames(X),
extra_class = c("tulpa_nested_laplace", "list"),
data = list(y = y, n_trials = n_trials,
model_matrix = X, family = family,
offset = offset, phi = phi)))
}
if (!is.null(likelihood)) {
stop("`likelihood` (model-supplied spec) needs a latent prior block. ",
"Pass `prior` (a block or list of blocks).", call. = FALSE)
}
if (is.null(prior$type)) {
stop("Single-block `prior` must have a `type` field, ",
"or pass a list of blocks for multi-block.", call. = FALSE)
}
type <- tolower(prior$type)
tm$mark("setup")
# Every full-grid solve of this fit -- the first one and the placement pass's
# refit -- is screened at the same tolerance and gated the same way, so no
# reported grid is a pruned one whose screen ranking was unreliable.
solve_grid <- function(prior_i) {
out <- .nl_prune_gate(
.nl_dispatch(type, cargs, prior_i, prior_i),
prune_tol_eff,
function() .nl_dispatch(type, utils::modifyList(cargs, list(prune_tol = 0)),
prior_i, prior_i))
.nl_refine_registry(
out, function(theta_mat, screen = NULL) {
blk2 <- .nl_registry_write_theta(
list(prior_i), theta_mat, colnames(theta_mat))[[1L]]
.nl_dispatch(type, .nl_refine_cargs(cargs_no_ckpt, screen, prune_tol_eff),
blk2, prior_i)
},
prior_i, FALSE, .prov$auto,
consistency = var_of_means_consistency,
adaptive_grid = adaptive_grid,
edge_thresh = adaptive_grid_edge_thresh,
max_passes = adaptive_grid_max_passes,
user = control$axis_refine)
}
res <- solve_grid(prior)
tm$mark("grid")
# Integrate exp(log_marginal) over the outer grid. `log_marginal` is the
# inner marginal in the internal (log-scale) hyperparameter parameterization
# the grid is laid out on, and each cell carries the prior mass of the node it
# sits at: its cell width on that coordinate times the axis's declared density
# (R/hyper_quadrature.R). On an evenly spaced grid the widths are equal and
# this reduces to the plain marginal softmax, so where the nodes fall changes
# the quadrature error and not the measure being integrated.
res <- .nl_attach_outer_integration(res, "outer grid")
res <- .nl_posterior_moments(res, type, within = within_cell)
if (isTRUE(keep_grid_hessians)) {
res <- .nl_attach_grid_hessians(
res, p_fixed, keep_joint_prec,
var_cols = .nl_role_cols(prior, NULL, res$Q_csc_n, "temporal"))
}
tm$mark("postproc")
res <- .nl_attach_pareto_k(res, prior, cargs_no_ckpt, "single", type, NULL,
k_samples, compute = diagnose_k,
tail_points = k_tail_points)
# Auto-recenter a railed default outer axis (the registry
# generalization of the joint-path fix, extended over the family table):
# a single mode-Hessian recenter-and-refit when the fixed default axis
# does not contain its own mode. A no-op (zero extra fit) unless it actually
# did; a PINNED axis always wins (`.nl_axis_is_pinned()`: a grid the caller
# named and neither marked with `auto_grid()` nor equal to the engine's own
# default axis). Which families are covered, and why each of the rest
# declines, is `.NL_REGISTRY_AXIS_FIELD`.
registry_rescue <- .nl_registry_grid_rescue(
res, type, prior,
refit = function(prior_i) {
r <- solve_grid(prior_i)
r <- .nl_attach_outer_integration(r, "outer grid")
r <- .nl_posterior_moments(r, type, within = within_cell)
if (isTRUE(keep_grid_hessians)) {
r <- .nl_attach_grid_hessians(
r, p_fixed, keep_joint_prec,
var_cols = .nl_role_cols(prior_i, NULL, r$Q_csc_n, "temporal"))
}
.nl_attach_pareto_k(r, prior_i, cargs_no_ckpt, "single", type, NULL,
k_samples, compute = diagnose_k,
tail_points = k_tail_points)
},
refit_log_marginal = function(prior_i, theta_mat) {
blk2 <- .nl_registry_write_theta(
list(prior_i), theta_mat, colnames(theta_mat))[[1L]]
.nl_dispatch(type, cargs_no_ckpt, blk2, prior_i)$log_marginal
},
auto = .prov$auto,
policy = .nl_recenter_mode(control$auto_recenter))
res <- registry_rescue$res
prior <- registry_rescue$prior
res$outer_grid_placement <- res$outer_grid_placement %||% "fixed"
res <- .nl_inner_skew_at_theta(res, prior, cargs_no_ckpt, "single", type,
NULL, p_fixed, skew_idx,
compute = diagnose_skew)
res <- .nl_skew_correction_attach(res, p_fixed, skew_correct)
res <- .nl_subspace_debias_attach(
res, sd_cfg,
redispatch = function(req) .nl_dispatch(
type, utils::modifyList(cargs_no_ckpt, list(debias = req)),
.nl_refined_block(res, prior), prior),
p_fixed = p_fixed, beta_names = colnames(X))
res <- .nl_cila_attach(
res, cila_cfg,
redispatch = function(req) .nl_dispatch(
type, utils::modifyList(cargs_no_ckpt, list(cila = req)),
.nl_refined_block(res, prior), prior),
p_fixed = p_fixed, beta_names = colnames(X),
remoments = function(r) .nl_posterior_moments(r, type, within = within_cell))
tm$mark("diagnostics")
res$prior <- prior
res$timing <- tm$timing()
.finalize_fit(res, backend = "nested_laplace",
n_fixed = p_fixed, fixed_names = colnames(X),
extra_class = c("tulpa_nested_laplace", "list"),
data = list(y = y, n_trials = n_trials, model_matrix = X,
family = family, offset = offset, phi = phi))
}
# --- cheap-pass screening knobs ---------------------------------------------
#
# The kernels read a tolerance in [0, 1) and a screening depth of at least one
# Newton step. A value outside either range is a caller error rather than
# something to clamp: a clamped tolerance reads back off the fit as a setting
# that was applied.
.nl_check_prune_tol <- function(x) {
if (length(x) != 1L || !is.numeric(x) || !is.finite(x) || x < 0 || x >= 1) {
stop("`prune_tol` must be a finite numeric in [0, 1).", call. = FALSE)
}
as.numeric(x)
}
.nl_check_screen_iters <- function(x) {
if (length(x) != 1L || !is.numeric(x) || !is.finite(x) || x < 1 ||
x > .Machine$integer.max || x != round(x)) {
stop("`screen_iters` must be a single integer >= 1.", call. = FALSE)
}
as.integer(x)
}
# Gate one screened grid solve. At a zero tolerance nothing was screened and
# the kernel result is returned untouched, so a fit with pruning off takes the
# same path it always did. Otherwise the shared cheap-pass gate decides whether
# the kept set can be trusted, re-running `resolve_full` (the same solve at
# `prune_tol = 0`) when it cannot.
.nl_prune_gate <- function(res, prune_tol, resolve_full) {
if (prune_tol <= 0) return(res)
.joint_prune_safety_gate(res, resolve_full = resolve_full,
fn = "tulpa_nested_laplace()")
}
# Positive-scale hyperparameter grid fields: one log-transform unconstrains
# them, so the outer Pareto-k-hat (`.nested_grid_pareto_k`) is well-defined and
# needs no bounded-parameter Jacobian. A block carrying a bounded axis (e.g. a
# correlation `rho_grid`) is absent from this set and is declined to the
# quadrature-ESS fallback rather than transformed with a guessed support.
.NL_POS_GRID <- c("sigma_grid", "sigma2_grid", "tau_grid", "phi_gp_grid",
"phi_grid", "range_grid", "scale_grid", "sigma_spatial_grid",
"lambda_grid", "kappa_grid")
# Attach res$pareto_k for a single-block, single positive-scale-axis nested fit.
# Leaves res$pareto_k = NA (the diagnostic layer then reports the grid's
# quadrature ESS) for multi-block, multi-axis, or bounded-parameter grids.
# `refit` re-evaluates the inner marginal at a substituted grid through the SAME
# driver `res` came from; run with the RNG restored so the fit is unperturbed.
.nl_attach_pareto_k <- function(res, prior, cargs, dispatch_kind, type,
likelihood, n_samples, compute = TRUE,
tail_points = NULL) {
res$pareto_k <- NA_real_
res$pareto_k_is_ess <- NA_real_
res$pareto_k_scope <- "outer (hyperparameter) Gaussian proposal"
# The outer-integration regime is read off the stored grid
# weights, so it costs nothing and is attached on every path -- including the
# ones that DECLINE the k-hat below, where knowing the grid collapsed (and
# whether against a boundary) is the only outer signal left.
res <- .joint_attach_pareto_k_regime(res)
# Every decline below says which one it was.
decline <- function(reason, detail = NULL) {
.k_attach_declined(res, .k_decline(reason, detail))
}
if (!compute) return(decline("not_requested")) # field present, not computed
blocks <- if (is.list(prior) && is.null(prior$type)) prior else list(prior)
if (length(blocks) != 1L) return(decline("not_applicable", "multi-block grid"))
blk <- blocks[[1L]]
# The multi-block dispatch does not carry a `type`, so resolve it from the block
# itself for the length-1 case (a single latent block wrapped in a list, or the
# model-supplied `likelihood` path). Without this the registry lookup below is
# `.NL_REGISTRY[[NULL]]`, an error rather than NULL.
type <- type %||% tolower(blk$type)
# A fit that relied on the default grid carries no `*_grid` field on the input
# prior, but the realised res$theta_grid is a valid single positive-scale axis
# regardless. Fill the type's defaults (the same ones .nl_dispatch applied) so
# eligibility keys on the actual grid, not on whether the user named it. A
# multi-axis type then exposes >1 grid field and still declines below.
spec <- .NL_REGISTRY[[type]]
if (!is.null(spec) && is.function(spec$defaults)) {
blk <- tryCatch(spec$defaults(blk, cargs), error = function(e) blk)
}
gfs <- grep("_grid$", names(blk), value = TRUE)
gfs <- gfs[vapply(gfs, function(f) is.numeric(blk[[f]]) && length(blk[[f]]) >= 2L,
logical(1))]
if (length(gfs) != 1L) {
return(decline("not_applicable",
if (length(gfs)) "multi-axis grid" else "no resolvable grid axis"))
}
if (!(gfs %in% .NL_POS_GRID)) {
# A bounded axis (a correlation `rho_grid`, say): declining beats
# transforming it with a guessed support -- and it will decline for every
# fit of this family, which is what the reason says.
return(decline("unguessable_axis", gfs))
}
tg <- as.numeric(res$theta_grid)
if (length(tg) != length(res$weights)) {
return(decline("internal_inconsistency", "weights do not match theta_grid"))
}
if (any(tg <= 0)) {
return(decline("internal_inconsistency",
paste0("non-positive node on the log axis ", gfs)))
}
# The batch is the diagnostic's own importance sample, not a grid the caller
# chose, so it runs with the internal-batch flag set: the cell-count advice
# would name a number `control$k_samples` sets and offer a remedy that does
# not reach it (gcol33/tulpa#614).
# Both dispatches read the integrated axes off the fit's own block, the
# multi-block one with the draws as its grid, so the draws are scored against
# the hyperprior the grid carries.
refit <- function(theta_mat) {
blk2 <- blk; blk2[[gfs]] <- as.numeric(theta_mat[, 1L])
lm <- tryCatch(
.nl_with_internal_batch(
if (dispatch_kind == "multi")
.nl_dispatch_multi(cargs, list(blk), likelihood = likelihood,
theta_grid_override = theta_mat)$log_marginal
else .nl_dispatch(type, cargs, blk2, blk)$log_marginal),
error = function(e) rep(-Inf, nrow(theta_mat)))
if (length(lm) != nrow(theta_mat)) rep(-Inf, nrow(theta_mat)) else lm
}
kd <- .with_preserved_seed(tryCatch(
.nested_grid_pareto_k(matrix(log(tg), ncol = 1L), res$weights, refit,
n_samples, tail_points = tail_points),
error = function(e) NULL))
if (is.null(kd)) return(decline("degenerate_proposal", "the scorer errored"))
res$pareto_k <- kd$pareto_k
res$pareto_k_is_ess <- kd$is_ess
res$pareto_k_proposal_source <- kd$proposal_source %||% NA_character_
res$pareto_k_first_pass <- kd$first_pass_k %||% NA_real_
# The GPD tail the shape was fitted on. A raised or capped tail was
# invisible on this path (gcol33/tulpa#692).
res$pareto_k_tail_points <- kd$tail_points %||% NA_integer_
.k_attach_declined(res, kd)
}
# --- Per-prior registry & dispatcher ---
#
# Each entry describes one supported prior `type`. To add a new prior:
# 1. Add a row here with cpp_fn, defaults, pack, theta.
# 2. Make sure cpp_nested_laplace_<variant> exists in RcppExports.
# No changes needed to .nl_dispatch() or the public driver.
# --- Required block fields, per dispatch path -------------------------------
#
# Each registry entry declares which fields its converter indexes, keyed by the
# path that consumes them. The C++ kernels read those vectors positionally with
# no bounds check, so a block naming a field wrongly -- a typo, a stale name
# after a rename, a block copied from a different family -- arrives with a
# zero-length vector and takes the session down. Checking
# presence at the R boundary turns that into a message naming the fields.
#
# Path keys:
# axis read by defaults() / theta() while the outer grid is built, so
# checked on every path.
# single extra fields pack() indexes for the single-block kernel.
# multi extra fields .nl_block_spec_for_cpp() indexes.
# joint extra fields .joint_block_spec_for_cpp() indexes.
# joint_single
# extra fields .joint_call_kernel_via_multi() indexes on the
# single-block joint path, where the per-observation index lives on
# each arm rather than on the block.
# projected
# the fields a projected block carries on a multi-block path: an areal
# block's graph or an NNGP block's neighbour structure, with
# `projector` in place of `spatial_idx`.
# A key absent from an entry means the type is not dispatchable on that path;
# the path's own converter raises its "only supported inside ..." message.
.NL_REQ_AREAL_GRAPH <- c("n_spatial_units", "adj_row_ptr", "adj_col_idx",
"n_neighbors")
.NL_REQ_AREAL <- c("spatial_idx", .NL_REQ_AREAL_GRAPH)
.NL_REQ_TEMPORAL <- c("temporal_idx", "n_times")
.NL_REQ_NNGP_FIELD <- c("coords", "nn_idx", "nn_dist", "n_spatial", "nn")
.NL_REQ_HSGP_ARM <- c("m_total", "phi", "n_obs_per_arm", "eigenvalues")
.NL_REQ_SPDE_MESH <- c("n_mesh", "A_x", "A_i", "A_p",
"C0_diag", "G1_x", "G1_i", "G1_p", "nu")
# `paths` is one or more keys above; their union is checked. A field that is
# present but empty counts as missing -- a zero-length vector is exactly what
# the kernel indexes out of bounds.
.nl_check_block_fields <- function(p, paths, block_index = NULL) {
type <- tolower(p$type %||% "")
spec <- .NL_REGISTRY[[type]]
if (is.null(spec) || is.null(spec$required)) return(invisible(p))
req <- unique(unlist(spec$required[paths], use.names = FALSE))
if (!length(req)) return(invisible(p))
bad <- req[vapply(req, function(f) length(p[[f]]) == 0L, logical(1))]
if (!length(bad)) return(invisible(p))
stop("prior block ", if (is.null(block_index)) "" else paste0(block_index, " "),
"'", type, "' is missing required field(s): ",
paste(bad, collapse = ", "), ".", call. = FALSE)
}
# Areal CSR adjacency block -- shared by icar/bym2/car_proper.
.nl_adj_args <- function(p) list(
spatial_idx = as.integer(p$spatial_idx),
n_spatial_units = as.integer(p$n_spatial_units),
adj_row_ptr = as.integer(p$adj_row_ptr),
adj_col_idx = as.integer(p$adj_col_idx),
n_neighbors = as.integer(p$n_neighbors)
)
# Default outer grid for a p-field separable MCAR, returned in column-major
# lower-triangle log-Cholesky coordinates of Sigma = L L' (the columns the
# integrator reads). For the common p = 2 case the nodes are placed in the
# interpretable (sigma_1, sigma_2, rho) space -- with nodes AT sensible
# correlation values -- and converted to log-Cholesky, so the rho axis is
# covered evenly (raw log-Cholesky off-diagonal nodes map to awkward, uneven
# correlations and can miss the mode under sharp likelihoods). For p > 2 a raw
# log-Cholesky tensor is used (coarser; the off-diagonal correlations are not
# separately interpretable).
.mcar_default_logchol_grid <- function(p) {
if (p == 2L) {
sig <- .nl_grid_axis("mcar_sd")
rho <- .nl_grid_axis("mcar_rho")
g <- expand.grid(s1 = sig, s2 = sig, rho = rho,
KEEP.OUT.ATTRS = FALSE)
# Sigma = [[s1^2, rho s1 s2], [rho s1 s2, s2^2]] -> lower Cholesky L:
# L11 = s1, L21 = rho s2, L22 = s2 sqrt(1 - rho^2).
out <- cbind(L11 = log(g$s1),
L21 = g$rho * g$s2,
L22 = log(g$s2 * sqrt(1 - g$rho^2)))
colnames(out) <- c("L11", "L21", "L22")
return(out)
}
diag_nodes <- log(.nl_grid_axis("mcar_logchol_diag"))
off_nodes <- .nl_grid_axis("mcar_logchol_off")
m <- p * (p + 1L) / 2L
axes <- vector("list", m)
nm <- character(m)
t <- 1L
for (j in seq_len(p)) for (i in j:p) {
axes[[t]] <- if (i == j) diag_nodes else off_nodes
nm[t] <- sprintf("L%d%d", i, j)
t <- t + 1L
}
g <- as.matrix(do.call(expand.grid, c(axes, list(KEEP.OUT.ATTRS = FALSE))))
colnames(g) <- nm
g
}
.NL_REGISTRY <- list(
icar = list(
cpp_fn = "cpp_nested_laplace_icar",
required = list(single = .NL_REQ_AREAL, multi = .NL_REQ_AREAL,
joint = .NL_REQ_AREAL,
projected = .NL_REQ_AREAL_GRAPH,
joint_single = .NL_REQ_AREAL_GRAPH),
defaults = function(p, a) .nl_fill_family_axes(p, "icar"),
pack = function(p) c(.nl_adj_args(p), list(
tau_grid = as.numeric(p$tau_grid)
)),
theta = function(p) list(grid = as.numeric(p$tau_grid), names = "tau")
),
bym2 = list(
cpp_fn = "cpp_nested_laplace_bym2",
required = list(single = .NL_REQ_AREAL, multi = .NL_REQ_AREAL,
joint = .NL_REQ_AREAL,
projected = .NL_REQ_AREAL_GRAPH,
joint_single = .NL_REQ_AREAL_GRAPH),
defaults = function(p, a) .nl_fill_family_axes(p, "bym2"),
pack = function(p) c(.nl_adj_args(p), list(
scale_factor = as.numeric(p$scale_factor %||% 1.0),
sigma_spatial_grid = as.numeric(p$sigma_grid),
rho_grid = as.numeric(p$rho_grid),
node_prec = if (!is.null(p$node_prec)) as.numeric(p$node_prec)
)),
theta = function(p) list(
grid = cbind(sigma = p$sigma_grid, rho = p$rho_grid),
names = c("sigma", "rho")
)
),
car_proper = list(
cpp_fn = "cpp_nested_laplace_car_proper",
required = list(single = .NL_REQ_AREAL, multi = .NL_REQ_AREAL,
joint = .NL_REQ_AREAL,
projected = .NL_REQ_AREAL_GRAPH,
joint_single = .NL_REQ_AREAL_GRAPH),
defaults = function(p, a) {
# `rho_car_grid` is the joint-API spelling of the correlation axis; accept
# it as an alias so a grid passed under either name is honoured.
p$rho_grid <- p$rho_grid %||% p$rho_car_grid
# Both axes supplied -> treat as the pre-paired integration grid (the
# theta builder pairs them with cbind), matching bym2 / ar1.
if (!is.null(p$tau_grid) && !is.null(p$rho_grid)) {
return(.nl_check_paired_axes(p, c("tau_grid", "rho_grid")))
}
# Otherwise default each missing axis independently and cross the two, so
# supplying only one axis keeps it instead of discarding it.
rb <- p$rho_bounds %||% .nl_grid_par("car_rho", "bounds")
# Margin in from the eigenvalue endpoints -- Q goes singular at the
# boundary, so anchoring the grid in (lower + eps, upper - eps) avoids
# NaN log-determinants from the per-grid-point Cholesky.
eps <- .nl_grid_par("car_rho", "margin") * (rb[2] - rb[1])
g_tau <- if (is.null(p$tau_grid))
.nl_grid_axis("car_tau") else sort(unique(p$tau_grid))
g_rho <- if (is.null(p$rho_grid))
seq(rb[1] + eps, rb[2] - eps,
length.out = .nl_grid_par("car_rho", "n")) else sort(unique(p$rho_grid))
gr <- expand.grid(tau = g_tau, rho = g_rho)
p$tau_grid <- gr$tau; p$rho_grid <- gr$rho
p
},
pack = function(p) c(.nl_adj_args(p), list(
tau_grid = as.numeric(p$tau_grid),
rho_grid = as.numeric(p$rho_grid)
)),
theta = function(p) list(
grid = cbind(tau = p$tau_grid, rho = p$rho_grid),
names = c("tau", "rho")
)
),
mcar = list(
# Separable multivariate CAR (p coupled areal fields sharing Sigma (x) Q^-1).
# Multi-block-only (like iid): the coupled inner solve lives in the joint
# driver. The outer axes are the p(p+1)/2 log-Cholesky coordinates of Sigma.
cpp_fn = NULL,
required = list(axis = "n_fields",
joint = c(.NL_REQ_AREAL, "n_fields", "field_weight")),
defaults = function(p, a) {
if (is.null(p$logchol_grid)) {
p$logchol_grid <- .mcar_default_logchol_grid(as.integer(p$n_fields))
}
p
},
pack = function(p) stop(
"MCAR is only supported inside a multi-block prior (a coupled areal ",
"field declared by spatial(graph, ~ ... | cell)).", call. = FALSE
),
theta = function(p) {
g <- p$logchol_grid
list(grid = g, names = colnames(g))
}
),
miid = list(
# Multivariate IID: p coupled per-group coefficient
# fields sharing a free Sigma -- the non-spatial sibling of mcar (Q = I, so
# no graph). Multi-block-only (like mcar / iid): the coupled inner solve
# lives in the joint driver. The outer axes are the same p(p+1)/2
# log-Cholesky coordinates of Sigma; only the per-field precision (I vs
# D - W) differs, which does not change the outer grid, so the default grid
# is reused from mcar.
cpp_fn = NULL,
required = list(axis = "n_fields",
joint = c("obs_idx", "n_groups", "n_fields", "field_weight")),
defaults = function(p, a) {
if (is.null(p$logchol_grid)) {
p$logchol_grid <- .mcar_default_logchol_grid(as.integer(p$n_fields))
}
p
},
pack = function(p) stop(
"miid is only supported inside a multi-block prior (a correlated ",
"random-slope term, e.g. (1 + x | group)).", call. = FALSE
),
theta = function(p) {
g <- p$logchol_grid
list(grid = g, names = colnames(g))
}
),
nngp = list(
cpp_fn = "cpp_nested_laplace_nngp",
# spatial_idx / nn_order / cov_type are defaulted on the single-block path,
# so they are not required there; the multi-block path reads spatial_idx.
required = list(single = .NL_REQ_NNGP_FIELD,
multi = c(.NL_REQ_NNGP_FIELD, "spatial_idx"),
joint = c(.NL_REQ_NNGP_FIELD, "spatial_idx"),
projected = .NL_REQ_NNGP_FIELD),
defaults = function(p, a) .nl_fill_family_axes(p, "nngp"),
pack = function(p) {
# spatial_idx (1-based, length N) maps obs -> spatial unit. Default to
# 1..n_spatial when N == n_spatial (the legacy one-obs-per-location case).
spatial_idx <- if (!is.null(p$spatial_idx)) {
as.integer(p$spatial_idx)
} else {
seq_len(as.integer(p$n_spatial))
}
c(list(spatial_idx = spatial_idx),
.nl_nngp_field_args(p),
list(sigma2_grid = as.numeric(p$sigma2_grid),
phi_gp_grid = as.numeric(p$phi_gp_grid)))
},
theta = function(p) list(
grid = cbind(sigma2 = p$sigma2_grid, phi_gp = p$phi_gp_grid),
names = c("sigma2", "phi_gp")
)
),
spde = list(
# Continuous Matern SPDE field on a shared FEM mesh. INDEXED_MULTI: each
# observation projects onto ~3 mesh nodes through the FEM matrix A, so the
# block is consumed via obs_indices (the A CSC slots), not a one-node
# spatial_idx. Same SpdeQBuilder and make_spde_block factory the
# single-Laplace occupancy SPDE path (cpp_nested_laplace_spde) uses.
cpp_fn = NULL, # multi-block only; reached via .nl_dispatch_multi
# The single-block converter reads a scalar `n_obs`; the joint one reads a
# per-arm `n_obs_per_arm` alongside per-arm A CSC triples.
required = list(multi = c(.NL_REQ_SPDE_MESH, "n_obs"),
joint = c(.NL_REQ_SPDE_MESH, "n_obs_per_arm")),
defaults = function(p, a) {
if (is.null(p$range_grid) || is.null(p$sigma_grid)) {
# Geometric grid centred on the PC-prior mode (prior_range[1] /
# prior_sigma[1]). The span is kept tight (mode / SPAN .. mode * SPAN)
# rather than the wide [mode*0.3, mode*3] window: a rectangular
# Cartesian grid (no CCD mode-find here) that runs to a very small
# range / large sigma corner lets the binary-occupancy field over-fit
# there, and the corner cell -- being the highest inner likelihood --
# dominates the weighted field even though the PC prior disfavours it.
# Anchoring the grid near the regularised prior mode keeps the
# integrated field in the well-identified region the single-Laplace
# SPDE path fits at the fixed mode, while still integrating modest
# hyperparameter uncertainty.
n_g <- as.integer(p$n_grid %||% .nl_grid_par("spde_registry", "n"))
span <- as.numeric(p$grid_span %||% .nl_grid_par("spde_registry", "span"))
rng_mode <- (.hp_range_anchor(p)$anchor %||%
.nl_grid_par("spde_registry", "centre_range"))[1]
sig_mode <- (p$prior_sigma %||% .nl_scale_anchor())[1]
rg <- exp(seq(log(rng_mode / span), log(rng_mode * span), length.out = n_g))
sg <- exp(seq(log(sig_mode / span), log(sig_mode * span), length.out = n_g))
gr <- expand.grid(range = rg, sigma = sg)
p$range_grid <- gr$range
p$sigma_grid <- gr$sigma
}
p
},
pack = function(p) stop(
"spde is only supported inside a multi-block latent prior on the ",
"nested-Laplace path. Pass it as a one-element list to ",
"tulpa_nested_laplace(prior = list(blk), ...).",
call. = FALSE
),
theta = function(p) list(
grid = cbind(range = as.numeric(p$range_grid),
sigma = as.numeric(p$sigma_grid)),
names = c("range", "sigma")
)
),
iid = list(
# Unstructured (iid) Gaussian latent, integrated over a single sigma axis.
# Multi-block-only: there is no single-arm iid kernel (a lone iid field is a
# plain random intercept). Reached via .nl_dispatch_multi -- the C++
# build_blocks_from_spec builds the block with d_fac(k) = sigma_k and an
# N(0, 1) prior. A one-point sigma_grid conditions the field at a FIXED
# sigma (an EM-estimated random-effect term on a nested block carries its
# own grid this way; see .fit_block_via_nested_laplace).
cpp_fn = NULL,
required = list(multi = c("obs_idx", "n_units"),
joint = c("obs_idx", "n_units")),
defaults = function(p, a) .nl_fill_family_axes(p, "iid"),
pack = function(p) stop(
"iid is only supported inside a multi-block latent prior on the ",
"nested-Laplace path. Pass it as one element of a block list to ",
"tulpa_nested_laplace(prior = list(...), ...).",
call. = FALSE
),
theta = function(p) list(grid = as.numeric(p$sigma_grid), names = "sigma")
),
hsgp = list(
cpp_fn = "cpp_nested_laplace_hsgp",
# The single-block kernel takes the basis as one matrix (phi_basis) with its
# eigenvalues (lambda_eig); the joint one takes a per-arm list (phi) with
# `eigenvalues` and the basis / row counts it needs to slice them.
required = list(single = c("phi_basis", "lambda_eig"),
multi = c("phi_basis", "lambda_eig"),
joint = .NL_REQ_HSGP_ARM),
defaults = function(p, a) .nl_fill_family_axes(p, "hsgp"),
pack = function(p) list(
phi_basis = as.matrix(p$phi_basis),
lambda_eig = as.numeric(p$lambda_eig),
sigma2_grid = as.numeric(p$sigma2_grid),
lengthscale_grid = as.numeric(p$lengthscale_grid)
),
theta = function(p) list(
grid = cbind(sigma2 = p$sigma2_grid, lengthscale = p$lengthscale_grid),
names = c("sigma2", "lengthscale")
)
),
hsgp_mo = list(
# Multi-output (co-regionalization) HSGP block. K = n_arms
# correlated latent fields share basis + eigenvalues, cross-output
# coefficients live in Sigma. First ship: K = 2 with raw axes
# (sigma_1, sigma_2, rho, ell) and the natural Cartesian product over
# paired per-axis grids. Multi-block-only -- there is no single-arm
# version of multi-output HSGP.
cpp_fn = NULL,
required = list(joint = .NL_REQ_HSGP_ARM),
defaults = function(p, a) {
# All four grids treated as paired axes (each row of theta_grid is
# one (sigma_1, sigma_2, rho, ell) tuple). When any grid is missing
# we fill in a small Cartesian default; when the user supplies grids
# they must already be aligned (use expand.grid() on the R side to
# form a Cartesian, then pass the four aligned columns).
.nl_fill_family_axes(p, "hsgp_mo")
},
pack = function(p) stop(
"hsgp_mo is only supported inside a multi-block joint prior. ",
"Wrap the block in list() and pass to tulpa_nested_laplace_joint(prior = list(blk), ...).",
call. = FALSE
),
theta = function(p) list(
grid = cbind(sigma_1 = as.numeric(p$sigma_1_grid),
sigma_2 = as.numeric(p$sigma_2_grid),
rho = as.numeric(p$rho_grid),
ell = as.numeric(p$lengthscale_grid)),
names = c("sigma_1", "sigma_2", "rho", "ell")
)
),
rw1 = list(
cpp_fn = "cpp_nested_laplace_temporal",
required = list(single = .NL_REQ_TEMPORAL, multi = .NL_REQ_TEMPORAL,
joint = .NL_REQ_TEMPORAL),
defaults = function(p, a) .nl_fill_family_axes(p, "rw1"),
pack = function(p) list(
temporal_idx = as.integer(p$temporal_idx),
n_times = as.integer(p$n_times),
temporal_type = "rw1",
tau_grid = as.numeric(p$tau_grid),
rho_grid = numeric(0),
cyclic = isTRUE(p$cyclic),
n_groups = as.integer(p$n_groups %||% 1L)
),
theta = function(p) list(grid = as.numeric(p$tau_grid), names = "tau")
),
rw2 = list(
cpp_fn = "cpp_nested_laplace_temporal",
required = list(single = .NL_REQ_TEMPORAL, multi = .NL_REQ_TEMPORAL,
joint = .NL_REQ_TEMPORAL),
defaults = function(p, a) .nl_fill_family_axes(p, "rw2"),
pack = function(p) list(
temporal_idx = as.integer(p$temporal_idx),
n_times = as.integer(p$n_times),
temporal_type = "rw2",
tau_grid = as.numeric(p$tau_grid),
rho_grid = numeric(0),
cyclic = FALSE,
n_groups = as.integer(p$n_groups %||% 1L)
),
theta = function(p) list(grid = as.numeric(p$tau_grid), names = "tau")
),
ar1 = list(
cpp_fn = "cpp_nested_laplace_temporal",
required = list(single = .NL_REQ_TEMPORAL, multi = .NL_REQ_TEMPORAL,
joint = .NL_REQ_TEMPORAL),
defaults = function(p, a) {
# The temporal kernel takes rho as given, so an AR1 correlation outside
# (-1, 1) -- a non-stationary process -- is refused here, by name.
.nl_check_axis_open_interval(p, "rho_grid", -1, 1,
"an AR1 correlation")
.nl_fill_family_axes(p, "ar1")
},
pack = function(p) list(
temporal_idx = as.integer(p$temporal_idx),
n_times = as.integer(p$n_times),
temporal_type = "ar1",
tau_grid = as.numeric(p$tau_grid),
rho_grid = as.numeric(p$rho_grid),
cyclic = FALSE,
n_groups = as.integer(p$n_groups %||% 1L)
),
theta = function(p) list(
grid = cbind(tau = p$tau_grid, rho = p$rho_grid),
names = c("tau", "rho")
)
),
lf = list(
# Latent factor block. F = 1 only. No outer-grid axes:
# identifiability is handled by tight Gaussian anchors on (u_1,
# lambda_1) inside the C++ factory; everything else is jointly
# optimized by the inner Newton. Returns a 1 x 0 grid so the joint
# Cartesian product treats this block as a "no axis" contributor.
cpp_fn = NULL,
required = list(joint = c("obs_idx", "n_latent")),
defaults = function(p, a) p,
pack = function(p) stop(
"lf is only supported inside a multi-block prior. ",
"Wrap the block in list(): tulpa_nested_laplace_joint(prior = list(blk), ...).",
call. = FALSE
),
theta = function(p) list(
grid = matrix(numeric(0), nrow = 1L, ncol = 0L),
names = character(0)
)
),
tgmrf = list(
# User-defined GMRF block (see ?tgmrf). Only available inside a multi-
# block prior. The outer theta-grid is the Cartesian product of
# per-axis grids; per-axis defaults are 5 equispaced points between
# block$bounds$lower[j] and block$bounds$upper[j], or block$init[j] +/-
# 2 if no bounds are supplied. Precomputed Q_k / logdet / log p(theta_k)
# at every grid row are packed into the C++ block spec by
# `.nl_block_spec_for_cpp`.
cpp_fn = NULL,
# The multi-block converter builds the per-grid Q / logdet / log-prior from
# the block's own closures; the joint one takes them already materialised.
required = list(
axis = c("theta_dim", "theta_names"),
multi = c("n_latent", "obs_idx"),
joint = c("n_latent", "obs_idx", "Q_csc_p_per_grid", "Q_csc_i_per_grid",
"Q_csc_x_per_grid", "logdet_Q_per_grid",
"log_prior_theta_per_grid")),
defaults = function(p, a) {
if (is.null(p$theta_grid_built)) {
d <- p$theta_dim
axes <- vector("list", d)
for (j in seq_len(d)) {
if (!is.null(p$bounds)) {
lo <- p$bounds$lower[j]
hi <- p$bounds$upper[j]
} else {
hw <- .nl_grid_par("tgmrf_axis", "half_width")
lo <- p$init[j] - hw
hi <- p$init[j] + hw
}
axes[[j]] <- seq(lo, hi,
length.out = .nl_grid_par("tgmrf_axis", "n"))
}
names(axes) <- p$theta_names
gr <- do.call(expand.grid, axes)
p$theta_grid_built <- as.matrix(gr)
}
p
},
pack = function(p) stop(
"tgmrf is only supported inside a multi-block prior. ",
"Wrap the block in list(): tulpa_nested_laplace(prior = list(blk), ...).",
call. = FALSE
),
theta = function(p) list(
grid = p$theta_grid_built,
names = p$theta_names
)
)
)
.nl_dispatch <- function(type, a, p, declared) {
spec <- .NL_REGISTRY[[type]]
if (is.null(spec)) {
stop("Unknown prior type: ", type,
". Supported: ", paste(names(.NL_REGISTRY), collapse = ", "), ".",
call. = FALSE)
}
# A multi-block-only type has no single-arm kernel; let its own pack() raise
# the "only supported inside a multi-block prior" message rather than a
# field-name list for fields the single path would never read.
if (!is.null(spec$cpp_fn)) .nl_check_block_fields(p, c("axis", "single"))
hyperprior <- .hp_choice(a$hyperprior)
a$hyperprior <- NULL
# A refinement round's screen record (`.hyper_refine_screen()`): the cells
# are a batch joining a solved grid. Read here, never handed to a kernel.
screen <- a$screen_refine
a$screen_refine <- NULL
p <- spec$defaults(p, a)
th <- spec$theta(p)
tg <- .nl_theta_matrix(list(theta_grid = th$grid, theta_names = th$names))
# `declared` is the block whose grid the fit integrates: `p` itself for a grid
# solve, the fit's block when `p` carries a stencil or a probe row written onto
# its grid fields. The integrated columns are read off that grid, since the
# rows written onto `p` vary a column the fit holds fixed or hold one it
# integrates at a single value (gcol33/tulpa#760).
dth <- spec$theta(spec$defaults(declared, a))
hp_declared <- .hp_declare(
.nl_theta_matrix(list(theta_grid = dth$grid, theta_names = dth$names)))
# The block's hyperprior, folded where every caller of the kernel reads it:
# the grid solve, the placement refit and its stencil, the k-hat refit. A
# screened solve ranks its cells with it too.
hp <- list(.nl_block_log_hyperprior(p, tg, hyperprior, declared = hp_declared))
if ((a$prune_tol %||% 0) > 0) {
a$screen_log_offset <- .nl_screen_log_offset(tg, hp, screen = screen)
}
out <- do.call(spec$cpp_fn, c(spec$pack(p), a))
out$theta_grid <- th$grid
out$theta_names <- th$names
.nl_fold_hyperprior(out, hp)
}
# Per-cell log weight the cheap screen ranks with: the hyperprior R folds into
# `log_marginal` after the kernel returns, plus the cell measure the grid's
# weights are built with (`.nl_attach_outer_integration()`), in the grid's cell
# order. Both are fixed by the grid before any cell is solved, and together with
# the kernel's log-marginal they are the posterior weight a cell ends up with,
# so the screen and the full pass rank the same quantity.
#
# `screen` is the record a refinement round hands its kernel
# (`.hyper_refine_screen()`): the cells are a batch joining a solved grid, so
# their measure is read on the merged grid (`log_measure`) rather than on the
# batch, and the grid's own mass (`log_ref`) travels as the offset's `log_ref`
# attribute, which the kernel's screen normalises against.
.nl_screen_log_offset <- function(tg, hp_parts, specs = NULL, screen = NULL) {
tg <- as.matrix(tg)
n <- nrow(tg)
lp <- numeric(n)
axes <- character(0)
for (hp in hp_parts) {
lp <- lp + hp$lp
axes <- union(axes, hp$axes)
}
lq <- if (!is.null(screen)) {
if (length(screen$log_measure) != n) {
stop("A refinement screen measures ", length(screen$log_measure),
" cells; the kernel call holds ", n, ".", call. = FALSE)
}
screen$log_measure
} else {
if (is.null(specs) && !is.null(colnames(tg))) {
specs <- .joint_axis_specs_from_grid(tg, folded_axes = axes)
}
.hyper_log_quad_weights(tg, specs)
}
off <- lp + (if (length(lq) == n) lq else 0)
off[is.na(off)] <- -Inf
if (!is.null(screen)) attr(off, "log_ref") <- as.numeric(screen$log_ref)
off
}
# The outer grid's cell measure, its weights and its log evidence, set together
# from the grid and the log-marginals the fit carries. `log_marginal` already
# holds the hyperprior density the dispatch folded in (`R/hyperprior_default.R`),
# so a folded axis is measured by its cell widths alone and its support is not
# clipped to the span the nodes were declared over.
.nl_attach_outer_integration <- function(res, what = "outer grid") {
tg <- .nl_theta_matrix(res)
# A node a refinement pass appended was not declared, so the specs are read
# off the base cells and the cells' tags carry the rest.
refining <- res$refining_axis
base <- !nzchar(.hyper_slice_home(refining, nrow(tg)))
specs <- .joint_axis_specs_from_grid(
tg[base, , drop = FALSE], folded_axes = res$log_hyperprior_axes)
res$log_quad <- .hyper_log_quad_weights(tg, specs, refining = refining)
res$axis_support <- .hyper_grid_supports(tg, specs, refining = refining)
res$weights <- .nl_normalise_weights_safe(res$log_marginal, what,
log_quad = res$log_quad)
res <- .nl_prune_merged_record(res)
.nl_attach_evidence(res, tg, specs)
}
# Normalise log-marginals to integration weights summing to 1. The single
# weight normaliser for every outer hyperparameter grid: it drops non-finite
# nodes (an inner Newton diverging in a grid corner returns +Inf / NaN) from
# the max-shift, zeroes their weight, and renormalises over the finite cells;
# all-NA (with a warning) when no node carries mass, never NaN. There is no
# unguarded variant -- a single `max(lm)` over a grid with one +Inf / NaN cell
# would collapse every weight to NaN and break the grid sampler. `what` names
# the grid in the degenerate-case warning.
# `log_quad` adds the per-cell log quadrature weight of the outer grid (see
# R/hyper_quadrature.R), so posterior weight is likelihood times prior mass
# rather than likelihood alone. Callers integrating over a hyperparameter grid
# pass it; callers weighting something that is not a quadrature rule leave it
# NULL.
#' Normalise outer-grid log-marginals to posterior cell weights
#'
#' The single weight normaliser for every outer hyperparameter grid: softmax
#' of `lm` (optionally plus a per-cell log quadrature weight) into weights
#' summing to 1. Non-finite nodes (an inner Newton diverging in a grid
#' corner) are dropped from the max-shift and zeroed before renormalising
#' over the finite cells; an all-non-finite grid returns all-`NA` with a
#' warning rather than `NaN`. Intended for reconstructing a fit's outer-grid
#' posterior weights (with [tulpa_theta_matrix()] and
#' [tulpa_grid_log_quad()]) when the fit doesn't already carry `fit$weights`.
#'
#' @param lm Numeric vector of per-cell log-marginal-likelihood values.
#' @param what Label for the grid, used only in the degenerate-case warning.
#' @param log_quad Optional per-cell log quadrature weight (prior mass) of
#' the outer grid, e.g. from [tulpa_grid_log_quad()]; `NULL` normalises the
#' likelihood alone.
#' @return Numeric vector of posterior cell weights summing to 1, the same
#' length as `lm`, or all-`NA` if no cell carries finite mass.
#' @seealso [tulpa_theta_matrix()], [tulpa_grid_log_quad()]
#' @examples
#' tulpa_normalise_weights_safe(c(-1, 0, -2))
#' # With per-cell log quadrature weights (cell prior masses):
#' tulpa_normalise_weights_safe(c(-1, 0, -2), log_quad = log(c(0.2, 0.5, 0.3)))
#' @export
tulpa_normalise_weights_safe <- function(lm, what = "grids / data",
log_quad = NULL) {
.nl_normalise_weights_safe(lm, what = what, log_quad = log_quad)
}
.nl_normalise_weights_safe <- function(lm, what = "grids / data",
log_quad = NULL) {
if (!is.null(log_quad)) {
if (length(log_quad) != length(lm)) {
stop(sprintf("`log_quad` has length %d but `lm` has length %d (%s).",
length(log_quad), length(lm), what), call. = FALSE)
}
lm <- lm + log_quad
lm[is.na(lm)] <- -Inf
}
finite_lm <- lm[is.finite(lm)]
if (length(finite_lm) == 0L) {
warning(sprintf("All grid points returned non-finite log_marginal -- check %s.", what),
call. = FALSE)
return(rep(NA_real_, length(lm)))
}
w <- exp(lm - max(finite_lm))
w[!is.finite(w)] <- 0
if (sum(w) == 0) return(rep(NA_real_, length(lm)))
w / sum(w)
}
# Log evidence of an outer hyperparameter grid.
#
# `log_marginal` carries a normalised hyperprior density on each axis's
# integration coordinate and `log_measure` is each cell's ABSOLUTE log volume on
# that coordinate (its widths, and a declared atom's probability), so
#
# log p(y) = log sum_k exp(log_marginal_k + log_measure_k)
#
# is the evidence under the prior itself, not under a copy of it restricted to
# the grid and renormalised there. Where the grid covers the posterior the
# value does not move with where the nodes were laid or how many there are.
#
# A cell whose inner solve failed (`log_marginal` NaN, or `+Inf` from a non-PD
# Hessian) contributes nothing: its likelihood is unknown and is not guessed.
# A cell screened out by the pruning pass carries `-Inf` and contributes
# nothing either, which is what screening it asserted.
.nl_outer_log_evidence <- function(log_marginal, log_measure) {
lq <- as.numeric(log_measure)
if (length(lq) != length(log_marginal)) {
stop(sprintf(paste0("Outer evidence: `log_marginal` has length %d but the ",
"measure has %d."),
length(log_marginal), length(lq)), call. = FALSE)
}
known <- !is.na(log_marginal) & log_marginal < Inf & !is.na(lq)
if (!any(known)) return(NA_real_)
ev <- .tulpa_logsumexp(log_marginal[known] + lq[known])
if (is.finite(ev)) ev else NA_real_
}
# Record a grid's evidence on `res`, or why it has none. The evidence exists
# only when every integrated axis carries a normalised density; otherwise it is
# `NA` with `log_evidence_declined = "improper_hyperprior"` and
# `log_evidence_declined_axes` naming each such axis and its reason.
.nl_attach_evidence <- function(res, tg, specs,
log_measure = .hyper_log_quad_weights(
tg, specs, absolute = TRUE,
refining = res$refining_axis)) {
uncovered <- .nl_evidence_uncovered(tg, specs, res)
res$log_evidence_declined_axes <- if (length(uncovered)) uncovered else NULL
if (length(uncovered)) {
res$log_evidence <- NA_real_
res$log_evidence_declined <- "improper_hyperprior"
return(res)
}
n <- length(res$log_marginal)
if (!is.null(tg) && nrow(as.matrix(tg)) == n &&
!length(.hp_integrated_axes(as.matrix(tg)))) {
# Every axis is pinned, so each is part of the model and the grid integrates
# nothing: its cells are copies of one point, which together carry measure 1.
log_measure <- rep(-log(n), n)
}
res$log_evidence <- .nl_outer_log_evidence(res$log_marginal, log_measure)
res$log_evidence_declined <- NULL
res
}
`%||%` <- function(x, y) if (is.null(x)) y else x
# ----------------------------------------------------------------------------
# Per-grid-point fixed-effects Hessian extraction.
#
# When tulpa_nested_laplace() is called with keep_grid_hessians = TRUE the
# inner Laplace solves are run with store_Q = TRUE, so each grid point's full
# posterior precision (over the joint latent vector x = (beta, RE, latent))
# is returned in CSC lower-triangle form. Here we
#
# 1. Reconstruct sparse symmetric Q_k from the CSC triple at grid point k.
# 2. Compute Sigma_{beta beta} = (Q_k^{-1})[1:p, 1:p] via a single sparse
# solve with the first p columns of the identity as RHS.
# 3. Invert that small p x p block to get the marginal fixed-effects
# Hessian H_beta = Sigma_{beta beta}^{-1}.
#
# Per-grid memory cost: O(p^2). The dominant compute is the sparse solve,
# already required for the integrated log-marginal -- so adding it here is
# essentially free.
# ----------------------------------------------------------------------------
# `keep_joint` additionally retains the MODAL cell's joint precision as
# `res$H_joint`, which is what the subspace-debias coupling closure reads to
# grow the corrected set over strongly coupled neighbours
# (`.subspace_closure()`). One cell, not `n_grid`: the closure is a SELECTION
# decision, and the selection is made once. Off by default, so a fit that did
# not ask for closure retains nothing extra (gcol33/tulpa#862).
# The MODAL cell's joint posterior precision, assembled from the per-cell CSC
# scratch the inner solves return under `store_Q`. This is what the
# subspace-debias coupling closure reads (`.subspace_closure()`): a coordinate
# strongly coupled to a member of the corrected set and left OUT of it is being
# carried linearly, which is precisely the error the correction removes, so the
# closure needs the precision graph to find it.
#
# ONE cell, not `n_grid`. The closure is a selection decision made once, and
# the modal cell is the one the reported mixture is centred on. Both drivers
# read this -- the single/multi-block grid through
# `.nl_attach_grid_hessians()`, the joint through `.joint_finalize_grid_fixed()`
# -- so the two cannot disagree about which cell the coupling came from
# (gcol33/tulpa#862). NULL when the fit retained no precision, which the caller
# records as a declined closure rather than as "no neighbours found".
.nl_modal_joint_precision <- function(res) {
Q_p <- res$Q_csc_p_per_grid
Q_i <- res$Q_csc_i_per_grid
Q_x <- res$Q_csc_x_per_grid
n_x <- res$Q_csc_n
if (is.null(Q_p) || is.null(Q_i) || is.null(Q_x) || is.null(n_x)) return(NULL)
lm <- res$log_marginal
if (is.null(lm) || !any(is.finite(lm))) return(NULL)
k <- which.max(replace(lm, !is.finite(lm), -Inf))
if (k > length(Q_p) || length(Q_p[[k]]) != n_x + 1L) return(NULL)
L <- Matrix::sparseMatrix(i = Q_i[[k]], p = Q_p[[k]], x = Q_x[[k]],
dims = c(n_x, n_x), index1 = FALSE)
L + Matrix::t(L) - Matrix::Diagonal(n_x, Matrix::diag(L))
}
#
# `var_cols` names latent coordinates whose per-cell MARGINAL variance,
# diag(Q_k^{-1}), is retained alongside (`grid_field_var`, n_grid x
# length(var_cols), with the columns in `grid_field_var_cols`). The same solve
# delivers it: those unit columns join E. It defaults to the temporal field (the
# `temporal =` block, or a latent(temporal_ar*()) block), so
# temporal() reports each cell's Gaussian -- mode AND curvature -- mixed over
# the grid rather than the cell modes alone, which carried only the spread
# between cells and covered the truth at 57% for a nominal 95%
# (gcol33/tulpa#904). Only a field whose read-back needs it is retained: one
# column per coordinate per cell, and a large spatial field is not asked for.
#
# The field's full within-cell covariance (`grid_field_cov`, a list of
# n_var x n_var matrices parallel to the cells, NULL where a cell was pruned)
# is retained too, so temporal()'s draws carry the correlation BETWEEN time
# points a cell's Gaussian implies and not only each point's marginal: a
# derived quantity over time (a trend, a difference between two years, a
# cumulative sum) read off draws with independent coordinates has the wrong
# spread. Capped at `.NL_FIELD_COV_MAX` doubles over the grid; past the cap only
# the marginals are kept and the draws fall back to independent coordinates.
.NL_FIELD_COV_MAX <- 4e6
.nl_attach_grid_hessians <- function(res, p_fixed, keep_joint = FALSE,
var_cols = c(
.nl_role_cols(res$blocks,
res$block_latent_offsets,
res$Q_csc_n, "temporal"),
.nl_temporal_latent_block(res)$cols)) {
Q_p <- res$Q_csc_p_per_grid
Q_i <- res$Q_csc_i_per_grid
Q_x <- res$Q_csc_x_per_grid
n_x <- res$Q_csc_n
if (is.null(Q_p) || is.null(Q_i) || is.null(Q_x) || is.null(n_x)) {
warning("keep_grid_hessians = TRUE but the inner solve did not return Q. ",
"This is a tulpa bug; please report.",
call. = FALSE)
return(res)
}
n_grid <- length(Q_p)
grid_hessians <- vector("list", n_grid)
grid_modes <- vector("list", n_grid)
var_cols <- as.integer(var_cols)
var_cols <- var_cols[var_cols > p_fixed & var_cols <= n_x]
n_var <- length(var_cols)
field_var <- if (n_var) matrix(NA_real_, n_grid, n_var) else NULL
# The field's whole within-cell covariance is the lower-right block of the
# same S, so retaining it costs no further solve -- only the storage, which
# grows as n_var^2 per cell and is capped by `.NL_FIELD_COV_MAX`.
keep_cov <- n_var > 0L && as.numeric(n_grid) * n_var^2 <= .NL_FIELD_COV_MAX
field_cov <- if (keep_cov) vector("list", n_grid)
E <- matrix(0, nrow = n_x, ncol = p_fixed + n_var)
E[cbind(c(seq_len(p_fixed), var_cols), seq_len(p_fixed + n_var))] <- 1
modes_mat <- res$modes # n_grid x n_x
# An intrinsic block (ICAR, BYM2's structured part, RW1 / RW2) stands for a
# field whose constant is a HARD sum-to-zero constraint. The inner solve
# augments that direction so the Hessian inverts, and centres it into the
# intercept, so the stored precision carries a finite 1 / (tau J) variance
# on the level -- and the unconditioned inverse puts it into the intercept's
# standard error (3-4x the exact posterior SD on a 50-node ICAR,
# gcol33/tulpa#901) and into an intrinsic temporal field's marginal variance.
# The driver names each intrinsic block's constraint group
# (`constraint_cols`), and both reads are conditioned on them by the kriging
# correction a joint fit's retained block and posterior draws use
# (cpp_joint_inner_vcov_blocks, Rue & Held 2005):
# Cov(x | A x = 0) = V - W (A W)^-1 W', V = Q^-1 E, W = Q^-1 A'.
# The mode needs no correction: the centring already put it on the
# constraint.
constr <- res$constraint_cols
n_con <- length(constr)
if (n_con) {
E <- cbind(E, vapply(constr, function(cc) {
a <- numeric(n_x)
a[as.integer(cc)] <- 1
a
}, numeric(n_x)))
}
sel <- seq_len(p_fixed + n_var)
for (k in seq_len(n_grid)) {
# A cell the cheap screen pruned is never solved, so it holds no precision
# to retain and its CSC slot is empty. Skip the write and the slot stays
# blank, which is what `.nested_fixed_moments()` already reads: it drops any
# cell whose block is absent, and a pruned cell carries log_marginal = -Inf
# and hence zero weight, so nothing the report reads was in it. Skipping
# rather than assigning NULL is what keeps `grid_modes` / `grid_hessians`
# parallel to `weights` -- assigning NULL removes the element and shortens
# the list, which NAs the whole coefficient table (gcol33/tulpa#345).
if (length(Q_p[[k]]) != n_x + 1L) next
grid_modes[[k]] <- if (!is.null(modes_mat)) {
as.numeric(modes_mat[k, seq_len(p_fixed)])
} else {
rep(NA_real_, p_fixed)
}
L <- Matrix::sparseMatrix(
i = Q_i[[k]], p = Q_p[[k]], x = Q_x[[k]],
dims = c(n_x, n_x), index1 = FALSE
)
diag_L <- Matrix::diag(L)
Qk <- L + Matrix::t(L) - Matrix::Diagonal(n_x, diag_L)
# V = Qk^{-1} E selects the fixed-effect (and read field) columns of the
# joint inverse. A dense GP / NNGP field precision at close locations makes
# the joint precision ill-conditioned, which trips Matrix's
# condition-number guard. Retry with a negligible diagonal jitter (scaled
# to the precision) so the marginal-SE extraction stays available rather
# than aborting the whole fit.
V <- tryCatch(
as.matrix(Matrix::solve(Qk, E)),
error = function(e) {
jit <- 1e-8 * mean(Matrix::diag(Qk))
as.matrix(Matrix::solve(Qk + Matrix::Diagonal(n_x, jit), E))
}
)
rows <- c(seq_len(p_fixed), var_cols)
S <- V[rows, sel, drop = FALSE]
if (n_con) {
W <- V[, p_fixed + n_var + seq_len(n_con), drop = FALSE]
AW <- crossprod(E[, p_fixed + n_var + seq_len(n_con), drop = FALSE], W)
Wr <- W[rows, , drop = FALSE]
S <- S - Wr %*% solve((AW + t(AW)) / 2, t(Wr))
}
Sigma_bb <- S[seq_len(p_fixed), seq_len(p_fixed), drop = FALSE]
grid_hessians[[k]] <- .spd_inverse(Sigma_bb)
if (n_var) {
fi <- p_fixed + seq_len(n_var)
field_var[k, ] <- diag(S)[fi]
if (keep_cov) {
Sf <- S[fi, fi, drop = FALSE]
field_cov[[k]] <- (Sf + t(Sf)) / 2
}
}
}
res$grid_hessians <- grid_hessians
res$grid_modes <- grid_modes
if (n_var) {
res$grid_field_var <- field_var
res$grid_field_var_cols <- var_cols
if (keep_cov) res$grid_field_cov <- field_cov
}
if (isTRUE(keep_joint)) res$H_joint <- .nl_modal_joint_precision(res)
# Strip the verbose CSC scratch fields once Hessians are assembled.
res$Q_csc_p_per_grid <- NULL
res$Q_csc_i_per_grid <- NULL
res$Q_csc_x_per_grid <- NULL
res$Q_csc_n <- NULL
res
}
# --- Multi-block dispatch ---
#
# A multi-block prior is `list(block1, block2, ...)` where each blockN is a
# list with a `type` field (same shape as a single-block prior). At each
# outer-grid point all blocks share a Newton solve; the joint grid is the
# Cartesian product of per-block axes.
# Block types whose field may reach the observations through a projector: arm
# k's predictor receives A_k z in place of the unit gather z[spatial_idx]
# (restricted spatial regression, A_k = P_k S_k). Both multi-block drivers read
# one; an areal block's level is placed from A_k 1, an NNGP block's prior
# carries its own.
.NL_PROJECTABLE_AREAL <- c("icar", "bym2", "car_proper")
.NL_PROJECTABLE_TYPES <- c(.NL_PROJECTABLE_AREAL, "nngp")
.nl_block_is_projected <- function(p) {
is.list(p) && tolower(p$type %||% "") %in% .NL_PROJECTABLE_TYPES &&
!is.null(p$projector)
}
# A block `projector` is read by the two multi-block drivers, the single-arm
# one behind tulpa_nested_laplace() and the joint one. The single-block joint
# path would gather the field through `spatial_idx` and drop the projector, so
# it refuses it.
.nl_refuse_projector <- function(prior, door) {
blocks <- if (is.list(prior) && !is.null(prior$type)) list(prior) else prior
if (!is.list(blocks)) return(invisible(NULL))
has <- vapply(blocks, function(b) is.list(b) && !is.null(b$projector),
logical(1))
if (any(has)) {
stop(door, " does not read a block `projector`; a projected areal field ",
"is fitted by tulpa_nested_laplace_joint() with a list of blocks, ",
"`prior = list(<block>)`.", call. = FALSE)
}
invisible(NULL)
}
# A projected block's per-arm projectors as the column-compressed triple
# the C++ block builders read (0-based `A_i` / `A_p`, the SPDE projector's
# layout). `p$projector` is a list of n_arms matrices, each N_k x n_units
# (dense or a Matrix sparse matrix), or one matrix shared by every arm; a NULL
# entry is an arm the field does not reach. The projector replaces the unit
# gather, so a block carrying `spatial_idx` as well, or a per-row
# `svc_weight` that belongs inside A_k, is refused.
.nl_block_projector <- function(p, n_arms, block_index, arm_n_obs) {
type <- tolower(p$type)
lab <- paste0("Block ", block_index, " (type '", type, "')")
if (!is.null(p$spatial_idx)) {
stop(lab, ": pass either `projector` or `spatial_idx`, not both; a ",
"projector already maps the observations to the units.",
call. = FALSE)
}
if (!is.null(p$svc_weight)) {
stop(lab, ": `projector` and `svc_weight` cannot be combined; scale ",
"the projector's rows by the weight.", call. = FALSE)
}
if (is.null(arm_n_obs)) {
stop(lab, ": a projector needs the arms' observation counts.",
call. = FALSE)
}
n_units <- as.integer(p$n_spatial_units %||% p$n_spatial)
projector <- p$projector
if (!is.list(projector)) projector <- rep(list(projector), n_arms)
if (length(projector) != n_arms) {
stop(lab, ": `projector` must be a list of length n_arms (", n_arms,
"), or one matrix shared by every arm.", call. = FALSE)
}
csc <- lapply(seq_len(n_arms), function(k) {
A <- projector[[k]]
if (is.null(A)) {
return(list(x = numeric(0), i = integer(0), p = integer(n_units + 1L)))
}
if (!(is.matrix(A) || methods::is(A, "Matrix")) ||
nrow(A) != arm_n_obs[k] || ncol(A) != n_units) {
stop(lab, ": `projector[[", k, "]]` must be a ", arm_n_obs[k],
" x ", n_units, " matrix (arm ", k, "'s observations by the ",
"block's units).", call. = FALSE)
}
A <- methods::as(methods::as(Matrix::Matrix(A, sparse = TRUE),
"generalMatrix"), "CsparseMatrix")
if (!all(is.finite(A@x))) {
stop(lab, ": `projector[[", k, "]]` holds a non-finite entry.",
call. = FALSE)
}
list(x = as.numeric(A@x), i = as.integer(A@i), p = as.integer(A@p))
})
list(A_x = lapply(csc, `[[`, "x"),
A_i = lapply(csc, `[[`, "i"),
A_p = lapply(csc, `[[`, "p"))
}
# An NNGP block's neighbour structure as the C++ block builders read it, on
# either the single-block kernel or the multi-block driver. nn_order is passed
# 0-based, the convention cpp_laplace_fit_gp also uses.
.nl_nngp_field_args <- function(p) {
n_spatial <- as.integer(p$n_spatial)
list(
coords = as.matrix(p$coords),
nn_idx = as.matrix(p$nn_idx),
nn_dist = as.matrix(p$nn_dist),
nn_order = as.integer((p$nn_order %||% seq_len(n_spatial)) - 1L),
n_spatial = n_spatial,
nn = as.integer(p$nn),
cov_type = as.integer(p$cov_type %||% 2L) # default Matern-5/2
)
}
.is_multi_block_prior <- function(p) {
is.list(p) && is.null(p$type) &&
all(vapply(p, function(x) is.list(x) && !is.null(x$type), logical(1)))
}
# What each block of a nested prior is FOR: "spatial", "temporal" or "other".
# tulpa() tags its `spatial =` / `temporal =` blocks with `role`, and on a tagged
# prior every untagged block (an s(x) smoother, a latent() block, a `(1 | g)`
# iid block) is "other" -- a type does not say what a block is for, since a
# smoother is an rw1 / rw2 block and a random intercept an iid one. A prior
# built by hand carries no tags, and its blocks are classified by type
# (`.SPATIAL_NL_TYPES` / `.TEMPORAL_NL_TYPES`). `prior` is a single block or a
# list of blocks; the result has one entry per block.
.nl_block_roles <- function(prior) {
if (is.null(prior) || !is.list(prior) || !length(prior)) return(character(0))
blocks <- if (!is.null(prior$type)) list(prior) else prior
tagged <- vapply(blocks, function(b) {
r <- if (is.list(b)) b$role else NULL
if (is.character(r) && length(r) == 1L) r else NA_character_
}, character(1))
if (any(!is.na(tagged))) return(ifelse(is.na(tagged), "other", tagged))
types <- vapply(blocks, function(b) tolower(b$type %||% ""), character(1))
ifelse(types %in% .SPATIAL_NL_TYPES, "spatial",
ifelse(types %in% .TEMPORAL_NL_TYPES, "temporal", "other"))
}
# Columns of a nested fit's latent vector (`modes`, the per-cell precision)
# holding the first block whose role is `role`, or NULL. A multi-block solve
# reports where each block starts (`block_latent_offsets`, from the C++ block
# builders, which are what decide a block's width); a single-block solve lays
# out [fixed, RE, field], so its field is the tail. `blocks` is the fit's block
# list (or single block), `n_x` the latent length.
.nl_role_cols <- function(blocks, offsets, n_x, role) {
roles <- .nl_block_roles(blocks)
k <- which(roles == role)[1L]
if (is.na(k) || is.null(n_x)) return(NULL)
if (length(roles) == 1L && !is.null(blocks$type)) {
w <- .nl_single_block_width(blocks)
if (is.null(w) || w < 1L || w > n_x) return(NULL)
return(as.integer(seq.int(n_x - w + 1L, n_x)))
}
if (length(offsets) != length(roles) + 1L) return(NULL)
lo <- offsets[k]; hi <- offsets[k + 1L]
if (!is.finite(lo) || !is.finite(hi) || hi <= lo || hi > n_x) return(NULL)
as.integer(seq.int(lo + 1L, hi))
}
# Width of a single temporal block's field (the only single-block field whose
# coefficients a read-back slices; the rest report through their own
# accessors). NULL when the block does not say.
.nl_single_block_width <- function(p) {
nt <- p$n_times
if (is.null(nt)) return(NULL)
as.integer(nt) * max(as.integer(p$n_groups %||% 1L), 1L)
}
# Structural spec for one block sent to cpp_nested_laplace_multi. Grid values
# go in theta_grid separately; this function returns only the type-specific
# data needed to build a LatentBlock (idx vector, sizes, adjacency, etc.).
#
# `block_joint_grid` is the joint-grid sub-matrix (n_joint x block$theta_dim)
# already realised from the Cartesian product. tgmrf needs it to precompute
# Q(theta_k) at every joint-grid row; built-in types ignore the argument
# because they read theta values directly from theta_grid in the C++ side.
.nl_block_spec_for_cpp <- function(p, block_joint_grid = NULL, n_obs = NULL,
block_index = NULL) {
type <- tolower(p$type)
projected <- .nl_block_is_projected(p)
.nl_check_block_fields(p, c("axis", if (projected) "projected" else "multi"),
block_index)
# Optional per-observation design weight, attached after the type branch has
# built the spec: it turns any block into a varying coefficient whose eta
# contribution is svc_weight[i] * <the block's field at i>, and the C++ side
# honours it on every contribution kind. Reading it per type would decide
# which fields can vary by where the read happens to be written.
.with_svc <- function(out) {
if (is.null(p$svc_weight)) return(out)
out$svc_weight <- as.numeric(p$svc_weight)
out
}
if (type %in% .NL_PROJECTABLE_AREAL) {
out <- list(
type = type,
spatial_idx = if (!projected) as.integer(p$spatial_idx),
n_spatial_units = as.integer(p$n_spatial_units),
adj_row_ptr = as.integer(p$adj_row_ptr),
adj_col_idx = as.integer(p$adj_col_idx),
n_neighbors = as.integer(p$n_neighbors)
)
if (projected) {
out$projector <- .nl_block_projector(p, 1L, block_index %||% 1L, n_obs)
}
if (type == "bym2") {
out$scale_factor <- as.numeric(p$scale_factor %||% 1.0)
if (!is.null(p$node_prec)) out$node_prec <- as.numeric(p$node_prec)
}
.with_svc(out)
} else if (type == "nngp") {
out <- c(list(type = "nngp",
spatial_idx = if (!projected) as.integer(p$spatial_idx)),
.nl_nngp_field_args(p))
if (projected) {
out$projector <- .nl_block_projector(p, 1L, block_index %||% 1L, n_obs)
}
.with_svc(out)
} else if (type == "hsgp") {
.with_svc(list(
type = "hsgp",
phi_basis = as.matrix(p$phi_basis),
lambda_eig = as.numeric(p$lambda_eig)
))
} else if (type %in% c("rw1", "rw2", "ar1")) {
out <- list(
type = type,
temporal_idx = as.integer(p$temporal_idx),
n_times = as.integer(p$n_times)
)
if (type %in% c("rw1", "rw2")) out$cyclic <- isTRUE(p$cyclic)
.with_svc(out)
} else if (type == "iid") {
.with_svc(list(
type = "iid",
obs_idx = as.integer(p$obs_idx),
n_units = as.integer(p$n_units)
))
} else if (type == "spde") {
# FEM mesh + projection. The C++ side wraps the single A CSC in a length-1
# per-arm list and builds the LatentBlock via make_spde_block. Optional
# rational coefficients (fractional nu) pass through when present.
out <- list(
type = "spde",
n_mesh = as.integer(p$n_mesh),
n_obs = as.integer(p$n_obs),
A_x = as.numeric(p$A_x),
A_i = as.integer(p$A_i),
A_p = as.integer(p$A_p),
C0_diag = as.numeric(p$C0_diag),
G1_x = as.numeric(p$G1_x),
G1_i = as.integer(p$G1_i),
G1_p = as.integer(p$G1_p),
nu = as.numeric(p$nu)
)
if (!is.null(p$rational_poles) && !is.null(p$rational_weights)) {
out$rational_poles <- as.numeric(p$rational_poles)
out$rational_weights <- as.numeric(p$rational_weights)
}
.with_svc(out)
} else if (type == "tgmrf") {
# Precompute Q(theta_k) for every joint-grid row and pack CSC triples +
# logdet + log p(theta_k). The C++ side reads them directly at grid
# point k - no R callback during the inner Newton loop. Iterating over
# block_joint_grid (rather than the block's own grid) is what lets
# tgmrf compose with other blocks in a multi-block prior: when other
# blocks vary across joint rows but tgmrf does not, the duplicated Q
# matrices land at the right joint index.
#
# The Q / log_prior(theta) evaluation either goes through the user's R
# closure (`backend == "r"`, the default) or through the registered C++
# spec (`backend == "cpp"`). Wire format downstream is identical in both
# cases.
if (is.null(block_joint_grid)) {
stop("Internal error: tgmrf block_spec called without block_joint_grid.",
call. = FALSE)
}
tg <- as.matrix(block_joint_grid)
n_cells <- nrow(tg)
n_latent <- p$n_latent
Q_csc_p <- vector("list", n_cells)
Q_csc_i <- vector("list", n_cells)
Q_csc_x <- vector("list", n_cells)
logdet_Q <- numeric(n_cells)
log_pi_th <- numeric(n_cells)
use_cpp <- isTRUE(identical(p$backend, "cpp"))
if (use_cpp) {
if (is.null(p$cpp_id) || !nzchar(p$cpp_id)) {
stop("Internal error: tgmrf block has backend = 'cpp' ",
"but no cpp_id.", call. = FALSE)
}
if (!cpp_tgmrf_registry_has(p$cpp_id)) {
stop("tgmrf C++ spec '", p$cpp_id, "' is not registered in the ",
"current R session. Was the user .cpp file sourceCpp'd? ",
"(Re-load via tgmrf_cpp(cpp_file = ..., id = '", p$cpp_id,
"').)", call. = FALSE)
}
}
for (k in seq_len(n_cells)) {
th <- as.numeric(tg[k, ])
names(th) <- p$theta_names
if (use_cpp) {
evk <- cpp_tgmrf_eval(p$cpp_id, th)
if (as.integer(evk$n) != n_latent) {
stop("Registered tgmrf_cpp spec '", p$cpp_id,
"' returned Q of size ", evk$n, " at grid row ", k,
" but the block was registered with n_latent = ", n_latent,
".", call. = FALSE)
}
Q_csc_p[[k]] <- as.integer(evk$p)
Q_csc_i[[k]] <- as.integer(evk$i)
Q_csc_x[[k]] <- as.numeric(evk$x)
Qk <- Matrix::sparseMatrix(
i = Q_csc_i[[k]], p = Q_csc_p[[k]], x = Q_csc_x[[k]],
dims = c(n_latent, n_latent), index1 = FALSE
)
lp <- as.numeric(evk$log_prior)
} else {
Qk <- p$Q(th)
Qk <- methods::as(methods::as(Qk, "generalMatrix"), "CsparseMatrix")
if (nrow(Qk) != n_latent || ncol(Qk) != n_latent) {
stop("`Q(theta)` returned ", nrow(Qk), " x ", ncol(Qk),
" at grid row ", k,
" but the block was registered with n_latent = ", n_latent, ".",
call. = FALSE)
}
Q_csc_p[[k]] <- as.integer(Qk@p)
Q_csc_i[[k]] <- as.integer(Qk@i)
Q_csc_x[[k]] <- as.numeric(Qk@x)
lp <- as.numeric(p$prior(th))
}
# sparse Cholesky log-determinant; falls back to dense det for tiny
# blocks where Matrix::Cholesky balks.
ld <- tryCatch({
ch <- Matrix::Cholesky(methods::as(Qk, "symmetricMatrix"), LDL = FALSE)
# log|Q| = 2 * sum(log(diag(L))) for the LL' factor.
2 * Matrix::determinant(ch, sqrt = TRUE, logarithm = TRUE)$modulus
}, error = function(e) {
as.numeric(determinant(Qk, logarithm = TRUE)$modulus)
})
logdet_Q[k] <- as.numeric(ld)
if (!is.finite(lp)) {
stop("`prior(theta)` returned a non-finite value at grid row ", k,
" (theta = ", paste(format(th, digits = 4), collapse = ", "),
").", call. = FALSE)
}
log_pi_th[k] <- lp
}
list(
type = "tgmrf",
n_latent = as.integer(n_latent),
obs_idx = as.integer(p$obs_idx),
Q_csc_p_per_grid = Q_csc_p,
Q_csc_i_per_grid = Q_csc_i,
Q_csc_x_per_grid = Q_csc_x,
logdet_Q_per_grid = as.numeric(logdet_Q),
log_prior_theta_per_grid = as.numeric(log_pi_th)
)
} else {
stop("Block type '", type, "' is not supported in multi-block priors.",
call. = FALSE)
}
}
# Per-block axis grid as a matrix `[n_block_cells x n_axes_for_block]`.
# Single-axis types (icar / rw1 / rw2 / iid) get a 1-column matrix; 2-axis
# types (bym2 / car_proper / ar1) get a 2-column matrix.
.nl_block_axis_grid <- function(p) {
spec <- .NL_REGISTRY[[tolower(p$type)]]
if (is.null(spec)) {
stop("Unknown block type: ", p$type, call. = FALSE)
}
.nl_check_block_fields(p, "axis")
p <- spec$defaults(p, list())
th <- spec$theta(p)
if (is.matrix(th$grid)) {
grid <- th$grid
if (is.null(colnames(grid))) colnames(grid) <- th$names
} else {
grid <- matrix(as.numeric(th$grid), ncol = 1)
colnames(grid) <- th$names
}
list(grid = grid, names = colnames(grid), prepared = p)
}
# Wall-clock threshold, in seconds, above which a multi-block outer-grid
# solve gets a "this was slow" warning. Cell count alone does not predict
# solve time -- a family/likelihood pair, N, or a block's own per-cell cost
# (gcol33/tulpa#638) can make the same cell count cheap or expensive, and the
# engine's own default per-block grids already exceed any static cell-count
# threshold that fires on genuinely slow custom grids (gcol33/tulpa#820) --
# so the trigger is measured cost, read off the "grid" bucket both multi-block
# dispatchers already mark with `.tulpa_timer()`.
.NL_MULTI_GRID_WARN_SECONDS <- 2
# Shared by the single-response multi-block dispatch (`tulpa_nested_laplace()`
# in this file) and the joint multi-block dispatch
# (`R/nested_laplace_joint_multi.R`): both solve their outer grid in one call
# and mark its elapsed time as the "grid" timing bucket, so both warn off that
# measurement through this one function rather than keeping two copies of the
# same message. `layout` is `.nl_grid_crossing()`'s account of the axes that
# produced the count, so the advice names which grids to shrink.
.nl_multi_grid_warn <- function(elapsed, n_cells, remedy, layout = "") {
if (isTRUE(elapsed > .NL_MULTI_GRID_WARN_SECONDS) && !.nl_internal_batch()) {
warning(sprintf(
"Multi-block outer grid (%s cells) took %s to solve.%s %s",
format(n_cells, scientific = FALSE, trim = TRUE),
.format_duration(elapsed), layout, remedy
), call. = FALSE)
}
invisible(NULL)
}
# Cell count above which a dense multi-block outer tensor is announced with a
# warning BEFORE its solve starts. Each cell is one inner Laplace solve, so the
# count is the run's size in solves; it does not know the per-solve cost, which
# is why it warns rather than refuses and why the post-solve warning above reads
# measured time.
.NL_DENSE_GRID_WARN_CELLS <- 2048L
# Shared by both multi-block dispatchers: the dense tensor about to be solved
# has `n_cells` inner solves, `layout` (`.nl_grid_crossing()`) names the axes
# that produced them and `remedy` is the dispatcher's advice. Gated like the
# timing warning, so an internal batch the caller did not size stays quiet.
.nl_dense_grid_warn <- function(n_cells, remedy, layout = "") {
if (isTRUE(n_cells > .NL_DENSE_GRID_WARN_CELLS) && !.nl_internal_batch()) {
warning(sprintf(
paste0("Multi-block outer grid is a dense tensor of %s cells, one inner ",
"solve each.%s %s"),
format(n_cells, scientific = FALSE, trim = TRUE), layout, remedy
), call. = FALSE)
}
invisible(NULL)
}
# Whether the dispatch currently running is an INTERNAL batch rather than the
# grid the caller asked for. The soft-cap warning above is advice addressed to
# whoever chose the cell count, and on an internal batch that is not the caller:
# the outer Pareto-k diagnostic re-evaluates `log_marginal` at `k_samples`
# importance draws by substituting them for the block's grid axis, so a 7-node
# fit warned about a cell count the diagnostic's own budget set and advised
# reducing per-block grid sizes, which does not reach the number in the message
# (gcol33/tulpa#614).
#
# Transported on an option rather than an argument because the two warning sites
# sit behind different dispatchers -- `.nl_dispatch_multi()` takes `cargs`, the
# joint multi-block driver does not -- which is the same reason
# `tulpa.nl_progress` and `tulpa.nl_checkpoint` are options.
.nl_internal_batch <- function() {
isTRUE(getOption("tulpa.nl_internal_batch", FALSE))
}
# Run `expr` with that flag set. Nests and restores, so an internal batch that
# itself re-dispatches stays internal.
.nl_with_internal_batch <- function(expr) {
op <- options(tulpa.nl_internal_batch = TRUE)
on.exit(options(op), add = TRUE)
force(expr)
}
# The caller's ceiling on a dense multi-block outer tensor, in inner solves.
# Unset there is none: a cell count does not measure the run's cost, so the
# default is the pre-solve warning above, and a caller who wants a hard stop
# (a batch of fits on a shared node, say) sets `control$max_grid_cells`. The
# front doors publish the resolved value on the scoped
# `tulpa.nl_max_grid_cells` option, so every grid built inside the fit -- the
# multi-block dispatch, the joint multi-block dispatch, the k-hat
# re-evaluations -- reads one ceiling without carrying it in each signature,
# the transport `tulpa.nl_progress` and `tulpa.nl_checkpoint` already use.
.nl_max_grid_cells <- function(control = list()) {
v <- control$max_grid_cells
if (is.null(v)) v <- getOption("tulpa.nl_max_grid_cells", NULL)
if (is.null(v)) return(Inf)
if (!is.numeric(v) || length(v) != 1L || is.na(v) || v < 1) {
stop("`control$max_grid_cells` must be a single number >= 1 ",
"(the cell-count ceiling on a multi-block outer grid).", call. = FALSE)
}
as.numeric(v)
}
# The single enforcement of that ceiling, shared by both multi-block
# dispatchers, before any inner solve. `layout` is `.nl_grid_crossing()`'s
# account of the axes that produced the count and `remedy` the dispatcher's
# advice.
.nl_check_grid_cap <- function(n_cells, max_cells, remedy, layout = "") {
if (n_cells <= max_cells) return(invisible(n_cells))
fmt <- function(x) format(x, scientific = FALSE, trim = TRUE)
stop(sprintf(
paste0("Multi-block outer grid has %s cells, past control$max_grid_cells ",
"= %s.%s %s"),
fmt(n_cells), fmt(max_cells), layout, remedy
), call. = FALSE)
}
# The axes a dense outer tensor crosses, as the sentence the timing warning
# carries: " It crosses b1 (...) x phi_pos 4." `block_grids` are the per-block
# axis grids, so a fit that set no grid learns which default axes produced its
# count (gcol33/tulpa#913); `phi_axes` are the active named per-arm dispersion
# axes crossed on top, each a factor on the number of inner solves.
.nl_grid_crossing <- function(block_grids = NULL, phi_axes = NULL) {
parts <- c(if (length(block_grids)) .nl_grid_layout(block_grids),
if (length(phi_axes))
paste0("phi_", names(phi_axes), " ", lengths(phi_axes)))
if (length(parts)) paste0(" It crosses ", paste(parts, collapse = " x "), ".")
else ""
}
# One block's grid as "b<k> (<rows> rows: <axis> <levels> x ...)", the levels
# being each axis's distinct node count. A block grid is a tensor of its axes
# for the grid-axis families and a paired node list for others (BYM2's
# (sigma, rho) nodes), so the row count is stated alongside the per-axis
# levels rather than derived from them.
.nl_grid_layout <- function(block_grids) {
parts <- vapply(seq_along(block_grids), function(b) {
g <- as.matrix(block_grids[[b]])
nm <- colnames(g) %||% paste0("axis", seq_len(ncol(g)))
lv <- vapply(seq_len(ncol(g)), function(j)
length(unique(as.numeric(g[, j]))), integer(1))
sprintf("b%d (%d rows: %s)", b, nrow(g),
paste(nm, lv, collapse = " x "))
}, character(1))
paste(parts, collapse = " x ")
}
.nl_dispatch_multi <- function(cargs, prior_list, likelihood = NULL,
progress = .nl_progress_args(list(progress = FALSE)),
within_cell = NULL,
# Single-cell re-dispatch for the inner-Laplace
# skewness diagnostic (.nl_inner_skew_at_theta()):
# when supplied, `theta_grid_override` (a 1-row
# matrix at the outer grid's MAP cell) replaces
# the Cartesian-product grid construction below,
# reusing the SAME blocks_spec machinery so the
# re-solve is at exactly the fitted MAP theta.
# The SPDE-prior fold / weight normalisation /
# posterior-moments tail is skipped in that case
# -- the caller only wants inner_skew, not a
# re-normalised single-cell "posterior".
theta_grid_override = NULL,
compute_skew = FALSE, skew_idx = NULL,
debias = NULL,
# Corrected integrated Laplace:
# the kernel-facing request list, off by default.
cila = NULL) {
# Inject a default obs_idx for any tgmrf block that didn't supply one.
# The C++ scatter needs obs_idx[i] -> latent slot for each observation;
# the canonical "one obs per latent slot" case has N == n_latent and the
# identity mapping. Heterogeneous (N != n_latent) layouts must supply
# obs_idx explicitly at tgmrf() construction.
N_obs <- length(cargs$y)
for (b in seq_along(prior_list)) {
blk <- prior_list[[b]]
if (identical(tolower(blk$type %||% ""), "tgmrf") && is.null(blk$obs_idx)) {
if (blk$n_latent != N_obs) {
stop("tgmrf block has n_latent = ", blk$n_latent,
" but data has N = ", N_obs, " observations. Provide an explicit ",
"`obs_idx` argument to tgmrf() mapping observations to latent slots.",
call. = FALSE)
}
prior_list[[b]]$obs_idx <- seq_len(N_obs)
}
}
per_block <- lapply(prior_list, .nl_block_axis_grid)
prepared <- lapply(per_block, function(x) x$prepared)
block_grids <- lapply(per_block, function(x) x$grid)
axis_counts <- vapply(block_grids, ncol, integer(1))
axis_offsets <- as.integer(c(0L, cumsum(axis_counts)))
# Block-prefix the axis names so they're unambiguous (e.g. block1.tau).
axis_names <- unlist(lapply(seq_along(block_grids), function(b) {
paste0("b", b, ".", colnames(block_grids[[b]]))
}))
# Set only on the Cartesian-product branch below, so the post-solve timing
# check (which reads these) is a no-op on a single-cell override re-dispatch
# (skew / k-hat / recenter probes -- already never the caller's own grid).
n_cells <- NA_integer_
grid_warn_remedy <- NULL
grid_layout <- ""
if (!is.null(theta_grid_override)) {
joint_grid <- as.matrix(theta_grid_override)
total_axes <- axis_offsets[length(axis_offsets)]
if (ncol(joint_grid) != total_axes) {
stop("theta_grid_override must have ", total_axes,
" columns (sum of block axes), got ", ncol(joint_grid), ".",
call. = FALSE)
}
colnames(joint_grid) <- axis_names
} else {
# Cartesian product of per-block row indices. This path has no CCD
# alternative (unlike the joint multi-block dispatch), so the only
# remedy for a slow or oversized grid is a smaller one.
row_counts <- vapply(block_grids, nrow, integer(1))
idx <- do.call(expand.grid, lapply(row_counts, seq_len))
n_cells <- nrow(idx)
grid_warn_remedy <- "Reduce per-block grid sizes."
grid_layout <- .nl_grid_crossing(block_grids)
.nl_check_grid_cap(n_cells, .nl_max_grid_cells(), grid_warn_remedy,
grid_layout)
.nl_dense_grid_warn(n_cells, grid_warn_remedy, grid_layout)
# Concatenate per-block axis grids into the joint theta_grid.
joint_grid <- do.call(cbind, lapply(seq_along(block_grids), function(b) {
block_grids[[b]][idx[[b]], , drop = FALSE]
}))
colnames(joint_grid) <- axis_names
}
blocks_spec <- lapply(seq_along(prepared), function(b) {
cols <- (axis_offsets[b] + 1L):axis_offsets[b + 1L]
block_joint <- joint_grid[, cols, drop = FALSE]
.nl_block_spec_for_cpp(prepared[[b]], block_joint,
n_obs = length(cargs$y), block_index = b)
})
# Every block's hyperprior over its own columns, folded on the override path
# too: a placement stencil re-evaluates the marginal through it and has to
# difference the target the grid integrates. A screened solve ranks with it.
# The integrated columns are read off the block's declared grid, never off the
# override, whose rows hold a fixed axis at other values or an integrated one
# at a single value (gcol33/tulpa#760).
hp_parts <- lapply(seq_along(prepared), function(b) {
cols <- (axis_offsets[b] + 1L):axis_offsets[b + 1L]
tg_b <- joint_grid[, cols, drop = FALSE]
colnames(tg_b) <- colnames(block_grids[[b]])
blk <- prepared[[b]]
blk$log_prior_theta_per_grid <- blocks_spec[[b]]$log_prior_theta_per_grid
.hp_prefix(.nl_block_log_hyperprior(blk, tg_b, .hp_choice(cargs$hyperprior),
declared = .hp_declare(block_grids[[b]])),
paste0("b", b, "."))
})
prune_tol <- as.numeric(cargs$prune_tol %||% 0)
grid_solve_start <- proc.time()[["elapsed"]]
out <- cpp_nested_laplace_multi(
y = cargs$y,
n = cargs$n,
X = cargs$X,
re_idx = cargs$re_idx,
n_re_groups = cargs$n_re_groups,
sigma_re = cargs$sigma_re,
blocks_spec = blocks_spec,
theta_grid = joint_grid,
axis_offsets = axis_offsets,
family = cargs$family,
phi = cargs$phi,
max_iter = cargs$max_iter,
tol = cargs$tol,
n_threads = cargs$n_threads,
x_init_nullable = cargs$x_init_nullable,
store_Q = isTRUE(cargs$store_Q),
# Cheap-pass screening tolerance; 0 (the default) solves every cell. A
# single-cell re-dispatch inherits it harmlessly -- the driver screens only
# a grid with more than one cell.
prune_tol = prune_tol,
likelihood = likelihood,
progress = isTRUE(progress$progress),
progress_every = as.integer(progress$progress_every),
progress_throttle = as.numeric(progress$progress_throttle),
progress_file = as.character(progress$progress_file),
checkpoint_path = as.character(cargs$checkpoint_path %||% ""),
compute_skew = isTRUE(compute_skew),
skew_idx = skew_idx,
debias = debias,
cila = cila,
offset_nullable = cargs$offset_nullable,
screen_log_offset = if (prune_tol > 0)
.nl_screen_log_offset(joint_grid, hp_parts,
screen = cargs$screen_refine)
)
if (!is.null(grid_warn_remedy)) {
.nl_multi_grid_warn(proc.time()[["elapsed"]] - grid_solve_start,
n_cells, grid_warn_remedy, grid_layout)
}
out$theta_grid <- joint_grid
out$theta_names <- axis_names
out$axis_offsets <- axis_offsets
out$blocks <- prepared
out <- .nl_fold_hyperprior(out, hp_parts)
if (!is.null(theta_grid_override)) return(out) # skew re-dispatch / the
# auto-recenter FD stencil
out <- .nl_attach_outer_integration(out, "multi-block outer grid")
out <- .nl_posterior_moments_multi(out, prepared, axis_offsets, joint_grid,
within = .nl_within_cell_mode(within_cell))
out
}
#' Build a `prior` list for [tulpa_nested_laplace()] from a tulpa spec object
#'
#' @description
#' Validates a `tulpa_temporal` or `tulpa_spatial` specification against
#' `data`, then converts it to the prior list shape consumed by
#' [tulpa_nested_laplace()]. Mainly an internal helper for callers that already
#' have a fitted spec; users typically pass `spec` + `data` directly to
#' `tulpa_nested_laplace()` instead.
#'
#' Supported spec types:
#' * `tulpa_temporal` with `type in {"rw1", "rw2", "ar1"}`
#' * `tulpa_spatial`, areal `type in {"car", "icar", "car_proper", "bym2"}`.
#' An areal spec returns the graph fields only -- `type`, `spatial_idx`,
#' `n_spatial_units`, `adj_row_ptr`, `adj_col_idx`, `n_neighbors` -- plus
#' `scale_factor` for bym2 and the eigenvalue-derived `rho_bounds` for
#' proper CAR. No grid field is returned for any areal type: the outer axes
#' are built later by the family's registry `defaults()` closure, and proper
#' CAR's `(tau, rho)` grid is built there from those bounds.
#' * `tulpa_gp` / `tulpa_hsgp`, continuous `type in {"gp", "nngp", "hsgp"}`:
#' validated against `data` and routed to the `nngp` / `hsgp` nested kernel
#' through the shared converter.
#'
#' SPDE is the one continuous field not built here -- it carries its own
#' (range, sigma) FEM integrator; call [fit_spde()] with the `spatial_spde()`
#' spec.
#'
#' @param spec A `tulpa_temporal` or `tulpa_spatial` object.
#' @param data Data frame the spec resolves time/group/site indices against.
#' @return A `prior` list ready for [tulpa_nested_laplace()].
#' @examples
#' d <- data.frame(year = rep(1:10, each = 3))
#' pr <- prior_from_spec(temporal_ar1("year"), d)
#' str(pr)
#' @export
prior_from_spec <- function(spec, data) {
if (inherits(spec, "tulpa_temporal")) {
return(.prior_from_temporal_spec(spec, data))
}
if (inherits(spec, "tulpa_spatial")) {
return(.prior_from_spatial_spec(spec, data))
}
stop("`spec` must inherit from 'tulpa_temporal' or 'tulpa_spatial'.",
call. = FALSE)
}
.prior_from_temporal_spec <- function(spec, data) {
spec <- validate_temporal(spec, data)
type <- tolower(spec$type)
if (!type %in% c("rw1", "rw2", "ar1")) {
stop("tulpa_nested_laplace() supports temporal types rw1, rw2, ar1; got '",
type, "'.", call. = FALSE)
}
out <- list(
type = type,
temporal_idx = as.integer(spec$time_index),
n_times = as.integer(spec$n_times)
)
if (type %in% c("rw1", "rw2")) out$cyclic <- isTRUE(spec$cyclic)
out
}
.prior_from_spatial_spec <- function(spec, data) {
type <- tolower(spec$type)
# Continuous coordinate-addressed fields (gp/nngp/hsgp). Validate to populate
# the spec's resolved coordinate structure (NNGP neighbour info / HSGP basis),
# then build the nested prior block through the shared converter
# `.spatial_spec_to_nl_prior` -- the single source of truth for spec -> nested
# prior, also used by the tulpa() front door. gp/nngp drive the `nngp` kernel,
# hsgp the `hsgp` kernel; both are single-block (cpp_fn present in .NL_REGISTRY).
if (type %in% c("gp", "nngp")) {
return(.spatial_spec_to_nl_prior(validate_gp(spec, data)))
}
if (type == "hsgp") {
return(.spatial_spec_to_nl_prior(validate_hsgp(spec, data)))
}
if (type == "spde") {
stop("tulpa_nested_laplace() does not host an SPDE field: SPDE carries its ",
"own (range, sigma) FEM integrator. Call fit_spde() with the ",
"spatial_spde() spec.", call. = FALSE)
}
if (!type %in% .NL_FRONTDOOR_AREAL) {
stop("Unknown spatial type '", type, "'. Areal: ",
paste(.NL_FRONTDOOR_AREAL, collapse = "/"),
"; continuous: gp/nngp/hsgp; SPDE via fit_spde().", call. = FALSE)
}
validate_spatial(spec, data)
# Areal field: resolve the per-observation spatial_idx (1-based row indices
# into the adjacency) onto the spec, then build the block through the shared
# converter. Empty cells (no data) are allowed -- the ICAR / CAR / BYM2 prior
# is well-defined on every graph node; unobserved cells contribute no
# likelihood term. See `.resolve_spatial_idx()` for the resolution rules.
adj <- spec$adjacency
n_spatial_units <- nrow(as.matrix(adj))
if (!is.null(spec$level) && spec$level == "group" &&
!is.null(spec$group_var)) {
if (!(spec$group_var %in% names(data))) {
stop("Spatial group variable '", spec$group_var,
"' not found in data.", call. = FALSE)
}
spec$spatial_idx <- .resolve_spatial_idx(
values = data[[spec$group_var]],
n_spatial_units = n_spatial_units,
adjacency = adj,
group_var = spec$group_var
)
} else {
if (nrow(data) != n_spatial_units) {
stop("Observation-level spatial spec requires nrow(data) == ",
"nrow(adjacency).", call. = FALSE)
}
spec$spatial_idx <- seq_len(n_spatial_units)
}
.spatial_spec_to_nl_prior(spec)
}
# Normalize the outer-grid progress knobs from a `control` list into the four
# scalars the C++ entry points accept. On by default; `progress = FALSE` turns
# off the flushed cell-k/n_grid + ETA reporter (see the GridProgress reporter in
# inst/include/tulpa/nested_progress.h). `progress.file`
# adds a heartbeat file for detached runs where Rcout flushing is unreliable.
# The inner refinement / EM / CCD-probe call sites pass `progress = FALSE`
# explicitly so only the top-level fit ticks.
.nl_progress_args <- function(control) {
# Control keys win; otherwise fall back to the scoped `tulpa.nl_progress`
# option a caller (e.g. tulpaObs) may have set for a whole fit, so progress
# reaches grids run through fixed-control internal paths (the EM per-block
# nested solve) without threading the knobs through every layer.
# `[[` (exact) not `$` on `progress`: it is a strict prefix of the other three
# keys, so a `$` read returns whichever ONE of them is set (and NULL when two
# are, the match being ambiguous). `isTRUE()` of a file path or a cadence is
# FALSE, which turned asking for the heartbeat file into switching the console
# bar off.
has_ctrl <- !is.null(control[["progress"]]) || !is.null(control$progress.every) ||
!is.null(control$progress.throttle) || !is.null(control$progress.file)
if (!has_ctrl) {
opt <- getOption("tulpa.nl_progress", NULL)
if (is.list(opt)) return(opt)
}
list(
progress = isTRUE(control[["progress"]] %||% TRUE),
progress_every = as.integer(control$progress.every %||% 0L),
progress_throttle = as.numeric(control$progress.throttle %||% 2),
progress_file = as.character(control$progress.file %||% "")
)
}
# Resolve a checkpoint spec to a normalized `list(path, resume)`.
# `control$checkpoint = list(path = , resume = TRUE)` enables it; `resume =
# FALSE` starts over (the caller removes any prior file before the first
# kernel call). Like `.nl_progress_args`, an absent control key falls back to
# the scoped `tulpa.nl_checkpoint` option (when `use_option`, the default) so
# the checkpoint reaches grids run through fixed-control internal paths;
# `tulpa_sample_glmm()`'s per-chain checkpoint is a distinct namespace from
# the nested-Laplace grid-cell one this option serves, so it reads
# `use_option = FALSE` and ignores it.
.nl_checkpoint_args <- function(control, use_option = TRUE) {
cp <- control$checkpoint
if (is.null(cp)) {
if (use_option) {
opt <- getOption("tulpa.nl_checkpoint", NULL)
if (is.list(opt)) return(opt)
}
return(list(path = "", resume = TRUE))
}
if (!is.list(cp) || is.null(cp$path) || !is.character(cp$path) ||
length(cp$path) != 1L || is.na(cp$path) || !nzchar(cp$path)) {
stop("`control$checkpoint` must be a list with a single non-empty ",
"character `path` (e.g. ",
"list(path = \"fit.ckpt\", resume = TRUE)).", call. = FALSE)
}
resume <- cp$resume %||% TRUE
if (!is.logical(resume) || length(resume) != 1L || is.na(resume)) {
stop("`control$checkpoint$resume` must be a single TRUE or FALSE.",
call. = FALSE)
}
list(path = path.expand(cp$path), resume = isTRUE(resume))
}
Try the tulpa package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.