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#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
11pub enum RidgeMatrixForm {
12 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#[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 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#[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#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
239pub enum StabilizationRule {
240 FixedConstant,
242 InertiaTarget { spd_floor: f64 },
244 Heuristic,
246 UserSpecified,
248 BackoffEscalation { attempts: usize },
250}
251
252#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
254pub enum StabilizationKind {
255 None,
256 SolverDampingOnly,
259 NumericalPerturbation,
264 ApproximationOnly,
268 ObjectiveStabilization,
271 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#[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 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 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}
520macro_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#[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 #[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#[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#[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#[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#[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 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 pub fn exact_exp(&self) -> Array1<f64> {
763 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}