Skip to main content

gam_solve/
model_types.rs

1//! # Shared model-estimation contract types
2//!
3//! Lower-layer types that both the `families` layer (which constructs penalty
4//! and dispersion specifications and propagates estimation errors) and the
5//! `solver` layer (which consumes them) need to name. Hosting them here breaks
6//! the `families → solver::estimate` back-edge that #1135 tracks: families now
7//! import these from `crate::model_types` instead of reaching *up* into
8//! `crate::solver::estimate`.
9//!
10//! ## Layering
11//! These types depend only on lower or sibling layers (`linalg`, `terms`,
12//! `families`' error types) — never on `solver`. `EstimationError` carries
13//! `#[from]` conversions for the family error types it wraps, which is the
14//! allowed downward direction.
15
16use ndarray::{Array1, Array2};
17
18pub use gam_problem::EstimationError;
19
20// ===========================================================================
21// Dispersion
22// ===========================================================================
23
24/// Dispersion contract used by inferential covariance and reference distributions.
25///
26/// This type lives in `gam-problem`; re-exported here so model result APIs and
27/// existing engine code name the same neutral contract.
28pub use gam_problem::Dispersion;
29
30// ===========================================================================
31// Constraint/KKT carriers
32// ===========================================================================
33
34/// Active row block of the joint linear inequality constraint matrix at the
35/// converged inner iterate.
36#[derive(Clone, Debug)]
37pub struct ActiveLinearConstraintBlock {
38    /// `k_active x p` matrix of active constraint rows.
39    pub a: Array2<f64>,
40}
41
42/// Subspace represented by a stored KKT residual.
43#[derive(Clone, Copy, Debug, PartialEq, Eq)]
44pub enum KktResidualSubspace {
45    /// Residual after active-constraint normal components have been stripped:
46    /// `r_A = P_T(Sbeta + Gamma beta - grad ell)`.
47    ActiveProjected,
48    /// Residual additionally projected into the retained identifiable range:
49    /// `r_R = R R^T r_A`.
50    ReducedRange,
51}
52
53/// KKT residual `r = grad_beta L_pen(beta_hat)` at the converged inner
54/// iterate, tagged with the exact represented subspace.
55#[derive(Clone, Debug)]
56pub struct ProjectedKktResidual {
57    /// The residual vector in the full coefficient coordinates. Active and
58    /// reduced-range projection zero out excluded directions rather than
59    /// shortening the vector, so its length remains `p`.
60    pub(crate) residual: Array1<f64>,
61    pub(crate) subspace: KktResidualSubspace,
62    /// The KKT-stationarity tolerance the inner solver compared the residual
63    /// against when the certificate fired.
64    pub(crate) residual_tol: Option<f64>,
65    /// `total_p - active_set_size` at the producing iterate.
66    pub(crate) free_rank: Option<usize>,
67}
68
69impl ProjectedKktResidual {
70    /// Construct from `r_A = P_T(Sbeta + Gamma beta - grad ell)`, with active
71    /// constraint multipliers removed but before any reduced-range projection.
72    pub fn from_active_projected(residual: Array1<f64>) -> Self {
73        Self {
74            residual,
75            subspace: KktResidualSubspace::ActiveProjected,
76            residual_tol: None,
77            free_rank: None,
78        }
79    }
80
81    /// Construct from `r_R = R R^T r_A`, where `R` is the actual reduced
82    /// identifiable basis used by the projected inverse kernel.
83    pub(crate) fn from_reduced_range(residual: Array1<f64>) -> Self {
84        Self {
85            residual,
86            subspace: KktResidualSubspace::ReducedRange,
87            residual_tol: None,
88            free_rank: None,
89        }
90    }
91
92    /// Attach the KKT tolerance and free-subspace rank to a previously
93    /// constructed residual.
94    pub fn with_metadata(mut self, residual_tol: f64, free_rank: usize) -> Self {
95        self.residual_tol = Some(residual_tol);
96        self.free_rank = Some(free_rank);
97        self
98    }
99
100    /// Attach the free-subspace rank alone, for producers that know which
101    /// subspace they projected onto but have no inner-solver tolerance to
102    /// report.
103    ///
104    /// `with_metadata` demands both numbers, and the standard REML assembly
105    /// builders have only one: the active set is empty by construction there
106    /// (see the guards at both call sites), so `free_rank` is exactly the
107    /// coefficient count, while the inner solve exposes no KKT tolerance at
108    /// that seam. Before this existed those producers passed neither, and a
109    /// residual that does not record its own projection cannot be audited --
110    /// `subspace=ActiveProjected free_rank=None` says "some directions were
111    /// removed" without saying which or how many. Supplying a made-up
112    /// tolerance just to reach `with_metadata` would be worse: a number where
113    /// a measurement belongs.
114    pub fn with_free_rank(mut self, free_rank: usize) -> Self {
115        self.free_rank = Some(free_rank);
116        self
117    }
118
119    /// Borrow the underlying free-space residual for the H^-1*r solve and its
120    /// rho-derivatives.
121    pub fn as_array(&self) -> &Array1<f64> {
122        &self.residual
123    }
124
125    pub fn subspace(&self) -> KktResidualSubspace {
126        self.subspace
127    }
128
129    /// Sup-norm of the projected residual — the scalar a stationarity verdict is
130    /// actually taken on. A refusal that reports only its cycle count cannot say
131    /// whether the solve was one order from the tolerance or ten, and those have
132    /// different causes; this is the number to put next to the count.
133    pub fn inf_norm(&self) -> f64 {
134        self.residual
135            .iter()
136            .fold(0.0_f64, |acc, value| acc.max(value.abs()))
137    }
138}
139
140// ===========================================================================
141// CoefficientPriorMean + PenaltySpec
142// ===========================================================================
143
144/// Programmatic prior mean for a coefficient penalty block.
145///
146/// This type now lives in the neutral `gam-problem` crate (with its inherent
147/// `evaluate` returning `gam_problem::PriorMeanError`); re-exported here so all
148/// existing `crate::estimate::CoefficientPriorMean` references keep resolving.
149/// Solver-side callers map `PriorMeanError` into `EstimationError::InvalidInput`.
150pub use gam_problem::CoefficientPriorMean;
151
152pub use gam_terms::penalty_spec::PenaltySpec;
153
154pub(crate) mod result_types;
155
156pub(crate) use result_types::CERTIFICATE_RAIL_MARGIN;
157pub use result_types::{
158    BlockRole, CertifiedRung, CoefficientCovarianceDefinition,
159    CovarianceDeclined, CurvatureAdmissibility, CurvatureEvidence, CurvatureFloorClearance,
160    DisplayCoefficientUncertainty,
161    ExecutionPath, FitArtifacts, FitGeometry, FitInference, FitOptions, FittedBlock,
162    FittedLinkState, InferenceCovarianceMode, NO_CRITERION_AT_EXACT_FIT, OuterCriterionCertificate,
163    OuterStationarityCertificate, RailCoordinate, RailTailEvidence, RailedCoordinateFact,
164    SmoothingCorrectionMethod, UnifiedFitResult, UnifiedFitResultParts, WorkingGeometry,
165    criterion_display, ensure_finite_scalar, is_zero_dispersion_boundary,
166    saved_latent_cloglog_state_from_fit, saved_mixture_state_from_fit, saved_sas_state_from_fit,
167    validate_all_finite, validate_dense_hessian_export,
168    validate_explicit_dense_hessian_for_whitening,
169};