gam_solve/estimate/
external_options.rs

1use super::*;
2
3pub struct ExternalOptimResult {
4    pub beta: Array1<f64>,
5    /// Canonical optimized log-strengths. Physical `lambdas` are derived from
6    /// this vector through the shared exact-domain conversion.
7    pub log_lambdas: Array1<f64>,
8    pub lambdas: Array1<f64>,
9    pub likelihood_family: LikelihoodSpec,
10    pub likelihood_scale: LikelihoodScaleMetadata,
11    pub log_likelihood_normalization: LogLikelihoodNormalization,
12    pub log_likelihood: f64,
13    /// Residual scale on the response scale.
14    ///
15    /// Contract: Gaussian identity models store the residual standard
16    /// deviation sigma here. Non-Gaussian families keep the response-scale
17    /// summary used by their explicit likelihood-scale metadata.
18    pub standard_deviation: f64,
19    pub iterations: usize,
20    pub finalgrad_norm: f64,
21    /// True iff the outer optimizer reached a stationary point (gradient
22    /// norm below tolerance), as reported by the optimizer itself. False
23    /// when the run exhausted its iteration budget without reaching the
24    /// gradient tolerance. Downstream consumers should NOT assume that a
25    /// fit with `outer_converged == false` is unusable — it may still be
26    /// the best basin reached given the budget — but they must not treat
27    /// it as certified-converged either.
28    pub outer_converged: bool,
29    pub pirls_status: crate::pirls::PirlsStatus,
30    pub deviance: f64,
31    /// Stable quadratic penalty term βᵀSβ, including any solver ridge quadratic.
32    pub stable_penalty_term: f64,
33    pub used_device: bool,
34    pub max_abs_eta: f64,
35    pub constraint_kkt: Option<crate::pirls::ConstraintKktDiagnostics>,
36    pub artifacts: FitArtifacts,
37    pub geometry: Option<FitGeometry>,
38    pub inference: Option<FitInference>,
39    /// Complete REML/LAML objective value used for smoothing selection, or
40    /// `None` when the converged fit sits on the zero-dispersion Gaussian
41    /// boundary and therefore has no finite criterion value at all. Same
42    /// contract as [`crate::estimate::UnifiedFitResult::reml_score`]; the
43    /// boundary is DETECTED here, where the dispersion is actually estimated,
44    /// rather than predicted by an entry-level shape rule (#2595).
45    pub reml_score: Option<f64>,
46    pub fitted_link: FittedLinkState,
47    /// Number of outer REML cost-only evaluations executed during the fit
48    /// (the count the outer optimizer's trust-region/line-search probes drive,
49    /// each paying an inner P-IRLS solve). Surfaced for regression guards on
50    /// outer work (#1575); not part of the statistical contract.
51    pub outer_cost_evals: usize,
52    /// Number of *actual* full-n inner P-IRLS solves performed (cache-missing
53    /// `prepare_eval_bundlewithkey` calls). This is the true cost driver the
54    /// #1575 slowdown is measured in ("~150 outer cost evals each running a
55    /// full n-sized P-IRLS"): unlike `outer_cost_evals`, it excludes single-slot
56    /// cache hits and prior short-circuits and includes every solve done by the
57    /// seed-grid prepass, screening, multistart, and finalize phases. Surfaced
58    /// for a regression guard that pins the warm-start / parsimony-waiver /
59    /// PSIS-optin economy (#1575); not part of the statistical contract.
60    pub inner_pirls_solves: usize,
61}
62
63#[derive(Clone)]
64pub struct ExternalOptimOptions {
65    pub family: gam_problem::LikelihoodSpec,
66    pub latent_cloglog: Option<LatentCLogLogState>,
67    pub mixture_link: Option<MixtureLinkSpec>,
68    pub optimize_mixture: bool,
69    pub sas_link: Option<SasLinkSpec>,
70    pub optimize_sas: bool,
71    pub compute_inference: bool,
72    /// Internal lifecycle knob for fits whose result will be immediately
73    /// superseded. Keeps ordinary inference work but skips the live-objective
74    /// rho posterior certificate/escalation until the returned model is known.
75    pub skip_rho_posterior_inference: bool,
76    pub max_iter: usize,
77    pub tol: f64,
78    pub nullspace_dims: Vec<usize>,
79    pub linear_constraints: Option<crate::pirls::LinearInequalityConstraints>,
80    /// Optional explicit Firth override for external fitting families that
81    /// support Jeffreys/Firth bias reduction.
82    /// - `Some(true)`: force Firth on
83    /// - `Some(false)`: force Firth off
84    /// - `None`: use family default behavior
85    pub firth_bias_reduction: Option<bool>,
86    /// Fixed prior on smoothing parameters for explicit joint HMC sampling
87    /// flows. Standard fitting stays on the REML/Laplace path.
88    pub rho_prior: gam_problem::RhoPrior,
89    /// Kronecker-factored penalty system for tensor-product smooth terms.
90    pub kronecker_penalty_system: Option<gam_terms::smooth::KroneckerPenaltySystem>,
91    /// Full Kronecker factored basis for P-IRLS factored reparameterization.
92    pub kronecker_factored: Option<gam_terms::basis::KroneckerFactoredBasis>,
93    /// Explicit cross-process warm-start capability for this fit. `None` is
94    /// disk-silent; clones share one caller-configured store handle.
95    pub persistent_warm_start_store: Option<gam_runtime::warm_start::ConfiguredWarmStartStore>,
96}
97
98pub(crate) fn resolve_external_family(
99    family: &gam_problem::LikelihoodSpec,
100    firth_override: Option<bool>,
101) -> Result<(GlmLikelihoodSpec, bool), EstimationError> {
102    let external_glm_supported = match (&family.response, family.link_function()) {
103        (ResponseFamily::Gaussian, LinkFunction::Identity)
104        | (ResponseFamily::Poisson, LinkFunction::Log)
105        | (ResponseFamily::Gamma, LinkFunction::Log)
106        | (ResponseFamily::Tweedie { .. }, LinkFunction::Log)
107        | (ResponseFamily::NegativeBinomial { .. }, LinkFunction::Log)
108        | (ResponseFamily::Binomial, LinkFunction::Logit)
109        | (ResponseFamily::Binomial, LinkFunction::Probit)
110        | (ResponseFamily::Binomial, LinkFunction::CLogLog)
111        // LogLog and Cauchit are ordinary state-less probability links: they
112        // narrow into `StandardLink` (gam-spec), carry a full 5-jet Fisher
113        // weight (`fisher_weight_jet5` → `component_fisher_weight_jet5` for
114        // LinkComponent::{LogLog,Cauchit}), and their inverse-link jets live
115        // in `mixture_link.rs` — the exact same external-design/P-IRLS
116        // machinery probit/cloglog ride. #2104 un-gated them at validation
117        // (`link_legal_for_family`); this is the fitting half of that wiring.
118        | (ResponseFamily::Binomial, LinkFunction::LogLog)
119        | (ResponseFamily::Binomial, LinkFunction::Cauchit)
120        | (ResponseFamily::Binomial, LinkFunction::Sas)
121        | (ResponseFamily::Binomial, LinkFunction::BetaLogistic) => true,
122        // Beta regression with a constant precision φ is a genuine-dispersion
123        // mean family on par with Gamma/Tweedie/Negative-Binomial: the inner
124        // P-IRLS carries its full fixed-φ Fisher information and the outer loop
125        // estimates φ by the Pearson moment estimator (`estimate_beta_phi_from_eta`,
126        // mirroring the Tweedie φ / Gamma shape / NegBin θ locks). A
127        // `noise_formula` upgrades it to a dispersion-location-scale model that
128        // smooths log φ; without one, the external GLM route fits the mean with
129        // a single estimated φ exactly as betareg does by default.
130        (ResponseFamily::Beta { .. }, LinkFunction::Logit) => true,
131        _ => false,
132    };
133    if !external_glm_supported {
134        crate::bail_invalid_estim!(
135            "optimize_external_design requires a supported standard GLM family/link; got {}. \
136             The external-design route supports Gaussian(identity), Binomial(logit/probit/cloglog/loglog/cauchit/SAS/Beta-Logistic), \
137             Beta(logit), and Poisson/Gamma/Tweedie/Negative-Binomial(log). For Beta precision modeling \
138             add a noise_formula to upgrade to the dispersion-location-scale route",
139            family.pretty_name(),
140        );
141    }
142
143    let supports_firth = family.supports_firth();
144    if firth_override == Some(true) && !supports_firth {
145        crate::bail_invalid_estim!(
146            "firth_bias_reduction requires a Binomial inverse link with a Fisher-weight jet; {} does not support it",
147            family.pretty_name(),
148        );
149    }
150
151    if let ResponseFamily::Tweedie { p } = &family.response {
152        if !gam_problem::is_valid_tweedie_power(*p) {
153            crate::bail_invalid_estim!("optimize_external_design requires a GLM family; Tweedie variance power must be finite and strictly between 1 and 2; use PoissonLog or GammaLog for boundary cases"
154                    .to_string(),);
155        }
156    }
157    Ok((
158        GlmLikelihoodSpec::canonical(family.clone()),
159        firth_override.unwrap_or(false) && supports_firth,
160    ))
161}
162
163#[inline]
164pub(crate) fn effective_sas_link_for_family(
165    family: &gam_problem::LikelihoodSpec,
166    sas_link: Option<SasLinkSpec>,
167) -> Option<SasLinkSpec> {
168    if (family.is_binomial_sas() || family.is_binomial_beta_logistic()) && sas_link.is_none() {
169        Some(SasLinkSpec {
170            initial_epsilon: 0.0,
171            initial_log_delta: 0.0,
172        })
173    } else {
174        sas_link
175    }
176}
177
178#[inline]
179pub(crate) fn resolved_external_inverse_link(
180    link: LinkFunction,
181    latent_cloglog: Option<LatentCLogLogState>,
182    mixture_link: Option<&MixtureLinkSpec>,
183    sas_link: Option<SasLinkSpec>,
184) -> Result<InverseLink, EstimationError> {
185    if let Some(state) = latent_cloglog {
186        return Ok(InverseLink::LatentCLogLog(state));
187    }
188    if let Some(spec) = mixture_link {
189        return Ok(InverseLink::Mixture(state_fromspec(spec).map_err(|e| {
190            EstimationError::InvalidInput(format!("invalid blended inverse link: {e}"))
191        })?));
192    }
193    if let Some(spec) = sas_link {
194        return Ok(match link {
195            LinkFunction::BetaLogistic => {
196                InverseLink::BetaLogistic(state_from_beta_logisticspec(spec).map_err(|e| {
197                    EstimationError::InvalidInput(format!("invalid Beta-Logistic link: {e}"))
198                })?)
199            }
200            _ => InverseLink::Sas(
201                state_from_sasspec(spec)
202                    .map_err(|e| EstimationError::InvalidInput(format!("invalid SAS link: {e}")))?,
203            ),
204        });
205    }
206    Ok(InverseLink::Standard(StandardLink::try_from(link).map_err(|e| {
207        EstimationError::InvalidInput(format!(
208            "inverse link resolution: {e}; supply `sas_link` or `latent_cloglog` configuration for state-bearing links"
209        ))
210    })?))
211}
212
213#[inline]
214pub(crate) fn resolved_external_config(
215    opts: &ExternalOptimOptions,
216) -> Result<(RemlConfig, Option<SasLinkSpec>), EstimationError> {
217    if opts.latent_cloglog.is_some() && (opts.mixture_link.is_some() || opts.sas_link.is_some()) {
218        crate::bail_invalid_estim!(
219            "latent_cloglog cannot be combined with mixture_link or sas_link"
220        );
221    }
222    if opts.mixture_link.is_some() && opts.sas_link.is_some() {
223        crate::bail_invalid_estim!("mixture_link and sas_link are mutually exclusive");
224    }
225    if opts.family.is_latent_cloglog() && opts.latent_cloglog.is_none() {
226        crate::bail_invalid_estim!("BinomialLatentCLogLog requires latent_cloglog state");
227    }
228    if opts.latent_cloglog.is_some() && !opts.family.is_latent_cloglog() {
229        crate::bail_invalid_estim!("latent_cloglog is only supported with BinomialLatentCLogLog");
230    }
231    let effective_sas_link = effective_sas_link_for_family(&opts.family, opts.sas_link);
232    let (likelihood, firth_active) =
233        resolve_external_family(&opts.family, opts.firth_bias_reduction)?;
234    let link = likelihood.link_function();
235    let mut cfg = RemlConfig::external(likelihood, opts.tol, firth_active);
236    cfg.link_kind = resolved_external_inverse_link(
237        link,
238        opts.latent_cloglog,
239        opts.mixture_link.as_ref(),
240        effective_sas_link,
241    )?;
242    Ok((cfg, effective_sas_link))
243}
244
245/// Shape/bounds validation for a single [`PenaltySpec`] against the total
246/// coefficient width `p`. The checks live in `gam-terms` beside the
247/// `PenaltySpec` definition they validate; re-exported here so `estimate`'s
248/// existing paths keep resolving and both crates emit identical diagnostics.
249pub(crate) use gam_terms::validate_penalty_spec_shape;
gam_solve/estimate/external_options.rs

gam_solve/estimate/
external_options.rs
