R/inference_modes.R

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

#' Inference mode system
#'
#' @description
#' tulpa uses an explicit tier system for inference that encodes
#' **epistemic guarantees**, not just runtime characteristics.
#'
#' This design makes the difference between inference methods **first-class
#' and unavoidable**, rather than hiding them as implementation details.
#'
#' @section The Three Tiers:
#'
#' **Tier 1 - Exact:**
#' Asymptotically correct posterior inference (up to Monte Carlo error).
#' Credible intervals are interpretable as posterior uncertainty.
#' This is the reference standard.
#'
#' **Tier 2 - Structured:**
#' Accurate inference *conditional on explicit structural assumptions*.
#' Typically requires latent Gaussian structure, conditional independence,
#' smooth posteriors. Very fast for the right model class, but can be wrong
#' outside that class. Failure modes are predictable and explainable.
#'
#' **Tier 3 - Optimized:**
#' No general correctness guarantee beyond empirical usefulness.
#' Point estimates usually good, but uncertainty often underestimated.
#' Tails and correlations unreliable. Failure is usually silent.
#' This is optimization, not sampling.
#'
#' @section Auto Mode:
#'
#' `mode = "auto"` chooses between Tier 1 and Tier 2 only.
#' It will **never** silently choose Tier 3 (Optimized).
#'
#' The contract: "Use the most reliable method that is expected to finish
#' for this model."
#'
#' Auto decisions are deterministic, explainable, and overrideable.
#'
#' @section Implementation Rules:
#'
#' 1. Modes change semantics, not just runtime.
#'    Intervals from Optimized do not mean the same as Exact.
#'
#' 2. The mode must always be visible in output.
#'
#' 3. No silent upgrading or downgrading.
#'    If Exact fails, we error - we do not switch to Structured.
#'
#' 4. Backends slot into tiers.
#'    Adding a backend never introduces a new epistemic promise.
#'
#' @name inference_modes
#' @noRd
NULL


# ==============================================================================
# Tier metadata + backend registry (single source of truth)
# ==============================================================================

#' Tier metadata: the epistemic guarantee attached to each tier
#'
#' Tier *membership* (which backends live in a tier) is not stored here -- it
#' is derived from `BACKEND_REGISTRY`. This table holds only the per-tier
#' description/guarantee so the two never drift.
#'
#' @noRd
TIER_META <- list(
  exact = list(
    tier = 1L,
    name = "Exact",
    description = "Asymptotically correct posterior inference",
    guarantee = "Credible intervals interpretable as posterior uncertainty",
    note = "Reference standard"
  ),
  structured = list(
    tier = 2L,
    name = "Structured",
    description = "Accurate inference conditional on structural assumptions",
    guarantee = "Correct if model meets structural assumptions",
    note = "Controlled approximation, not heuristics"
  ),
  optimized = list(
    tier = 3L,
    name = "Optimized",
    description = "Optimization-based approximate inference",
    guarantee = "No general correctness guarantee",
    note = "Uncertainty often underestimated; requires explicit opt-in"
  )
)


#' Backend registry -- single source of truth for the inference backends
#'
#' @description
#' One entry per backend. Adding a backend is a single entry here; tier
#' membership, family support, the input contract, and R-level reachability
#' all derive from this list -- the single source of truth for per-tier
#' `backends` and family support.
#'
#' Fields:
#' * `emits`    -- the kind of posterior representation the backend's draws
#'     hold, which fixes whether MCMC chain diagnostics apply (independent of
#'     `tier`: an exact SMC sampler emits `"iid"` particles, a Tier-3 VI fit
#'     also emits `"iid"`, while a nested-Laplace Tier-2 fit emits `"iid"` too):
#'     - `"chain"` : autocorrelated MCMC output -- Rhat and autocorrelation-ESS
#'                   are meaningful (`diagnostics()` computes them).
#'     - `"iid"`   : exchangeable draws from a deterministic approximation, or
#'                   resampled particles -- split-Rhat is vacuous and ESS = n by
#'                   construction, so they say nothing about approximation bias.
#'     - `"point"` : no posterior sample (mode + covariance only).
#' * `tier`     -- tier key (`"exact"`, `"structured"`, `"optimized"`).
#' * `input`    -- the input contract the backend's fitter consumes:
#'     - `"design"`    : design-matrix bundle (`y`, `n_trials`, `X`, RE structure).
#'     - `"logpost"`   : a `log_posterior(theta)` closure plus dimension/init.
#'     - `"modeldata"` : a tulpa `ModelData` / `LikelihoodSpec` (C-ABI NUTS path).
#'     - `"nested"`    : a design-matrix bundle plus one or more latent prior
#'                       blocks (`latent(tgmrf(...))`), integrated over the block
#'                       hyperparameters by the nested-Laplace driver.
#'     - `"spde"`      : a design bundle plus a self-contained SPDE spec
#'                       (mesh + FEM matrices); fit_spde()'s own CCD / grid
#'                       hyperparameter engine.
#' * `fitter`   -- name of the R function implementing the backend, or `NULL`
#'     when only a C++ kernel exists with no R entry point yet. Stored as a
#'     *string* (not the function object) so the registry is independent of
#'     source-load order; resolved lazily via [resolve_backend_fitter()].
#'     `NULL` => not selectable from R; dispatch fails loudly.
#' * `families` -- character vector of supported family identifiers, or `NULL`
#'     for an unrestricted backend.
#' * `cabi`     -- the registered C-ABI callable backing the backend (the symbol
#'     a model package reaches via `LinkingTo: tulpa`, and the one an R wrapper
#'     would call), or `NULL`.
#' * `carries_weights` -- `TRUE` when a `"design"` backend threads observation
#'     `weights` into its likelihood (non-spatially); absent otherwise. Every
#'     `"logpost"` backend carries them through the log-posterior builder, and
#'     `.backend_carries_weights()` is the one predicate reading both.
#' * `hyperprior` -- `TRUE` when the backend's fitter takes the outer
#'     `hyperprior = c("proper", "flat")` choice `tulpa()` forwards; absent
#'     otherwise.
#' * `note`     -- optional human-readable note.
#'
#' Family identity (for `families`) is checked against `family$name`,
#' `family$distribution`, and `family$numerator$distribution` (tulpaRatio-style
#' ratio families nest a per-process distribution).
#'
#' @noRd
BACKEND_REGISTRY <- list(
  # ---- Tier 1: Exact ----
  hmc = list(
    emits = "chain",
    tier = "exact", input = "modeldata", fitter = "tulpa_sample_glmm",
    families = NULL,
    # tulpa_sample_glmm() always runs run_hmc_parallel_chains_cpp (n_chains
    # defaults to 4, but the call shape is the same at n_chains = 1), which is
    # what "tulpa_run_nuts_chains" is registered against; "tulpa_run_nuts_generic"
    # is the single-chain callable a model package reaches on its own and this
    # R door never calls (gcol33/tulpa#775).
    cabi = "tulpa_run_nuts_chains",
    note = paste("Generic NUTS over a fixed-effect GLM (tulpa_sample_glmm);",
                 "model packages also drive NUTS through the C ABI")
  ),
  ess = list(
    emits = "chain",
    tier = "exact", input = "modeldata", fitter = "tulpa_sample_glmm",
    families = NULL, cabi = "tulpa_run_ess_sampler"
  ),
  gibbs = list(
    emits = "chain",
    # The Polya-Gamma spatial sweep updates ONE random-intercept block
    # alongside the field, so a second (1 | g) term is refused at dispatch and
    # auto must not pick it for such a model (gcol33/tulpa#681).
    max_re_terms = 1L,
    tier = "exact", input = "design", fitter = "tulpa_gibbs",
    families = c("binomial", "neg_binomial_2"),
    carries_offset = FALSE,
    cabi = c("tulpa_pg_binomial_gibbs", "tulpa_pg_negbin_gibbs",
             "tulpa_pg_negbin_spatial_gibbs"),
    note = paste("Polya-Gamma (PG) Gibbs (binomial / negbin) via tulpa_gibbs();",
                 "base plus spatial (icar/bym2/rsr/gp/multiscale_gp) and temporal",
                 "dispatch through the cpp_pg_binomial_gibbs* / cpp_pg_negbin_gibbs",
                 "kernels, also reachable from model packages through the C ABI",
                 "symbols above")
  ),
  re_cov_gibbs = list(
    emits = "chain",
    tier = "exact", input = "design", fitter = "tulpa_re_cov_gibbs",
    # The four restrictions its own dispatch enforces, declared here so the auto
    # selector can see them: auto sent every RE-bearing formula straight to this
    # backend, which then refused the call (gcol33/tulpa#666).
    families = c("binomial", "poisson", "gaussian", "neg_binomial_2"),
    carries_offset = FALSE, cabi = NULL,
    note = paste("Correlated random-slope term (1 + x | g): exact",
                 "Metropolis-within-Gibbs debias of the RE covariance Sigma.",
                 "Auto-selected from the Laplace path with control$re_cov = 'gibbs'")
  ),
  # SGHMC / SGLD carry discretization + minibatch bias at a fixed step size, and
  # MCLMC's default is the UNADJUSTED sampler (mclmc_adjusted = 0), so none of
  # the three targets the exact posterior: they are approximate MCMC, tier
  # "optimized" (no general correctness guarantee, explicit opt-in). They still
  # emit chains, so Rhat/ESS measure mixing -- but mixing does not certify the
  # (biased) stationary distribution, which is exactly what "exact" would imply.
  sghmc = list(
    emits = "chain",
    tier = "optimized", input = "modeldata", fitter = "tulpa_sample_glmm",
    families = NULL, cabi = "tulpa_sghmc_fit"
  ),
  sgld = list(
    emits = "chain",
    tier = "optimized", input = "modeldata", fitter = "tulpa_sample_glmm",
    families = NULL, cabi = "tulpa_sgld_fit"
  ),
  mclmc = list(
    emits = "chain",
    tier = "optimized", input = "modeldata", fitter = "tulpa_sample_glmm",
    families = NULL, cabi = "tulpa_mclmc_fit",
    note = "Microcanonical Langevin Monte Carlo over a fixed-effect GLM (unadjusted by default; approximate)"
  ),
  smc = list(
    emits = "iid",
    tier = "exact", input = "modeldata", fitter = "tulpa_sample_glmm",
    families = NULL, cabi = "tulpa_smc_fit",
    note = "Sequential Monte Carlo over a fixed-effect GLM"
  ),
  imh_laplace = list(
    emits = "chain",
    tier = "exact", input = "logpost", fitter = "imh_laplace",
    families = NULL, cabi = NULL
  ),
  mala = list(
    emits = "chain",
    tier = "exact", input = "logpost", fitter = "mala",
    families = NULL, cabi = NULL
  ),
  # ---- Tier 2: Structured ----
  laplace = list(
    emits = "iid",
    tier = "structured", input = "design", fitter = "tulpa_laplace",
    families = NULL, cabi = "tulpa_laplace_spec_dense", carries_weights = TRUE
  ),
  re_cov_nested = list(
    emits = "iid",
    tier = "structured", input = "design", fitter = "tulpa_re_cov_nested",
    families = NULL, cabi = NULL, hyperprior = TRUE,
    # Its inner solve is tulpa_laplace()'s non-spatial kernel, which scales each
    # row's log-likelihood by its weight, so a weighted (1 | g) model keeps its
    # scale INTEGRATED rather than being routed to a conditional backend
    # (gcol33/tulpa#874).
    carries_weights = TRUE,
    note = paste("Correlated random-slope term (1 + x | g): nested-Laplace",
                 "integration over the RE covariance Sigma (CCD design + PC/LKJ",
                 "prior). Auto-selected from the Laplace path for a single",
                 "correlated term")
  ),
  pathfinder = list(
    emits = "iid",
    tier = "structured", input = "logpost", fitter = "pathfinder",
    families = NULL, cabi = NULL
  ),
  agq = list(
    emits = "iid",
    tier = "structured", input = "design", fitter = "agq_fit",
    families = c("binomial", "poisson", "gaussian"), cabi = NULL,
    carries_offset = TRUE,
    note = "Intercept-only adaptive Gauss-Hermite (the shared GLMM oracle)"
  ),
  ep = list(
    emits = "iid",
    tier = "structured", input = "design", fitter = "ep_fit",
    families = NULL, cabi = NULL,
    carries_offset = TRUE, phi2 = TRUE,
    note = paste("Expectation Propagation over a fixed-effect GLM with a",
                 "mean-zero Gaussian coefficient prior; matches marginal moments",
                 "(exact for a Gaussian likelihood), no random effects / fields")
  ),
  eb = list(
    emits = "iid",
    tier = "structured", input = "design", fitter = "tulpa_eb",
    families = NULL, cabi = NULL, hyperprior = TRUE,
    note = paste("Empirical Bayes over the random-effect covariances: the mode",
                 "of the same outer objective re_cov_nested integrates, with",
                 "the fixed effects reported conditional on it. Tier 2 because",
                 "the plug-in drops the hyperparameter uncertainty -- the",
                 "intervals are conditional on Sigma_hat, not marginal over",
                 "Sigma -- so it is opt-in by name and auto never selects it.")
  ),
  nested_laplace = list(
    emits = "iid",
    tier = "structured", input = "nested", fitter = "tulpa_nested_laplace",
    families = NULL,
    # cpp_nested_laplace_multi is an Rcpp export (an R-level entry point), not
    # a registered C-ABI callable -- a model package reaches this driver's
    # per-field kernels individually instead, through the six symbols below
    # (gcol33/tulpa#775).
    cabi = c("tulpa_nested_laplace_icar", "tulpa_nested_laplace_bym2",
             "tulpa_nested_laplace_car_proper", "tulpa_nested_laplace_temporal",
             "tulpa_nested_laplace_nngp", "tulpa_nested_laplace_hsgp"),
    hyperprior = TRUE,
    note = "Single-arm nested Laplace; integrates latent-block hyperparameters"
  ),
  nested_laplace_joint = list(
    emits = "iid",
    tier = "structured", input = "nested", fitter = "tulpa_nested_laplace_joint",
    families = NULL,
    # cpp_nested_laplace_joint_multi is an Rcpp export, not a registered C-ABI
    # callable, and no callable backs the joint multi-block driver at all
    # (gcol33/tulpa#775).
    cabi = NULL, hyperprior = TRUE,
    # Has an R-level fitter (reachable, callable directly) but is NOT one
    # `tulpa(mode =)` can dispatch to: it needs multiple response arms, which
    # a single tulpa() formula cannot express (gcol33/tulpa#786). Read by the
    # `input == "nested"` dispatch in tulpa() and by inference_mode_info()'s
    # printed backend list, so the two cannot drift apart the way the
    # hardcoded `backend != "nested_laplace"` check and the `[R]` tag did.
    dispatchable = FALSE,
    note = paste("Joint multi-arm nested Laplace; driven by model packages, not the",
                 "single-response tulpa() formula (cannot express multiple arms)")
  ),
  spde = list(
    emits = "iid",
    tier = "structured", input = "spde", fitter = "fit_spde",
    families = c("binomial", "poisson", "neg_binomial_2", "gaussian"),
    # cpp_nested_laplace_spde is the Rcpp export fit_spde() calls; the
    # registered C-ABI symbol a model package reaches is
    # "tulpa_nested_laplace_spde" (gcol33/tulpa#775).
    cabi = "tulpa_nested_laplace_spde", hyperprior = TRUE,
    note = paste("Continuous Matern SPDE field; nested-Laplace integration over",
                 "(range, sigma) via fit_spde(). Uses its own CCD / grid",
                 "hyperparameter engine (the FEM Q-builder rebuilds the precision",
                 "per node), not the generic registry grid that tulpa_nested_laplace",
                 "drives -- so an SPDE field selecting nested_laplace is redirected",
                 "here in tulpa().")
  ),
  # ---- Tier 3: Optimized ----
  vi = list(
    emits = "iid",
    tier = "optimized", input = "modeldata", fitter = "tulpa_sample_glmm",
    families = NULL, cabi = "tulpa_fit_vi",
    note = "Generic VI over a fixed-effect GLM (mean-field / low-rank / full-rank)"
  )
)


#' Backend names belonging to a tier (derived from the registry)
#' @noRd
.tier_backends <- function(tier_key) {
  names(Filter(function(e) e$tier == tier_key, BACKEND_REGISTRY))
}

# Backend names whose fitter takes the outer `hyperprior` choice.
.hyperprior_backends <- function() {
  names(Filter(function(e) isTRUE(e$hyperprior), BACKEND_REGISTRY))
}


#' Inference tiers with their properties
#'
#' Derived from `TIER_META` and `BACKEND_REGISTRY`.
#' Preserves the `INFERENCE_TIERS$<tier>$backends` interface relied on by
#' downstream code and tests.
#' @noRd
INFERENCE_TIERS <- local({
  out <- vector("list", length(TIER_META))
  names(out) <- names(TIER_META)
  for (tk in names(TIER_META)) {
    out[[tk]] <- c(TIER_META[[tk]], list(backends = .tier_backends(tk)))
  }
  out
})


#' Backend -> supported families, derived from the registry
#' @noRd
BACKEND_FAMILY_SUPPORT <- local({
  out <- list()
  for (b in names(BACKEND_REGISTRY)) {
    fam <- BACKEND_REGISTRY[[b]]$families
    if (!is.null(fam)) out[[b]] <- fam
  }
  out
})


#' All registered backends across every tier (derived from the registry)
#' @noRd
ALL_BACKENDS <- names(BACKEND_REGISTRY)


# Areal spatial field types the nested-Laplace front door supports. The
# prior block is buildable from type + adjacency + a resolved spatial_idx
# alone, and a spatial(col) term names the per-obs unit. Intrinsic CAR ("car")
# shares the ICAR precision.
.NL_FRONTDOOR_AREAL <- c("icar", "car", "bym2", "car_proper")

# Continuous (coordinate-addressed) spatial field types the nested-Laplace
# front door supports via the generic registry grid (tulpa_nested_laplace).
# Addressed by coordinate columns in a spatial_gp() / spatial_gp(approx = "hsgp", ) spec
# (no spatial(col) term); tulpa() validates via validate_gp() to derive the
# unique locations, the obs -> location map, and the NNGP neighbour structure.
# spatial_gp() emits type "gp" (the nested kernel block is "nngp");
# spatial_gp(approx = "hsgp", ) emits type "hsgp" (the nested kernel block consumes a
# Laplacian basis built by cpp_hsgp_basis_2d).
.NL_FRONTDOOR_CONTINUOUS <- c("gp", "nngp", "hsgp")

# Field shapes a restricted-spatial-regression projection can be applied to:
# the ones a Polya-Gamma Gibbs kernel carries the field's own prior precision
# for. Areal goes to `cpp_pg_binomial_gibbs_rsr()` (an adjacency), the NNGP
# spellings to `cpp_pg_binomial_gibbs_gp_rsr()` (Vecchia factors at the current
# range). An HSGP basis and an SPDE mesh have neither and are refused at
# `spatial_rsr()` (gcol33/tulpa#815, gcol33/tulpa#848).
.RSR_FIELDS <- c(.NL_FRONTDOOR_AREAL, "gp", "nngp")

# Of those, the ones whose projection is applied on the CONTINUOUS kernel.
.RSR_CONTINUOUS <- c("gp", "nngp")

# The `spatial_type` tulpa() hands the backend selector for a restricted field:
# "rsr" for an areal one and "gp_rsr" for a continuous one. Both also ride the
# nested-Laplace multi-block driver, which reads the unit projector as the
# block's `projector`.
.rsr_spatial_type <- function(field_type) {
  if (field_type %in% .RSR_CONTINUOUS) "gp_rsr" else "rsr"
}

# SPDE is also coordinate-addressed (no spatial(col) term) and nested-integrated,
# but it carries its OWN nested engine -- fit_spde() rebuilds the Matern precision
# Q(range, sigma) per node via the FEM Q-builder and integrates (range, sigma)
# with a CCD / grid design in R, not the generic registry grid above. The spec
# (spatial_spde() / spatial_spde_custom()) is self-contained (mesh + FEM matrices
# built at construction), so there is no validate_* step. Kept separate from
# .NL_FRONTDOOR_CONTINUOUS so the grammar branch knows it is coord-addressed yet
# does NOT route it through validate_gp() / tulpa_nested_laplace(); the dedicated
# `spde` backend (redirected from nested_laplace in tulpa()) drives fit_spde().
.NL_FRONTDOOR_SPDE <- "spde"

# The multi-scale (local + regional) NNGP field is coordinate-addressed like
# .NL_FRONTDOOR_CONTINUOUS and derives its two neighbour structures through the
# same validate_gp(), but it has no nested-Laplace kernel: the exact-NUTS
# sampler (compute_multiscale_gp_prior + the msgp_nc_apply transform) is its
# only front-door integrator. Kept out of .NL_FRONTDOOR_NESTED so auto /
# structured leave it on the conditional Laplace path instead of routing it to
# a nested kernel that does not exist.
.FRONTDOOR_MULTISCALE <- "multiscale"

# Spatial field types tulpa() routes to a nested-Laplace integrator (areal +
# continuous + SPDE). Single source of truth for the auto/structured router
# (auto_select_mode, select_backend_for_mode). The areal + .NL_FRONTDOOR_CONTINUOUS
# subset goes through the generic nested_laplace backend (and the spec -> prior
# converter .spatial_spec_to_nl_prior); SPDE is redirected to the `spde` backend.
# Types outside this set stay on the conditional Laplace path under auto/structured.
.NL_FRONTDOOR_NESTED <- c(.NL_FRONTDOOR_AREAL, .NL_FRONTDOOR_CONTINUOUS,
                          .NL_FRONTDOOR_SPDE, "rsr", "gp_rsr")

# Spatial types that can share the nested-Laplace multi-block stack with a
# random-effect term, a temporal field, smoothers or latent blocks: the areal
# and continuous fields, restricted or not.
.NL_FRONTDOOR_STACKABLE <- c(.NL_FRONTDOOR_AREAL, .NL_FRONTDOOR_CONTINUOUS,
                             "rsr", "gp_rsr")

# The subset of .NL_FRONTDOOR_NESTED the exact ModelData NUTS sampler (`hmc`
# and its siblings) also threads directly -- read off dispatch_glmm_modeldata()
# (R/tulpa.R): icar / bym2 / car_proper via the areal spec, gp / nngp / hsgp via
# .gp_sampler_spec() / .hsgp_sampler_spec(). Plain intrinsic 'car' and 'spde'
# are NOT in this set -- the modeldata dispatch refuses 'car' by name (areal
# there is icar/bym2 only) and has no SPDE spec builder at all -- so auto must
# not offer either as a fallback when nested_laplace cannot carry a call
# (gcol33/tulpa#769).
.NL_FRONTDOOR_HMC_SPATIAL <- c("icar", "bym2", "car_proper",
                               .NL_FRONTDOOR_CONTINUOUS)


#' Test whether a backend supports the given family
#' @noRd
backend_supports_family <- function(backend, family) {
  supported <- BACKEND_REGISTRY[[backend]]$families
  if (is.null(supported)) return(TRUE)  # unrestricted backend

  fam_name <- family$name %||% ""
  fam_dist <- family$distribution %||% (family$numerator$distribution %||% "")
  fam_name %in% supported || fam_dist %in% supported
}


#' Redirect a backend selection: swap the backend and restamp mode / tier /
#' tier name from the registry, recording the reason shown on the fit
#'
#' A redirect off an EXPLICIT request (anything but `mode = "auto"`) is an
#' override, not a resolution: the fit does what the model structure requires
#' instead of what the caller asked for. That is recorded on the selection --
#' `sel$overridden` for callers, and a clause appended to the reason the fit
#' reports -- so an overridden fit is distinguishable from one that was never
#' asked for a mode at all. Only the FIRST override is recorded: with several
#' redirects in a chain, the request the user actually made is the one worth
#' naming, not the intermediate backend a previous redirect chose.
#'
#' `notify` says whether tulpa() should additionally WARN. It is for the case
#' where the requested backend could have fitted the model and the redirect takes
#' a capability away -- a smoother sending an `eb` / `agq` request to the nested
#' kernels loses the random-effect SD those two would have estimated. It is off
#' where the requested mode is not expressible for the structure at all and the
#' redirect is the documented resolution rather than a loss: a random slope has no
#' scalar `sigma_re` for `mode = "laplace"` to condition on, and an SPDE field
#' redirected from `nested_laplace` to `spde` is the same mode and tier reaching
#' its own integrator. Those are recorded, not warned about, so a documented route
#' does not warn on every fit.
#' @noRd
.sel_redirect <- function(sel, backend, reason, notify = TRUE) {
  if (isTRUE(sel$explicit) && is.null(sel$overridden) &&
      !identical(sel$backend, backend)) {
    sel$overridden <- list(requested = sel$requested %||% sel$mode,
                           backend   = sel$backend,
                           notify    = isTRUE(notify))
  }
  # Appended to whatever reason this selection ends up carrying, not only to the
  # reason of the redirect that recorded the override: each redirect REPLACES
  # `sel$reason`, so in a chain the later one would otherwise drop the statement
  # and the fit would again look like it was never asked for a mode.
  if (!is.null(sel$overridden)) {
    # A tier mode ("structured") names a different thing than the backend it
    # resolved to ("laplace"), and both are worth stating; a request that IS a
    # backend name would otherwise print it twice.
    req <- if (identical(sel$overridden$requested, sel$overridden$backend)) {
      sprintf("mode = '%s'", sel$overridden$requested)
    } else {
      sprintf("mode = '%s' (backend '%s')", sel$overridden$requested,
              sel$overridden$backend)
    }
    reason <- sprintf("%s -- overrides the requested %s", reason, req)
  }
  sel$backend <- backend
  ti <- get_backend_tier(backend)
  sel$mode      <- ti$mode
  sel$tier      <- ti$tier
  sel$tier_name <- ti$name
  sel$reason    <- reason
  sel
}

#' Get tier for a backend
#' @param backend Character string naming the backend
#' @return List with tier information
#' @noRd
get_backend_tier <- function(backend) {
  entry <- BACKEND_REGISTRY[[backend]]
  if (is.null(entry)) {
    stop(sprintf("Unknown backend: '%s'", backend), call. = FALSE)
  }
  meta <- TIER_META[[entry$tier]]
  list(
    tier = meta$tier,
    name = meta$name,
    mode = entry$tier,
    description = meta$description,
    guarantee = meta$guarantee,
    note = meta$note
  )
}


# ==============================================================================
# Reachability + dispatch spine
# ==============================================================================

#' Is a backend reachable from R (does it have an R-level fitter)?
#' @noRd
backend_is_reachable <- function(backend) {
  entry <- BACKEND_REGISTRY[[backend]]
  !is.null(entry) && !is.null(entry$fitter)
}


#' Error if a selected backend has no R-level fitter
#'
#' @description
#' Enforces the registry honesty contract: a backend may ship a C++ kernel
#' reachable from model packages (via `LinkingTo: tulpa`) yet have no R entry
#' point. Such a backend must never be *silently* selectable from R --
#' selecting it errors with a precise message naming the C-ABI symbol, rather
#' than pretending to dispatch.
#'
#' @noRd
assert_backend_reachable <- function(backend) {
  entry <- BACKEND_REGISTRY[[backend]]
  if (is.null(entry)) {
    stop(sprintf("Unknown backend: '%s'", backend), call. = FALSE)
  }
  if (is.null(entry$fitter)) {
    meta <- TIER_META[[entry$tier]]
    cabi <- paste(entry$cabi %||% "(none)", collapse = ", ")
    stop(sprintf(paste0(
      "Backend '%s' is registered (Tier %d, %s) but has no R-level fitter.\n",
      "Its C++ kernel is reachable from model packages through the C ABI\n",
      "(%s), but it cannot be driven directly from R yet. Choose a reachable\n",
      "backend, or wire an R wrapper over the kernel."),
      backend, meta$tier, meta$name, cabi),
      call. = FALSE)
  }
  invisible(TRUE)
}


#' Resolve the R fitter function for a backend (errors if unreachable)
#'
#' Looks up the fitter *name* string in the registry and resolves it lazily,
#' so the registry stays independent of source-file load order.
#' @noRd
resolve_backend_fitter <- function(backend) {
  assert_backend_reachable(backend)
  match.fun(BACKEND_REGISTRY[[backend]]$fitter)
}


#' Dispatch a fit to the backend chosen by the mode/tier system
#'
#' @description
#' The routing spine. Resolves `mode` to a concrete backend via
#' [select_inference_mode()], asserts the backend is R-reachable (fails loudly
#' otherwise), then calls its fitter with `fitter_args`. The caller supplies
#' `fitter_args` matching the backend's input contract
#' (`BACKEND_REGISTRY$<backend>$input`).
#'
#' The selected mode/tier/backend are stamped onto the returned fit (without
#' overwriting any the fitter already set), so the inference contract is always
#' visible in the output.
#'
#' @param mode User-specified mode (`"auto"`, a tier, or a backend name).
#' @param fitter_args Named list of arguments forwarded to the backend fitter.
#' @param family,n_obs,has_spatial,has_temporal,has_latent,spatial_type,temporal
#'   Model characteristics forwarded to [select_inference_mode()].
#' @return The fitter's result, with `inference_mode`, `inference_tier`,
#'   `backend`, and `selection_reason` ensured.
#' @noRd
tulpa_dispatch <- function(mode,
                           fitter_args = list(),
                           family = NULL,
                           n_obs = NULL,
                           has_spatial = FALSE,
                           has_temporal = FALSE,
                           has_latent = FALSE,
                           spatial_type = NULL,
                           temporal = NULL) {
  sel <- select_inference_mode(
    mode,
    family = family, n_obs = n_obs,
    has_spatial = has_spatial, has_temporal = has_temporal,
    has_latent = has_latent, spatial_type = spatial_type, temporal = temporal
  )

  fitter <- resolve_backend_fitter(sel$backend)
  fit <- do.call(fitter, fitter_args)

  if (is.list(fit)) {
    fit$inference_mode <- fit$inference_mode %||% sel$mode
    fit$inference_tier <- fit$inference_tier %||% sel$tier
    fit$selection_reason <- fit$selection_reason %||% sel$reason
    fit <- .finalize_fit(fit, backend = sel$backend)
  }
  fit
}


#' Attach the standard `tulpa_fit` contract to a fitter's result
#'
#' The dispatch layer and the fitters that return generic-accessor-facing
#' objects route their return value through this helper (some special-purpose
#' fitters -- tulpa_ep, the categorical drivers -- still
#' stamp their class by hand), so a directly-called fitter and a
#' `tulpa()`-dispatched one yield the same
#' enriched object: the `tulpa_fit` class (so the generic S3 methods --
#' `coef` / `summary` / `vcov` / `confint` / `tidy` / `glance` / `ranef` --
#' dispatch), the fixed-effect layout (`n_fixed` / `fixed_names` / `param_names`,
#' which the `.fit_fixed_table` summary path reads), and an explicit
#' posterior-provenance tag (`draws_kind`) the chain-vs-iid diagnostic gate
#' (`.tulpa_draws_kind()`) reads to decide whether Rhat/ESS apply. Each field is
#' filled only when the fitter did not already set it, so a fitter that knows
#' better wins and the helper is idempotent under `tulpa_dispatch()`.
#'
#' `draws_kind` precedence is: a value already on the fit, then the explicit
#' `draws_kind` argument, then the registry `emits` property for `backend`.
#' Backends that are not registry keys (the `tgmrf_*` fitters) must pass
#' `draws_kind` explicitly, since their `emits` cannot be looked up.
#'
#' @param fit The fitter result (a list); returned unchanged if not a list.
#' @param backend Backend key (sets `$backend`; supplies the default
#'   `draws_kind` via the registry `emits` property when it is a registry key).
#' @param draws_kind Explicit `"chain"` / `"iid"` / `"point"` tag; used when the
#'   backend is absent from `BACKEND_REGISTRY` or to override the registry.
#' @param n_fixed,fixed_names,param_names Fixed-effect layout, each filled only
#'   when the fitter left it unset. An unnamed `means` vector of the matching
#'   length is named by `param_names`.
#' @param data Observation-level fields the fitter was called with --
#'   `list(y=, n_trials=, model_matrix=, family=, offset=, phi=, phi2=)` --
#'   stamped onto the fit so a directly-called fitter carries the same
#'   `nobs()` / `fitted()` / `predict()` / `posterior_predict()` surface a
#'   `tulpa()`-dispatched one does (gcol33/tulpa#781). Any subset may be
#'   supplied; each entry fills only when the fitter did not already set the
#'   matching field, so a fitter that already stamped a richer value (e.g. the
#'   full design bundle `tulpa()` attaches post-dispatch) wins.
#' @param extra_class Subclass(es) to prepend before `tulpa_fit`.
#' @return The enriched fit, classed `c(extra_class, ..., "tulpa_fit")`.
#' @noRd
.finalize_fit <- function(fit, backend = NULL, draws_kind = NULL,
                          n_fixed = NULL, fixed_names = NULL,
                          param_names = NULL, extra_class = NULL,
                          data = NULL) {
  if (!is.list(fit)) return(fit)

  if (!is.null(backend)) fit$backend <- fit$backend %||% backend
  reg <- if (!is.null(fit$backend)) BACKEND_REGISTRY[[fit$backend]] else NULL
  fit$draws_kind <- fit$draws_kind %||% draws_kind %||% reg$emits

  if (!is.null(data)) {
    fit$y            <- fit$y %||% data$y
    fit$n_trials     <- fit$n_trials %||% data$n_trials
    fit$model_matrix <- fit$model_matrix %||% data$model_matrix
    fit$family       <- fit$family %||% data$family
    fit$offset       <- fit$offset %||% data$offset
    fit$phi          <- fit$phi %||% data$phi
    fit$phi2         <- fit$phi2 %||% data$phi2
    # Was that `phi` estimated or conditioned on? Only `tulpa_eb()` set the
    # field, so on every other door a consumer could not tell a fitted
    # dispersion from the default (gcol33/tulpa#849). FALSE is the answer
    # everywhere else, and it is set only for the families that read a
    # dispersion at all -- elsewhere there is no quantity to have estimated.
    if (is.null(fit$phi_estimated) &&
        .family_base(fit$family) %in% .PHI_FAMILIES) {
      fit$phi_estimated <- FALSE
    }
  }

  # Axis fields the fit's own path could not read, published
  # for the duration of the fit by the nested-Laplace front doors.
  fit$axis_fields_dropped <- fit$axis_fields_dropped %||%
    getOption("tulpa.nl_axis_dropped", NULL)

  # A zero-length name vector (`character(0)`, from `colnames(X) %||%
  # positional` when `X` has NULL colnames) is not information -- treat it
  # the same as NULL rather than latching it in, which previously left
  # `fixed_names` at `character(0)` and turned `fixed_names[idx]` into NA for
  # every coefficient (gcol33/tulpa#780). `%||%` only substitutes for NULL, so
  # both sides of the assignment are length-checked instead: an already-set
  # `character(0)` (from a prior finalize pass, or a fitter that pre-set the
  # field) does not block a later, real name vector either.
  if (!is.null(n_fixed) && is.null(fit$n_fixed)) fit$n_fixed <- n_fixed
  if (length(fixed_names) && !length(fit$fixed_names)) fit$fixed_names <- fixed_names
  if (length(param_names) && !length(fit$param_names)) fit$param_names <- param_names

  # A blank ("") name survives every length/NULL check above but is refused
  # by R's OWN `[.matrix` when later used to select a column -- `m[, ""]`
  # errors "subscript out of bounds" even when a column is literally named
  # "" (`cbind(1, x)` names its first column ""), which is what
  # `.tulpa_point_linpred()` does with `names(coef(fit))` at predict time
  # (gcol33/tulpa#780). Repaired positionally as soon as `fixed_names` holds
  # SOME real content; a fit carrying no names at all (mala(), pathfinder(),
  # a raw `logpost` sampler with no fixed/random split) is left alone here --
  # synthesizing a placeholder for those would block the later, real name a
  # front door adds through a SECOND `.finalize_fit()` call (mode = "mala"
  # wrapped by tulpa()); this only ever narrows an already-nonempty vector.
  if (!is.null(fit$n_fixed) && length(fit$fixed_names)) {
    fit$fixed_names <- .tulpa_fill_names(fit$fixed_names, fit$n_fixed, "beta")
  }
  # The design a standalone fitter stores verbatim (`$model_matrix <- X`) is
  # what predict() / fitted() / residuals() later index BY NAME against
  # `names(coef(fit))`; if `X` itself had no or blank column names, stamp the
  # SAME repaired names onto its leading `fixed_names` columns so the two
  # stay in lock-step rather than one saying "beta1" and the other "".
  # `tulpa()`'s own dispatch overwrites `$model_matrix` unconditionally with
  # the (always fully-named) design bundle afterward, so this never fights a
  # later real design.
  if (length(fit$fixed_names) && is.matrix(fit$model_matrix) &&
      ncol(fit$model_matrix) >= length(fit$fixed_names)) {
    mmcn <- colnames(fit$model_matrix)
    p <- length(fit$fixed_names)
    if (is.null(mmcn)) mmcn <- rep("", ncol(fit$model_matrix))
    if (!identical(mmcn[seq_len(p)], fit$fixed_names)) {
      mmcn[seq_len(p)] <- fit$fixed_names
      colnames(fit$model_matrix) <- mmcn
    }
  }
  fit <- .name_means_by_parameter(fit)

  cls <- oldClass(fit)
  if (!is.null(extra_class)) cls <- c(setdiff(extra_class, cls), cls)
  if (!("tulpa_fit" %in% cls)) cls <- c(cls, "tulpa_fit")
  class(fit) <- cls
  fit
}


# Fill NULL, blank ("") or short/absent names positionally with `<prefix><i>`.
# The one place a name vector reaching a boundary that cannot tolerate NA, "",
# or a length mismatch is repaired -- `row.names<-` (print/summary), a
# `draws[, ""]` subset (the per-parameter plots), and `dimnames<-` on an
# unnamed array all error opaquely on such a vector rather than reporting "no
# name" (gcol33/tulpa#780). Deliberately NOT written back onto the fit: a
# fitter's raw result is often re-finalized by a front door that supplies the
# real names afterward (`tulpa()` re-finalizing a backend's own fit with the
# design's own `param_names`), and persisting a synthesized "param3" would
# block that later, authoritative assignment the same way a stale
# `character(0)` used to. Every reader that needs a complete name vector
# (`.tulpa_draw_names()`, `.tulpa_pooled_draws()`) calls this at read time
# instead, off whatever the fit carries at that moment.
#' @keywords internal
.tulpa_fill_names <- function(nm, n, prefix = "param") {
  out <- rep(NA_character_, n)
  if (!is.null(nm) && length(nm)) {
    k <- min(length(nm), n)
    out[seq_len(k)] <- as.character(nm)[seq_len(k)]
  }
  blank <- is.na(out) | !nzchar(out)
  if (any(blank)) out[blank] <- paste0(prefix, which(blank))
  out
}


# Name the fit's `means` by `param_names`, the full-parameter vector both are
# laid out in. A producer that samples a bare numeric vector (the log-posterior
# samplers) cannot name it itself, so the names arrive with the layout. A vector
# already named, or of another length, is left as it is.
#' @keywords internal
.name_means_by_parameter <- function(fit) {
  pn <- fit$param_names
  if (is.character(pn) && length(pn) && is.numeric(fit$means) &&
      is.null(names(fit$means)) && length(fit$means) == length(pn)) {
    names(fit$means) <- pn
  }
  fit
}


# Rename the first `p` columns of a NUTS producer's `draws` matrix (and
# `means`, if named in step) from the C++ kernel's positional `beta[j]` to the
# design matrix's own column names -- the same names the Laplace path reports
# coefficients under. Declines (leaves `beta[j]`) rather than renaming under a
# duplicated or blank name.
#' @keywords internal
.nuts_name_fixed_effects <- function(res, X_names, p) {
  if (is.null(X_names) || length(X_names) != p ||
      anyDuplicated(X_names) || !all(nzchar(X_names))) {
    return(res)
  }
  dn <- colnames(res$draws)
  dn[seq_len(p)] <- X_names
  colnames(res$draws) <- dn
  if (!is.null(names(res$means)) && length(res$means) == length(dn)) {
    names(res$means)[seq_len(p)] <- X_names
  }
  res
}


#' Map mode to valid backends
#' @param mode Character: "auto", "exact", "structured", or "optimized"
#' @return Character vector of valid backends for this mode
#' @noRd
get_mode_backends <- function(mode) {
  mode <- tolower(mode)

  if (mode == "auto") {
    # Auto can choose Tier 1 or Tier 2, never Tier 3
    return(c(INFERENCE_TIERS$exact$backends, INFERENCE_TIERS$structured$backends))
  }

  if (mode %in% names(INFERENCE_TIERS)) {
    return(INFERENCE_TIERS[[mode]]$backends)
  }

  stop(sprintf(
    "Unknown mode: '%s'. Use one of: %s",
    mode,
    paste(c("auto", names(INFERENCE_TIERS)), collapse = ", ")
  ), call. = FALSE)
}


# ==============================================================================
# Mode Selection Logic
# ==============================================================================

#' Select inference mode and backend
#'
#' @description
#' Implements the mode selection logic for tulpa. Accepts either tier names
#' (auto, exact, structured, optimized) or backend names (hmc, ess, gibbs, laplace, vi).
#'
#' When mode is "auto", selects between Tier 1 (Exact) and Tier 2 (Structured)
#' based on model characteristics. Never selects Tier 3 (Optimized) automatically.
#'
#' @param mode User-specified mode or backend name
#' @param family Model family object
#' @param n_obs Number of observations
#' @param has_spatial Whether model has spatial effects
#' @param has_temporal Whether model has temporal effects
#' @param has_latent Whether model has latent factors
#' @param has_re Whether the model has any random-effect term (`(1 | g)` or
#'   `(1 + x | g)`), scalar or slope alike
#'
#' @return List with:
#'   - mode: The selected mode name
#'   - backend: The selected backend
#'   - tier: The tier number (1, 2, or 3)
#'   - tier_name: The tier name
#'   - reason: Explanation for the selection
#'
#' @noRd
select_inference_mode <- function(mode,
                                  family,
                                  n_obs,
                                  has_spatial = FALSE,
                                  has_temporal = FALSE,
                                  has_latent = FALSE,
                                  spatial_type = NULL,
                                  temporal = NULL,
                                  has_re = FALSE,
                                  feat = list()) {

  mode <- tolower(mode)

  # Check if mode is actually a backend name (derived from INFERENCE_TIERS)
  if (mode %in% ALL_BACKENDS) {
    # User specified a backend directly
    backend <- mode
    tier_info <- get_backend_tier(backend)

    return(list(
      mode = tier_info$mode,
      backend = backend,
      tier = tier_info$tier,
      tier_name = tier_info$name,
      reason = sprintf("User-specified backend: %s", backend),
      # The literal request, kept so a later redirect can name what it overrode.
      requested = backend,
      explicit = TRUE
    ))
  }

  # Validate tier name
  valid_modes <- c("auto", names(INFERENCE_TIERS))
  if (!mode %in% valid_modes) {
    stop(sprintf(
      "Invalid mode: '%s'. Use one of: %s, or a backend: %s",
      mode,
      paste(valid_modes, collapse = ", "),
      paste(ALL_BACKENDS, collapse = ", ")
    ), call. = FALSE)
  }

  # Auto-selection logic. Marked non-explicit: a redirect off an auto selection
  # is the router doing its job, not an override of anything the caller asked for.
  if (mode == "auto") {
    sel <- auto_select_mode(family, n_obs, has_spatial, has_temporal, has_latent, temporal,
                            spatial_type, has_re = has_re, feat = feat)
    sel$requested <- "auto"
    sel$explicit  <- FALSE
    return(sel)
  }

  # Explicit tier mode: select best backend within that tier
  backend <- select_backend_for_mode(mode, family, n_obs, has_spatial, has_temporal,
                                     has_latent, spatial_type, has_re = has_re, feat = feat)
  tier_info <- get_backend_tier(backend)

  return(list(
    mode = mode,
    backend = backend,
    tier = tier_info$tier,
    tier_name = tier_info$name,
    reason = sprintf("User-specified mode: %s", mode),
    requested = mode,
    explicit = TRUE
  ))
}


#' Auto-select mode (Tier 1 or Tier 2 only)
#'
#' @description
#' Implements the "auto" mode selection. Chooses the most reliable method
#' that is expected to finish for the given model.
#'
#' **Critical rule**: Auto never selects Tier 3 (Optimized).
#'
#' @noRd
auto_select_mode <- function(family, n_obs, has_spatial, has_temporal, has_latent, temporal = NULL,
                             spatial_type = NULL, has_re = FALSE,
                             feat = list()) {
  # `feat` carries the per-call features a backend can refuse at dispatch; see
  # .auto_backend_ok(). The `spatial` flag is set from has_spatial so a caller
  # cannot pass the two inconsistently.
  feat$spatial <- isTRUE(has_spatial)

  # Latent prior blocks (`latent(tgmrf(...))`) integrate their hyperparameters
  # via nested Laplace -- the designed Tier 2 hot path for latent Gaussian
  # structure. This takes precedence over the dataset-size heuristics below:
  # a latent-block model goes nested regardless of n.
  if (has_latent) {
    return(list(
      mode = "structured",
      backend = "nested_laplace",
      tier = 2L,
      tier_name = "Structured",
      reason = "latent prior block(s); nested-Laplace hyperparameter integration"
    ))
  }

  # Spatially- / temporally-varying coefficients are carried only by the exact
  # ModelData NUTS backend (`hmc`); there is no nested-Laplace SVC/TVC front
  # door. Any other backend errors at the tulpa() varying-coefficient guard, so
  # auto must pick the modeldata sampler directly rather than fall through to a
  # Laplace/size heuristic that dispatch would reject.
  if (identical(tolower(spatial_type %||% ""), "svc") ||
      (!is.null(temporal) && identical(tolower(temporal$type %||% ""), "tvc"))) {
    vc <- if (identical(tolower(spatial_type %||% ""), "svc"))
      "spatially" else "temporally"
    return(list(
      mode = "exact", backend = "hmc", tier = 1L, tier_name = "Exact",
      reason = sprintf("%s-varying coefficients (exact ModelData NUTS)", vc)
    ))
  }

  # A continuous-time temporal GP and the multi-scale temporal field are
  # sampler-path only: their hyperparameters are sampled jointly with the field
  # and there is no nested-Laplace kernel laying a grid over a dense T x T
  # Gaussian. auto used to fall through to the temporal branch below, which has
  # arms for rw1 / rw2 / ar1 only, and tulpa() then redirected to nested Laplace
  # -- the one door that refuses them (gcol33/tulpa#672).
  t_type <- tolower((temporal$type %||% "")[1])
  if (t_type %in% c("gp", "multiscale", "temporal_multiscale")) {
    return(list(
      mode = "exact", backend = "hmc", tier = 1L, tier_name = "Exact",
      reason = sprintf("continuous-time temporal field (%s); exact ModelData NUTS",
                       t_type)
    ))
  }

  # A plain temporal field (rw1 / rw2 / ar1; gp / multiscale are the sampler-
  # path-only arm just above) integrates through the nested-Laplace temporal
  # kernel, same as an areal spatial field below. With no arm here, auto fell
  # through unconditionally to the gradient-sampler default (picked on family
  # and dataset size alone) and relied on a redirect further downstream in
  # tulpa() to move a temporal-only model to nested_laplace -- but a ziformula
  # / weights / phi2 the default backend could not carry then errored at the
  # front door's feature guards BEFORE that redirect ever ran (gcol33/tulpa#769).
  # Deciding it here means the call's actual features are visible before a
  # backend is picked at all, and the fallback mirrors the spatial arm below:
  # the exact ModelData sampler consumes a temporal field directly and carries
  # ziformula / phi2 where the nested kernel does not.
  if (has_temporal && !has_spatial) {
    if (.auto_backend_ok("nested_laplace", family, feat)) {
      return(list(
        mode = "structured", backend = "nested_laplace", tier = 2L,
        tier_name = "Structured",
        reason = sprintf("temporal %s field; nested-Laplace integration", t_type)
      ))
    }
    if (.auto_backend_ok("hmc", family, feat)) {
      return(list(
        mode = "exact", backend = "hmc", tier = 1L, tier_name = "Exact",
        reason = sprintf(paste0(
          "temporal %s field; nested-Laplace cannot carry this call, exact ",
          "ModelData NUTS instead"), t_type)
      ))
    }
  }

  # The multi-scale (local + regional) field is carried only by the exact
  # ModelData NUTS backend too: it has no nested-Laplace kernel, and
  # dispatch_laplace_spatial() rejects it, so the same rule applies -- pick the
  # sampler that fits it rather than fall through to a Laplace heuristic whose
  # dispatch errors.
  if (identical(tolower(spatial_type %||% ""), .FRONTDOOR_MULTISCALE)) {
    return(list(
      mode = "exact", backend = "hmc", tier = 1L, tier_name = "Exact",
      reason = "multi-scale spatial field (exact ModelData NUTS)"
    ))
  }

  # Spatial latent Gaussian field. Like a latent block, its hyperparameter
  # (spatial precision / mixing) is integrated -- not conditioned at a fixed
  # scale, which is the explicit mode = "laplace". Three treatments, in
  # preference order:
  #  * binomial areal (icar/bym2/rsr): exact component-wise Polya-Gamma Gibbs
  #    (Tier 1). The PG samplers update the field component-wise, avoiding
  #    HMC's curse of dimensionality. A restricted field (spatial_rsr()) goes
  #    here too, areal or continuous; a continuous one has no other backend.
  #    An areal one under another family reaches nested Laplace below.
  #    dispatch_gibbs_spatial() also handles gp / nngp / multiscale_gp, but auto
  #    deliberately routes those to the more general nested path below (their
  #    Gibbs samplers need one observation per location); they stay reachable via
  #    the explicit back door tulpa_gibbs(spatial=). auto must never pick a
  #    backend that errors at dispatch.
  #  * areal (icar/car/bym2/car_proper), continuous gp/nngp/hsgp, and SPDE:
  #    nested Laplace (Tier 2) integrates the spatial hyperparameter through the
  #    tulpa() front door (.NL_FRONTDOOR_NESTED). The areal + gp/nngp/hsgp subset
  #    runs the generic nested_laplace backend; SPDE is redirected in tulpa() to
  #    the dedicated `spde` backend (its own fit_spde CCD / grid engine).
  #  * other continuous field types: the conditional Laplace path.
  if (has_spatial && !is.null(spatial_type)) {
    fam_nm <- family$name %||% family$distribution %||% ""
    # A temporal field routes through the joint nested-Laplace path in tulpa()
    # (the design-input Gibbs backend is redirected there), so Gibbs is reachable
    # from auto only when there is no temporal field. spatial_car()'s exported
    # "car" is the same intrinsic ICAR field dispatch_gibbs_spatial() dispatches
    # under "icar" (gcol33/tulpa#819).
    if (.areal_gibbs_type(spatial_type) %in% c("icar", "bym2", "rsr", "gp_rsr") &&
        identical(fam_nm, "binomial") &&
        !has_temporal && .auto_backend_ok("gibbs", family, feat)) {
      return(list(
        mode = "exact", backend = "gibbs", tier = 1L, tier_name = "Exact",
        reason = sprintf("%s spatial model (binomial Polya-Gamma Gibbs)", spatial_type)
      ))
    }
    if (spatial_type %in% .NL_FRONTDOOR_NESTED) {
      # The generic nested-Laplace converter carries no ziformula / (rarely)
      # other feature the call needs (gcol33/tulpa#769: an icar field +
      # ziformula picked this backend and then errored at the front door's
      # zero-inflation guard). The exact ModelData NUTS sampler threads areal
      # (icar/bym2/car_proper) and continuous (gp/nngp/hsgp) fields directly
      # and DOES carry ziformula/phi2 -- but not SPDE or plain intrinsic 'car'
      # (.NL_FRONTDOOR_HMC_SPATIAL), so only that subset gets the fallback;
      # elsewhere this keeps returning nested_laplace exactly as before and the
      # downstream guard produces its usual (unchanged) refusal.
      if (!.auto_backend_ok("nested_laplace", family, feat) &&
          spatial_type %in% .NL_FRONTDOOR_HMC_SPATIAL &&
          .auto_backend_ok("hmc", family, feat)) {
        return(list(
          mode = "exact", backend = "hmc", tier = 1L, tier_name = "Exact",
          reason = sprintf(paste0(
            "%s spatial field; nested-Laplace cannot carry this call, exact ",
            "ModelData NUTS instead"), spatial_type)
        ))
      }
      return(list(
        mode = "structured", backend = "nested_laplace", tier = 2L,
        tier_name = "Structured",
        reason = sprintf("%s spatial field; nested-Laplace hyperparameter integration",
                         spatial_type)
      ))
    }
    return(list(
      mode = "structured", backend = "laplace", tier = 2L, tier_name = "Structured",
      reason = sprintf("continuous spatial field (%s); conditional Laplace", spatial_type)
    ))
  }

  # Non-spatial size heuristic: very large designs go to Laplace (Tier 2).
  VERY_LARGE <- 50000
  if (n_obs > VERY_LARGE) {
    return(list(
      mode = "structured", backend = "laplace", tier = 2L, tier_name = "Structured",
      reason = sprintf("large dataset (n=%s)", format(n_obs, big.mark = ","))
    ))
  }

  # A random-effect term ((1 | g) or (1 + x | g)) has no dataset-implied scalar
  # sigma_re for auto to condition on -- conditioning at an unsupplied value
  # (silently fixed at 1) is not "the most reliable method that is expected to
  # finish for this model" when the covariance can instead be inferred. Route
  # to the exact Metropolis-within-Gibbs debias of the RE covariance
  # (re_cov_gibbs, still Tier 1) rather than falling through to the plain
  # fixed-effect default below. re_cov_gibbs already treats a scalar `(1 | g)`
  # term as the degenerate c = 1 covariance block, so a random-intercept-only
  # model is routed here exactly like a random-slope model -- no
  # special-casing by term shape. An explicit
  # mode = "laplace" / "mala" / ... still conditions on sigma_re (defaulting
  # to 1, with a warning) when the caller names it directly; only the auto
  # default changes.
  if (has_re) {
    # The exact debias first, the deterministic integration second: both infer
    # the covariance, and the second is what carries an offset, a ziformula, a
    # second dispersion and the families the Gibbs sweep does not.
    if (.auto_backend_ok("re_cov_gibbs", family, feat)) {
      return(list(
        mode = "exact", backend = "re_cov_gibbs", tier = 1L, tier_name = "Exact",
        reason = paste("random-effect term(s); RE covariance(s) integrated via",
                       "exact Metropolis-within-Gibbs debias")
      ))
    }
    if (.auto_backend_ok("re_cov_nested", family, feat)) {
      return(list(
        mode = "structured", backend = "re_cov_nested", tier = 2L,
        tier_name = "Structured",
        reason = paste("random-effect term(s); RE covariance(s) integrated via",
                       "nested Laplace (the exact debias does not carry this call)")
      ))
    }
  }

  # Default: Tier 1 (Exact). Both MALA and NUTS (backend "hmc", via
  # tulpa_sample_glmm) are reachable from the R front door; auto keeps MALA
  # as the default gradient sampler for its lower per-iteration cost on the
  # plain-GLM designs that reach this branch. mode = "exact" selects NUTS.
  #
  # MALA carries no ziformula (gcol33/tulpa#769: `tulpa(y ~ x, ..., ziformula =
  # ~ 1)` picked this branch unconditionally and errored at the front door's
  # zero-inflation guard on the very call that picked it). Gate on
  # .auto_backend_ok() and fall back to the exact ModelData sampler, which
  # carries ziformula / phi2 / offset at once, before refusing. It does NOT
  # carry weights (the ModelData likelihood has no per-row multiplier), and a
  # weighted call reaching here refuses rather than fits unweighted.
  if (.auto_backend_ok("mala", family, feat)) {
    return(list(
      mode = "exact", backend = "mala", tier = 1L, tier_name = "Exact",
      reason = "default (MALA gradient sampler; mode = 'exact' selects NUTS)"
    ))
  }
  if (.auto_backend_ok("hmc", family, feat)) {
    return(list(
      mode = "exact", backend = "hmc", tier = 1L, tier_name = "Exact",
      reason = "MALA cannot carry this call; exact ModelData NUTS instead"
    ))
  }
  stop(sprintf(paste0(
    "auto: no backend fits this call's feature combination (%s). Use an ",
    "explicit mode naming a backend that carries it, or drop the ",
    "conflicting feature."),
    paste(names(feat)[vapply(feat, isTRUE, logical(1))], collapse = ", ")),
    call. = FALSE)
}


# Can `backend` take this CALL, not just this model?
#
# `backend_supports_family()` answers the family half. The other half is the
# per-call features a backend refuses at dispatch -- an offset, weights, a
# ziformula, a second dispersion, more random-effect terms than its sweep
# updates -- and the auto selector could not see any of them, so it picked
# backends that then errored on the very call that selected them
# (gcol33/tulpa#666, #681). Every restriction is read from the registry or from
# the registry-derived capability sets, never restated here, so declaring one is
# a registry edit rather than another condition in the selector.
#
# `feat` is a named list of logical flags plus `n_re_terms`; an absent flag is
# FALSE.
#' @keywords internal
.auto_backend_ok <- function(backend, family, feat = list()) {
  if (!backend_supports_family(backend, family)) return(FALSE)
  reg <- BACKEND_REGISTRY[[backend]]
  if (is.null(reg)) return(FALSE)
  on <- function(k) isTRUE(feat[[k]])
  if (on("offset") && identical(reg$carries_offset, FALSE)) return(FALSE)
  if (on("ziformula") && !backend %in% .zi_backends()) return(FALSE)
  if (on("phi2") && !backend %in% .phi2_backends()) return(FALSE)
  # The spatial qualification is a property of the CALL, so it travels in
  # `feat` rather than in the registry.
  if (on("weights") &&
      !.backend_carries_weights(backend, spatial = on("spatial"))) return(FALSE)
  mx <- reg$max_re_terms
  if (!is.null(mx) && isTRUE((feat$n_re_terms %||% 0L) > mx)) return(FALSE)
  TRUE
}


# Does `backend` carry observation weights on a call with (`spatial = TRUE`) or
# without a spatial field? Every log-posterior backend does, through the R
# log-posterior builder; a design backend does when the registry declares
# `carries_weights`, and then non-spatially only -- the spatial Laplace solvers
# have no weight channel. The one predicate the auto selector and the fitter
# argument assembly both read, so what auto picks for a weighted call is what
# dispatch accepts.
#' @keywords internal
.backend_carries_weights <- function(backend, spatial = FALSE) {
  reg <- BACKEND_REGISTRY[[backend]]
  if (is.null(reg)) return(FALSE)
  if (identical(reg$input, "logpost")) return(TRUE)
  isTRUE(reg$carries_weights) && !isTRUE(spatial)
}


# Backends that DETERMINE the random-effect scale rather than conditioning on a
# supplied one.
#
# `?tulpa` says a `sigma_re` handed to one of these is ignored "with a warning",
# and only the EB and random-slope routes warned: agq / gibbs / hmc dropped the
# value silently (gcol33/tulpa#669). Kept as one list, read by the warning, so
# the doc's claim has a single referent.
#' @keywords internal
.re_scale_estimating_backends <- function() {
  c("eb", "agq", "gibbs", "re_cov_nested", "re_cov_gibbs", "re_aghq",
    names(BACKEND_REGISTRY)[vapply(
      BACKEND_REGISTRY,
      function(b) identical(b$input, "modeldata"),
      logical(1))])
}


#' Select best backend within a mode
#' @noRd
select_backend_for_mode <- function(mode, family, n_obs, has_spatial, has_temporal,
                                    has_latent = FALSE, spatial_type = NULL,
                                    has_re = FALSE, feat = list()) {

  if (mode == "exact") {
    # HMC is the default for Exact - most general
    return("hmc")
  }

  if (mode == "structured") {
    # Latent prior blocks need the nested-Laplace path; plain (conditional)
    # Laplace does not consume them. A spatial field is itself a latent Gaussian
    # block whose hyperparameter the designed Tier-2 path integrates, so
    # structured routes the nested-wired types (.NL_FRONTDOOR_NESTED: areal +
    # gp/nngp/hsgp + spde) to nested_laplace too (conditioning at a fixed scale
    # is the explicit mode = "laplace"). SPDE is redirected from nested_laplace
    # to the `spde` backend in tulpa() (its own fit_spde engine). Field types
    # with no nested kernel (multiscale) stay on the conditional Laplace path;
    # `auto` routes those to exact NUTS instead, but an explicit
    # mode = "structured" is honoured as asked. Otherwise Laplace is the Tier 2
    # default.
    if (has_latent) return("nested_laplace")
    if (has_spatial && !is.null(spatial_type) && spatial_type %in% .NL_FRONTDOOR_NESTED) {
      return("nested_laplace")
    }
    # A random-effect term has no dataset-implied scalar sigma_re to condition
    # on within the Structured tier either -- auto's own RE arm already
    # integrates it via re_cov_nested / re_cov_gibbs (see auto_select_mode()),
    # and the Structured tier carries the deterministic (nested-Laplace) half
    # of that pair. Conditioning at sigma_re = 1 is the explicit
    # mode = "laplace" (gcol33/tulpa#787). Guarded on !has_spatial: a spatial
    # field reaching here is one of the non-nested types (multiscale, plain
    # 'car', ...) that the Structured tier deliberately keeps on conditional
    # Laplace -- re_cov_nested carries no spatial field at all.
    if (!has_spatial && has_re && .auto_backend_ok("re_cov_nested", family, feat)) {
      return("re_cov_nested")
    }
    return("laplace")
  }

  if (mode == "optimized") {
    # VI is currently the only Tier 3 backend
    return("vi")
  }

  # Fallback
  return("hmc")
}


# ==============================================================================
# Mode Validation
# ==============================================================================

#' Validate that a fit used the expected mode
#'
#' @description
#' Ensures the fit object was created with the expected inference mode.
#' Useful for enforcing mode requirements in downstream analysis.
#'
#' A fit matches a mode when the mode names either the backend that ran
#' (`fit$backend`, e.g. `"laplace"`, `"mala"`, `"re_cov_gibbs"`) or the tier it
#' ran in (`fit$inference_mode`: `"exact"`, `"structured"`, `"optimized"`). So a
#' fit made with `mode = "laplace"` validates against both `"laplace"` and
#' `"structured"`, and a fit made with `mode = "auto"` validates against the
#' backend and the tier auto chose. A request that `tulpa()` overrode (see
#' `fit$mode_overridden`) validates against what ran, not what was asked for.
#'
#' @param fit A tulpa_fit object
#' @param expected_mode Mode(s) that should have been used: a character vector
#'   of backend and / or tier names, any one of which is accepted.
#' @param error If TRUE (default), error on mismatch. If FALSE, return logical.
#'
#' @return If error = FALSE, returns TRUE/FALSE. Otherwise errors on mismatch.
#'
#' @examples
#' \donttest{
#' set.seed(1)
#' d <- data.frame(x = rnorm(60))
#' d$y <- rpois(60, exp(0.5 + 0.3 * d$x))
#' fit <- tulpa(y ~ x, data = d, family = "poisson", mode = "hmc",
#'              control = list(n_iter = 400L, warmup = 200L, n_chains = 2L,
#'                             seed = 1L))
#' validate_mode(fit, "exact")                      # the tier hmc runs in
#' validate_mode(fit, "structured", error = FALSE)  # FALSE
#' }
#' @export
validate_mode <- function(fit, expected_mode, error = TRUE) {
  if (!inherits(fit, "tulpa_fit")) {
    stop("fit must be a tulpa_fit object", call. = FALSE)
  }
  if (!is.character(expected_mode) || !length(expected_mode) ||
      anyNA(expected_mode)) {
    stop("`expected_mode` must be a non-empty character vector of mode names.",
         call. = FALSE)
  }

  actual_backend <- tolower(fit$backend %||% NA_character_)
  actual_mode    <- tolower(fit$inference_mode %||% NA_character_)
  actual_tier    <- fit$inference_tier %||% NA

  # `fit$inference_mode` holds the TIER name, so comparing against it alone
  # rejected the backend a fit was run with (gcol33/tulpa#879). A name that is
  # neither a tier, a backend, nor this fit's own backend is a typo rather
  # than a mismatch, and says so.
  expected_mode <- tolower(expected_mode)
  known <- c(names(INFERENCE_TIERS), ALL_BACKENDS, actual_backend)
  unknown <- setdiff(expected_mode, known)
  if (length(unknown)) {
    stop(sprintf(paste0(
      "`expected_mode` names unknown mode(s): %s. Use a tier (%s) or a ",
      "backend (%s)."),
      paste0("'", unknown, "'", collapse = ", "),
      paste(names(INFERENCE_TIERS), collapse = ", "),
      paste(ALL_BACKENDS, collapse = ", ")), call. = FALSE)
  }

  expected_tiers <- vapply(intersect(expected_mode, names(INFERENCE_TIERS)),
                           function(m) INFERENCE_TIERS[[m]]$tier, integer(1))
  match <- any(expected_mode %in% c(actual_backend, actual_mode)) ||
    (!is.na(actual_tier) && actual_tier %in% expected_tiers)

  if (!match && isTRUE(error)) {
    stop(sprintf(
      "Mode mismatch: fit used backend '%s' (%s, Tier %s) but expected %s.\n%s",
      fit$backend %||% "unknown",
      fit$inference_mode %||% "unknown tier",
      ifelse(is.na(actual_tier), "?", actual_tier),
      paste0("'", expected_mode, "'", collapse = " or "),
      "Refit the model with the correct mode, or adjust your expectations."
    ), call. = FALSE)
  }

  match
}


# ==============================================================================
# User-Facing Information
# ==============================================================================
#' Print inference mode information
#'
#' @description
#' Displays information about the inference modes available in tulpa,
#' including their tiers, guarantees, and appropriate use cases.
#'
#' @return `NULL`, invisibly. Called for the side effect of printing the
#'   available inference modes, their tiers, guarantees, and use cases to the
#'   console.
#'
#' @examples
#' inference_mode_info()
#' @export
inference_mode_info <- function() {
  cat("tulpa Inference Modes\n")
  cat("========================\n\n")

  cat("Modes encode EPISTEMIC GUARANTEES, not just runtime characteristics.\n\n")

  for (tier_name in c("exact", "structured", "optimized")) {
    tier <- INFERENCE_TIERS[[tier_name]]

    cat(sprintf("TIER %d: %s\n", tier$tier, toupper(tier$name)))
    cat(sprintf("  %s\n", tier$description))
    cat(sprintf("  Guarantee: %s\n", tier$guarantee))
    cat(sprintf("  Backends: %s\n", paste(tier$backends, collapse = ", ")))
    cat(sprintf("  Note: %s\n", tier$note))
    cat("\n")
  }

  cat("AUTO MODE\n")
  cat("  Selects between Tier 1 and Tier 2 only.\n")
  cat("  NEVER silently chooses Tier 3 (Optimized).\n")
  cat("  Contract: 'Use the most reliable method expected to finish.'\n\n")

  cat("Usage (by tier):\n")
  cat("  tulpa(..., mode = 'auto')       # Tier 1 or 2 (default)\n")
  cat("  tulpa(..., mode = 'exact')      # Tier 1\n")
  cat("  tulpa(..., mode = 'structured') # Tier 2\n")
  cat("  tulpa(..., mode = 'optimized')  # Tier 3 (explicit opt-in)\n\n")

  cat(paste0(
    "Backends ([R] = tulpa(mode =) dispatch, [direct] = has an R fitter but ",
    "tulpa(mode =) refuses it -- call the fitter itself, [C-ABI] = ",
    "model-package kernel only):\n"))
  for (tk in c("exact", "structured", "optimized")) {
    for (b in .tier_backends(tk)) {
      entry <- BACKEND_REGISTRY[[b]]
      if (is.null(entry$fitter)) {
        tag <- "[C-ABI]"
        usage <- sprintf("mode = '%s'", b)
      } else if (isTRUE(entry$dispatchable %||% TRUE)) {
        tag <- "[R]    "
        usage <- sprintf("mode = '%s'", b)
      } else {
        tag <- "[direct]"
        usage <- sprintf("%s()", entry$fitter)
      }
      cat(sprintf("  %s %-26s # Tier %d (%s)\n",
                  tag, usage, TIER_META[[tk]]$tier, TIER_META[[tk]]$name))
    }
  }

  invisible(NULL)
}