Skip to main content

gam_problem/
types.rs

1use ndarray::{Array1, ArrayView1};
2use serde::{Deserialize, Serialize};
3use std::ops::{Deref, DerefMut};
4
5pub use gam_linalg::RidgePolicy;
6
7pub use gam_spec::*;
8
9/// Storage form of the ridge penalty matrix.
10#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
11pub enum RidgeMatrixForm {
12    /// Ridge matrix is `delta * I`.
13    ScaledIdentity,
14}
15
16#[derive(Debug, Clone, PartialEq, Eq)]
17pub struct InvalidStabilization {
18    reason: String,
19}
20
21impl InvalidStabilization {
22    fn new(reason: impl Into<String>) -> Self {
23        Self {
24            reason: reason.into(),
25        }
26    }
27
28    pub fn reason(&self) -> &str {
29        &self.reason
30    }
31}
32
33impl std::fmt::Display for InvalidStabilization {
34    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
35        write!(f, "invalid stabilization metadata: {}", self.reason)
36    }
37}
38
39impl std::error::Error for InvalidStabilization {}
40
41#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
42struct RidgePassportWire {
43    delta: f64,
44    matrix_form: RidgeMatrixForm,
45    policy: RidgePolicy,
46}
47
48/// Validated ridge metadata stamped into a fitted PIRLS result.
49///
50/// Construction and deserialization both reject non-finite or negative
51/// magnitudes; fields are private so invalid state cannot be assembled with a
52/// literal.
53#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
54#[serde(try_from = "RidgePassportWire", into = "RidgePassportWire")]
55pub struct RidgePassport {
56    delta: f64,
57    matrix_form: RidgeMatrixForm,
58    policy: RidgePolicy,
59}
60
61impl RidgePassport {
62    pub fn scaled_identity(delta: f64, policy: RidgePolicy) -> Result<Self, InvalidStabilization> {
63        if !(delta.is_finite() && delta >= 0.0) {
64            return Err(InvalidStabilization::new(format!(
65                "ridge delta must be finite and non-negative, got {delta:?}"
66            )));
67        }
68        Ok(Self {
69            delta: if delta == 0.0 { 0.0 } else { delta },
70            matrix_form: RidgeMatrixForm::ScaledIdentity,
71            policy,
72        })
73    }
74
75    /// Exact zero-ridge passport; this fixed sentinel has no unchecked input.
76    pub const fn zero(policy: RidgePolicy) -> Self {
77        Self {
78            delta: 0.0,
79            matrix_form: RidgeMatrixForm::ScaledIdentity,
80            policy,
81        }
82    }
83
84    #[inline]
85    pub const fn delta(self) -> f64 {
86        self.delta
87    }
88
89    #[inline]
90    pub const fn matrix_form(self) -> RidgeMatrixForm {
91        self.matrix_form
92    }
93
94    #[inline]
95    pub const fn policy(self) -> RidgePolicy {
96        self.policy
97    }
98
99    #[inline]
100    pub const fn penalty_logdet_ridge(self) -> f64 {
101        if self.policy.accounts_for_objective() {
102            self.delta
103        } else {
104            0.0
105        }
106    }
107
108}
109
110impl TryFrom<RidgePassportWire> for RidgePassport {
111    type Error = InvalidStabilization;
112
113    fn try_from(wire: RidgePassportWire) -> Result<Self, Self::Error> {
114        let mut passport = Self::scaled_identity(wire.delta, wire.policy)?;
115        passport.matrix_form = wire.matrix_form;
116        Ok(passport)
117    }
118}
119
120impl From<RidgePassport> for RidgePassportWire {
121    fn from(passport: RidgePassport) -> Self {
122        Self {
123            delta: passport.delta,
124            matrix_form: passport.matrix_form,
125            policy: passport.policy,
126        }
127    }
128}
129
130// ============================================================================
131// StabilizationLedger: canonical accounting for every fixed/heuristic ridge
132// added anywhere in the solver, linear-algebra, or family code paths.
133//
134// Five semantically distinct ridge uses must NEVER be conflated:
135//   1. SolverDampingOnly      — Levenberg/trust-region damping; never enters
136//                               objective, gradient, logdet, Hessian, or any
137//                               saved/serialized model artifact.
138//   2. NumericalPerturbation  — added strictly so a linear solve is well-
139//                               posed (e.g. Cholesky of a near-singular
140//                               matrix). Carries an optional backward-error
141//                               bound. Does NOT change the objective.
142//   3. ExplicitPrior          — model-level `delta * I` (or block-diagonal)
143//                               prior. Appears in quadratic, log normalizer,
144//                               Laplace Hessian, serialization, diagnostics.
145//   4. ApproximationOnly      — changes a named downstream approximation
146//                               (for example sigma-point cubature covariance)
147//                               but not the fitted model or its objective.
148//   5. ObjectiveStabilization — algorithm-selected ridge consistently included
149//                               in the fitted objective, preserving exact versus
150//                               approximate determinant provenance.
151//
152// `RidgePassport` above already encodes the inclusion-flag matrix for the
153// PIRLS Laplace ridge specifically; this ledger is the broader sibling for
154// every declared solver, approximation, and model ridge, so a downstream consumer can ask
155// `ledger.quadratic_delta()` rather than rediscovering the policy. The three
156// inclusion bits were lifted into the `StabilizationKind` discriminant so the
157// (kind, inclusion-flags) invariant is enforced statically — heterogeneous
158// combinations like "ExplicitPrior with quadratic excluded" no longer typecheck.
159// ============================================================================
160
161/// Inertia of a symmetric matrix (count of positive / zero / negative
162/// eigenvalues). Used by `bump_with_matrix` and other indefinite-aware
163/// stabilization rules to drive δ from spectral evidence rather than a
164/// condition-number heuristic.
165#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
166#[serde(try_from = "InertiaWire", into = "InertiaWire")]
167pub struct Inertia {
168    positive: usize,
169    zero: usize,
170    negative: usize,
171}
172
173#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
174struct InertiaWire {
175    positive: usize,
176    zero: usize,
177    negative: usize,
178}
179
180impl Inertia {
181    pub fn new(
182        positive: usize,
183        zero: usize,
184        negative: usize,
185    ) -> Result<Self, InvalidStabilization> {
186        let total = positive
187            .checked_add(zero)
188            .and_then(|value| value.checked_add(negative))
189            .ok_or_else(|| InvalidStabilization::new("inertia count sum overflows usize"))?;
190        if total == 0 {
191            return Err(InvalidStabilization::new(
192                "inertia must describe a non-empty matrix",
193            ));
194        }
195        Ok(Self {
196            positive,
197            zero,
198            negative,
199        })
200    }
201
202    pub const fn positive(self) -> usize {
203        self.positive
204    }
205
206    pub const fn zero(self) -> usize {
207        self.zero
208    }
209
210    pub const fn negative(self) -> usize {
211        self.negative
212    }
213
214    pub fn total(self) -> usize {
215        self.positive + self.zero + self.negative
216    }
217}
218
219impl TryFrom<InertiaWire> for Inertia {
220    type Error = InvalidStabilization;
221
222    fn try_from(wire: InertiaWire) -> Result<Self, Self::Error> {
223        Self::new(wire.positive, wire.zero, wire.negative)
224    }
225}
226
227impl From<Inertia> for InertiaWire {
228    fn from(inertia: Inertia) -> Self {
229        Self {
230            positive: inertia.positive,
231            zero: inertia.zero,
232            negative: inertia.negative,
233        }
234    }
235}
236
237/// Why a stabilization δ was chosen at this site.
238#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
239pub enum StabilizationRule {
240    /// δ is a hard-coded constant in the source.
241    FixedConstant,
242    /// δ chosen so the SPD floor τ is met: δ = max(0, τ - λ_min(H)).
243    InertiaTarget { spd_floor: f64 },
244    /// δ chosen via a condition-number / sqrt-ratio heuristic.
245    Heuristic,
246    /// User- or family-specified prior precision.
247    UserSpecified,
248    /// δ derived from a back-off escalation after a factorization failure.
249    BackoffEscalation { attempts: usize },
250}
251
252/// Semantically distinct flavours a ridge δ can have.
253#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
254pub enum StabilizationKind {
255    None,
256    /// LM/TR damping. NEVER enters the objective, gradient, logdet, Hessian,
257    /// or any saved model artifact. Lives only inside the trust-region step.
258    SolverDampingOnly,
259    /// Added strictly so a linear solve succeeds. The objective/Hessian the
260    /// caller sees is unchanged; the perturbation is a property of the
261    /// solver, not the model. Its optional backward-error bound lives on the
262    /// enclosing ledger.
263    NumericalPerturbation,
264    /// An explicit part of a downstream approximation, not of the fitted
265    /// model. Unlike `NumericalPerturbation`, consumers must not report the
266    /// result as if the unperturbed estimand had been evaluated.
267    ApproximationOnly,
268    /// Algorithm-selected ridge consistently included in the fitted objective.
269    /// The ledger's objective policy preserves determinant provenance.
270    ObjectiveStabilization,
271    /// Part of the model. Enters quadratic, log normalizer, Hessian,
272    /// serialization, and user-visible summaries.
273    ExplicitPrior,
274}
275
276#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
277struct StabilizationLedgerWire {
278    kind: StabilizationKind,
279    delta: f64,
280    matrix_form: RidgeMatrixForm,
281    chosen_by: StabilizationRule,
282    objective_policy: Option<RidgePolicy>,
283    backward_error_bound: Option<f64>,
284    inertia_before: Option<Inertia>,
285    inertia_after: Option<Inertia>,
286}
287
288/// Canonical validated record of one stabilization applied at one site.
289#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
290#[serde(try_from = "StabilizationLedgerWire", into = "StabilizationLedgerWire")]
291pub struct StabilizationLedger {
292    kind: StabilizationKind,
293    delta: f64,
294    matrix_form: RidgeMatrixForm,
295    chosen_by: StabilizationRule,
296    objective_policy: Option<RidgePolicy>,
297    backward_error_bound: Option<f64>,
298    inertia_before: Option<Inertia>,
299    inertia_after: Option<Inertia>,
300}
301
302impl StabilizationLedger {
303    /// "No stabilization applied at this site" sentinel.
304    pub const fn none() -> Self {
305        Self {
306            kind: StabilizationKind::None,
307            delta: 0.0,
308            matrix_form: RidgeMatrixForm::ScaledIdentity,
309            chosen_by: StabilizationRule::FixedConstant,
310            objective_policy: None,
311            backward_error_bound: None,
312            inertia_before: None,
313            inertia_after: None,
314        }
315    }
316
317    fn try_new(
318        kind: StabilizationKind,
319        delta: f64,
320        chosen_by: StabilizationRule,
321        backward_error_bound: Option<f64>,
322    ) -> Result<Self, InvalidStabilization> {
323        if matches!(kind, StabilizationKind::None) {
324            return Err(InvalidStabilization::new(
325                "None stabilization must be constructed with StabilizationLedger::none",
326            ));
327        }
328        if !(delta.is_finite() && delta >= 0.0) {
329            return Err(InvalidStabilization::new(format!(
330                "stabilization delta must be finite and non-negative, got {delta:?}"
331            )));
332        }
333        Self::validate_rule(chosen_by)?;
334        if let Some(bound) = backward_error_bound
335            && !(bound.is_finite() && bound >= 0.0)
336        {
337            return Err(InvalidStabilization::new(format!(
338                "backward-error bound must be finite and non-negative, got {bound:?}"
339            )));
340        }
341        if !matches!(kind, StabilizationKind::NumericalPerturbation)
342            && backward_error_bound.is_some()
343        {
344            return Err(InvalidStabilization::new(
345                "only a numerical perturbation may carry a backward-error bound",
346            ));
347        }
348        Ok(Self {
349            kind,
350            delta: if delta == 0.0 { 0.0 } else { delta },
351            matrix_form: RidgeMatrixForm::ScaledIdentity,
352            chosen_by,
353            objective_policy: None,
354            backward_error_bound,
355            inertia_before: None,
356            inertia_after: None,
357        })
358    }
359
360    fn validate_rule(rule: StabilizationRule) -> Result<(), InvalidStabilization> {
361        match rule {
362            StabilizationRule::InertiaTarget { spd_floor }
363                if !(spd_floor.is_finite() && spd_floor > 0.0) =>
364            {
365                Err(InvalidStabilization::new(format!(
366                    "inertia-target SPD floor must be finite and strictly positive, got {spd_floor:?}"
367                )))
368            }
369            StabilizationRule::BackoffEscalation { attempts } if attempts == 0 => Err(
370                InvalidStabilization::new("backoff escalation must record at least one attempt"),
371            ),
372            _ => Ok(()),
373        }
374    }
375
376    pub fn with_inertia(
377        mut self,
378        before: Option<Inertia>,
379        after: Option<Inertia>,
380    ) -> Result<Self, InvalidStabilization> {
381        if before.is_some() != after.is_some() {
382            return Err(InvalidStabilization::new(
383                "inertia diagnostics must record both the pre- and post-stabilization matrix",
384            ));
385        }
386        if let (Some(before), Some(after)) = (before, after)
387            && before.total() != after.total()
388        {
389            return Err(InvalidStabilization::new(format!(
390                "inertia dimensions disagree: before={}, after={}",
391                before.total(),
392                after.total()
393            )));
394        }
395        if matches!(self.chosen_by, StabilizationRule::InertiaTarget { .. }) {
396            let Some(after) = after else {
397                return Err(InvalidStabilization::new(
398                    "inertia-target stabilization must record post-stabilization inertia",
399                ));
400            };
401            if after.zero() != 0 || after.negative() != 0 {
402                return Err(InvalidStabilization::new(format!(
403                    "inertia-target stabilization did not certify SPD curvature: zero={}, negative={}",
404                    after.zero(),
405                    after.negative()
406                )));
407            }
408        }
409        self.inertia_before = before;
410        self.inertia_after = after;
411        Ok(self)
412    }
413
414    pub const fn kind(self) -> StabilizationKind {
415        self.kind
416    }
417
418    pub const fn delta(self) -> f64 {
419        self.delta
420    }
421
422    pub const fn matrix_form(self) -> RidgeMatrixForm {
423        self.matrix_form
424    }
425
426    pub const fn chosen_by(self) -> StabilizationRule {
427        self.chosen_by
428    }
429
430    /// Exact determinant/objective provenance for an explicit prior or
431    /// objective-accounted algorithmic stabilization. `None` for every
432    /// solver-only, numerical, and approximation-only perturbation.
433    pub const fn objective_policy(self) -> Option<RidgePolicy> {
434        self.objective_policy
435    }
436
437    pub const fn backward_error_bound(self) -> Option<f64> {
438        self.backward_error_bound
439    }
440
441    pub const fn inertia_before(self) -> Option<Inertia> {
442        self.inertia_before
443    }
444
445    pub const fn inertia_after(self) -> Option<Inertia> {
446        self.inertia_after
447    }
448
449}
450
451impl TryFrom<StabilizationLedgerWire> for StabilizationLedger {
452    type Error = InvalidStabilization;
453
454    fn try_from(wire: StabilizationLedgerWire) -> Result<Self, Self::Error> {
455        if matches!(wire.kind, StabilizationKind::None) {
456            if wire.delta != 0.0
457                || wire.chosen_by != StabilizationRule::FixedConstant
458                || wire.objective_policy.is_some()
459                || wire.backward_error_bound.is_some()
460                || wire.inertia_before.is_some()
461                || wire.inertia_after.is_some()
462            {
463                return Err(InvalidStabilization::new(
464                    "None stabilization must have zero delta and no diagnostic payload",
465                ));
466            }
467            return Ok(Self::none());
468        }
469        let mut ledger = Self::try_new(
470            wire.kind,
471            wire.delta,
472            wire.chosen_by,
473            wire.backward_error_bound,
474        )?;
475        if matches!(wire.kind, StabilizationKind::ExplicitPrior)
476            && wire.chosen_by != StabilizationRule::UserSpecified
477        {
478            return Err(InvalidStabilization::new(
479                "an explicit prior must be recorded as user specified",
480            ));
481        }
482        match (wire.kind, wire.objective_policy) {
483            (
484                StabilizationKind::ExplicitPrior | StabilizationKind::ObjectiveStabilization,
485                Some(policy),
486            ) if policy.accounts_for_objective() => {
487                ledger.objective_policy = Some(policy);
488            }
489            (StabilizationKind::ExplicitPrior | StabilizationKind::ObjectiveStabilization, _) => {
490                return Err(InvalidStabilization::new(
491                    "objective-accounted stabilization must preserve its ridge policy",
492                ));
493            }
494            (_, Some(_)) => {
495                return Err(InvalidStabilization::new(
496                    "only objective-accounted stabilization may carry objective ridge provenance",
497                ));
498            }
499            (_, None) => {}
500        }
501        ledger.matrix_form = wire.matrix_form;
502        ledger.with_inertia(wire.inertia_before, wire.inertia_after)
503    }
504}
505
506impl From<StabilizationLedger> for StabilizationLedgerWire {
507    fn from(ledger: StabilizationLedger) -> Self {
508        Self {
509            kind: ledger.kind,
510            delta: ledger.delta,
511            matrix_form: ledger.matrix_form,
512            chosen_by: ledger.chosen_by,
513            objective_policy: ledger.objective_policy,
514            backward_error_bound: ledger.backward_error_bound,
515            inertia_before: ledger.inertia_before,
516            inertia_after: ledger.inertia_after,
517        }
518    }
519}
520/// Generate a `#[repr(transparent)]` `Array1<f64>` newtype with the
521/// `new`/`Deref`/`DerefMut`/`AsRef`/`From` boilerplate used by unconstrained
522/// numeric vectors in this module.
523macro_rules! array1_f64_newtype {
524    ($name:ident) => {
525        #[repr(transparent)]
526        #[derive(Clone, Debug, PartialEq)]
527        pub struct $name(pub Array1<f64>);
528
529        impl $name {
530            #[inline]
531            pub fn new(values: Array1<f64>) -> Self {
532                Self(values)
533            }
534
535            #[inline]
536            pub fn zeros(len: usize) -> Self {
537                Self(Array1::zeros(len))
538            }
539        }
540
541        impl Deref for $name {
542            type Target = Array1<f64>;
543            #[inline]
544            fn deref(&self) -> &Self::Target {
545                &self.0
546            }
547        }
548
549        impl DerefMut for $name {
550            #[inline]
551            fn deref_mut(&mut self) -> &mut Self::Target {
552                &mut self.0
553            }
554        }
555
556        impl AsRef<Array1<f64>> for $name {
557            #[inline]
558            fn as_ref(&self) -> &Array1<f64> {
559                &self.0
560            }
561        }
562
563        impl From<Array1<f64>> for $name {
564            #[inline]
565            fn from(values: Array1<f64>) -> Self {
566                Self(values)
567            }
568        }
569
570        impl From<$name> for Array1<f64> {
571            #[inline]
572            fn from(values: $name) -> Self {
573                values.0
574            }
575        }
576    };
577}
578
579array1_f64_newtype!(Coefficients);
580array1_f64_newtype!(LinearPredictor);
581
582/// Index into `TermCollectionSpec::smooth_terms` (and the parallel
583/// `TermCollectionDesign::smooth.terms` slice produced from it).
584///
585/// This is **not** a penalty/ρ index, **not** a column index, and **not** a
586/// coefficient-offset index. Keeping it behind a `#[repr(transparent)]`
587/// newtype makes those confusables a compile error: a `SmoothTermIdx` cannot
588/// be silently used to index `rho`, `beta`, or a design column.
589#[repr(transparent)]
590#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord, Serialize, Deserialize)]
591pub struct SmoothTermIdx(usize);
592
593impl SmoothTermIdx {
594    #[inline]
595    pub const fn new(idx: usize) -> Self {
596        Self(idx)
597    }
598
599    /// Sentinel used by transient builders that must allocate a coord config
600    /// before the smooth term it references has been positioned in the spec.
601    /// Every code path that constructs a sentinel must overwrite it before
602    /// the value escapes the builder.
603    #[inline]
604    pub const fn placeholder() -> Self {
605        Self(usize::MAX)
606    }
607
608    #[inline]
609    pub const fn get(self) -> usize {
610        self.0
611    }
612
613}
614
615impl std::fmt::Display for SmoothTermIdx {
616    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
617        write!(f, "{}", self.0)
618    }
619}
620
621/// Index into the canonical penalty list `&[CanonicalPenalty]` — equivalently,
622/// the position of a smoothing parameter in the ρ / λ vector.
623///
624/// Penalty/ρ indices are not interchangeable with `SmoothTermIdx` (a smooth
625/// term can carry multiple canonical penalties — e.g. tensor-product double
626/// penalties — and structural penalties don't correspond to any smooth term).
627/// Keeping them as separate newtypes makes the historical bug pattern
628/// "indexed `rho` with a smooth-term ordinal" impossible to express.
629#[repr(transparent)]
630#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord, Serialize, Deserialize)]
631pub struct PenaltyIdx(usize);
632
633impl PenaltyIdx {
634    #[inline]
635    pub const fn new(idx: usize) -> Self {
636        Self(idx)
637    }
638
639    #[inline]
640    pub const fn get(self) -> usize {
641        self.0
642    }
643}
644
645impl std::fmt::Display for PenaltyIdx {
646    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
647        write!(f, "{}", self.0)
648    }
649}
650
651/// Index into a single smooth term's set of basis functions — i.e. the `k`
652/// in "the k-th basis function `B_k(x)` of this term".
653///
654/// Distinct from:
655///   * [`SmoothTermIdx`] — selects *which* smooth term in the spec.
656///   * [`PenaltyIdx`]    — selects *which* ρ/λ entry / canonical penalty.
657///   * A design-matrix column index — which lives in the *combined* layout
658///     after intercept/parametric blocks and per-term offsets are applied;
659///     a `BasisIdx` is term-local, a column index is model-global.
660///
661/// Keeping this as its own `#[repr(transparent)]` newtype makes the
662/// historically-easy confusion "indexed a global column slice with a
663/// term-local basis ordinal" (or vice versa) a compile error.
664#[repr(transparent)]
665#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord, Serialize, Deserialize)]
666pub struct BasisIdx(usize);
667
668impl BasisIdx {
669    #[inline]
670    pub const fn new(idx: usize) -> Self {
671        Self(idx)
672    }
673
674    #[inline]
675    pub const fn get(self) -> usize {
676        self.0
677    }
678}
679
680impl std::fmt::Display for BasisIdx {
681    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
682        write!(f, "{}", self.0)
683    }
684}
685
686/// Index into the user-facing design matrix `data: Array2<f64>` — i.e. the
687/// position of a covariate column in the raw input frame, *before* any
688/// per-family basis expansion or intercept/parametric layout is applied.
689///
690/// Distinct from:
691///   * [`BasisIdx`] — term-local basis-function ordinal `k` of `B_k(x)`.
692///   * [`SmoothTermIdx`] — position in `TermCollectionSpec::smooth_terms`.
693///   * A coefficient-vector offset `β[i]` — spans the combined design after
694///     expansion, which is much wider than the user-facing data matrix.
695///
696/// Keeping this as its own `#[repr(transparent)]` newtype rules out the easy
697/// confusion of indexing the raw data frame with an expanded-column offset.
698#[repr(transparent)]
699#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord, Serialize, Deserialize)]
700pub struct ColIdx(usize);
701
702impl ColIdx {
703    #[inline]
704    pub const fn new(idx: usize) -> Self {
705        Self(idx)
706    }
707
708    #[inline]
709    pub const fn get(self) -> usize {
710        self.0
711    }
712}
713
714impl std::fmt::Display for ColIdx {
715    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
716        write!(f, "{}", self.0)
717    }
718}
719
720/// Index of an observation (row) in the user-facing data frame / design
721/// matrix — i.e. the `i` in "the i-th observation".
722///
723/// Distinct from every column-type index in this module ([`ColIdx`],
724/// [`BasisIdx`], [`SmoothTermIdx`], [`PenaltyIdx`]) and from coefficient
725/// offsets. Keeping rows behind their own `#[repr(transparent)]` newtype
726/// makes the classic `data[[col, row]]` transposition a compile error.
727#[repr(transparent)]
728#[derive(Clone, Copy, Debug, PartialEq, Eq, Hash, PartialOrd, Ord, Serialize, Deserialize)]
729pub struct RowIdx(usize);
730
731impl RowIdx {
732    #[inline]
733    pub const fn new(idx: usize) -> Self {
734        Self(idx)
735    }
736
737    #[inline]
738    pub const fn get(self) -> usize {
739        self.0
740    }
741}
742
743impl std::fmt::Display for RowIdx {
744    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
745        write!(f, "{}", self.0)
746    }
747}
748
749#[repr(transparent)]
750#[derive(Clone, Copy, Debug)]
751pub struct LogSmoothingParamsView<'a>(ArrayView1<'a, f64>);
752
753impl<'a> LogSmoothingParamsView<'a> {
754    /// Borrow a smoothing vector only after every coordinate satisfies the
755    /// exact shared logarithmic-strength contract.
756    pub fn new(values: ArrayView1<'a, f64>) -> Result<Self, crate::IndexedLogStrengthDomainError> {
757        crate::validate_log_strengths(values.iter().copied())?;
758        Ok(Self(values))
759    }
760
761    /// Exact physical strengths for this already-validated vector.
762    pub fn exact_exp(&self) -> Array1<f64> {
763        // `new` established the private invariant; the borrow prevents the
764        // source array from being mutated for this view's lifetime.
765        self.0.mapv(f64::exp)
766    }
767}
768
769impl<'a> Deref for LogSmoothingParamsView<'a> {
770    type Target = ArrayView1<'a, f64>;
771
772    fn deref(&self) -> &Self::Target {
773        &self.0
774    }
775}
776
777#[cfg(test)]
778mod newtype_tests {
779    use super::*;
780    use ndarray::array;
781
782    #[test]
783    fn smooth_term_idx_ordering() {
784        let a = SmoothTermIdx::new(1);
785        let b = SmoothTermIdx::new(2);
786        assert!(a < b);
787        assert_eq!(a, SmoothTermIdx::new(1));
788    }
789
790    #[test]
791    fn coefficients_zeros_and_deref() {
792        let c = Coefficients::zeros(3);
793        assert_eq!(c.len(), 3);
794        assert!(c.iter().all(|&v| v == 0.0));
795    }
796
797    #[test]
798    fn coefficients_from_array1() {
799        let arr = array![1.0, 2.0, 3.0];
800        let c = Coefficients::from(arr.clone());
801        assert_eq!(*c, arr);
802    }
803
804    #[test]
805    fn log_smoothing_params_view_is_validated_and_exponentiates_exactly() {
806        let arr = array![crate::LOG_STRENGTH_MIN, 0.0, crate::LOG_STRENGTH_MAX];
807        let rho = LogSmoothingParamsView::new(arr.view()).expect("closed domain");
808        for (actual, expected) in rho.exact_exp().iter().zip(arr.iter()) {
809            assert_eq!(actual.to_bits(), expected.exp().to_bits());
810        }
811
812        let invalid = array![0.0, crate::LOG_STRENGTH_MAX + 1.0];
813        let error = LogSmoothingParamsView::new(invalid.view()).unwrap_err();
814        assert_eq!(error.coordinate, 1);
815        assert_eq!(error.value, crate::LOG_STRENGTH_MAX + 1.0);
816    }
817
818    #[test]
819    fn linear_predictor_zeros_and_deref() {
820        let lp = LinearPredictor::zeros(4);
821        assert_eq!(lp.len(), 4);
822        assert!(lp.iter().all(|&v| v == 0.0));
823    }
824}
825
826#[cfg(test)]
827mod ridge_policy_tests {
828    use super::{RidgePassport, RidgePolicy, StabilizationLedger};
829    use serde_json::json;
830
831    #[test]
832    fn serde_cannot_bypass_passport_validation() {
833        let negative = json!({
834            "delta": -1.0,
835            "matrix_form": "ScaledIdentity",
836            "policy": "SolverOnly"
837        });
838        assert!(serde_json::from_value::<RidgePassport>(negative).is_err());
839
840        let passport = RidgePassport::scaled_identity(
841            2.5e-7,
842            RidgePolicy::exact_full_objective(),
843        )
844        .expect("valid ridge");
845        let roundtrip: RidgePassport =
846            serde_json::from_value(serde_json::to_value(passport).expect("serialize passport"))
847                .expect("deserialize validated passport");
848        assert_eq!(roundtrip, passport);
849    }
850
851    #[test]
852    fn serde_cannot_bypass_ledger_semantics() {
853        let invalid_none = json!({
854            "kind": "None",
855            "delta": 1.0,
856            "matrix_form": "ScaledIdentity",
857            "chosen_by": "FixedConstant",
858            "objective_policy": null,
859            "backward_error_bound": null,
860            "inertia_before": null,
861            "inertia_after": null
862        });
863        assert!(serde_json::from_value::<StabilizationLedger>(invalid_none).is_err());
864
865        let invalid_prior_rule = json!({
866            "kind": "ExplicitPrior",
867            "delta": 1.0,
868            "matrix_form": "ScaledIdentity",
869            "chosen_by": "Heuristic",
870            "objective_policy": "ExactFullObjective",
871            "backward_error_bound": null,
872            "inertia_before": null,
873            "inertia_after": null
874        });
875        assert!(serde_json::from_value::<StabilizationLedger>(invalid_prior_rule).is_err());
876
877        let bound_on_wrong_kind = json!({
878            "kind": "ApproximationOnly",
879            "delta": 1.0,
880            "matrix_form": "ScaledIdentity",
881            "chosen_by": "FixedConstant",
882            "objective_policy": null,
883            "backward_error_bound": 1.0e-10,
884            "inertia_before": null,
885            "inertia_after": null
886        });
887        assert!(serde_json::from_value::<StabilizationLedger>(bound_on_wrong_kind).is_err());
888    }
889
890}