gam_terms/
isotropic_scale.rs

1//! A scalar coordinate scale for isotropic Euclidean smooths.
2//!
3//! Isotropic kernels admit one uniform change of coordinate units.  Encoding
4//! that scale as a vector made anisotropic states representable in the frozen
5//! model even though the kernel and its scale contract require one value in
6//! every direction.  `IsotropicScale` makes the geometric invariant explicit:
7//! anisotropy belongs to the separate ARD parameters, never to this frame.
8
9use ndarray::Array2;
10use serde::{Deserialize, Serialize};
11use std::fmt;
12
13/// A coordinate-valued length in the user's ORIGINAL covariate units.
14///
15/// A basis that auto-standardizes its Euclidean input divides the coordinates
16/// by an [`IsotropicScale`] before building anything, so its `centers` — and
17/// every radius the kernel evaluates — live in the standardized frame while
18/// the range the user asked for lives here.  Keeping the two frames in
19/// distinct types is what stops a consumer from pairing them: the kernel
20/// range and the radii it is compared against must agree, and getting that
21/// wrong is a silent O(1) relative error in the kernel bandwidth that no
22/// shape or count check can see.
23#[repr(transparent)]
24#[derive(Clone, Copy, Debug, PartialEq, PartialOrd, Serialize, Deserialize)]
25#[serde(transparent)]
26pub struct OriginalUnits(f64);
27
28/// A coordinate-valued length in the STANDARDIZED frame, i.e. the frame the
29/// auto-standardized `centers` and the kernel radii already live in.
30///
31/// This is the frame every radial kernel evaluation must be denominated in.
32#[repr(transparent)]
33#[derive(Clone, Copy, Debug, PartialEq, PartialOrd, Serialize, Deserialize)]
34#[serde(transparent)]
35pub struct StandardizedUnits(f64);
36
37impl OriginalUnits {
38    pub const fn new(value: f64) -> Self {
39        Self(value)
40    }
41
42    /// The scalar, named for its frame so a bare extraction that feeds
43    /// standardized-frame math reads wrong at the call site.
44    pub const fn original_value(self) -> f64 {
45        self.0
46    }
47}
48
49impl StandardizedUnits {
50    pub const fn new(value: f64) -> Self {
51        Self(value)
52    }
53
54    /// The scalar, named for its frame; see [`OriginalUnits::original_value`].
55    pub const fn standardized_value(self) -> f64 {
56        self.0
57    }
58}
59
60impl fmt::Display for OriginalUnits {
61    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
62        fmt::Display::fmt(&self.0, formatter)
63    }
64}
65
66impl fmt::Display for StandardizedUnits {
67    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
68        fmt::Display::fmt(&self.0, formatter)
69    }
70}
71
72/// A positive, finite scalar whose reciprocal is also representable.
73///
74/// The field is private so construction, deserialization, and every frozen
75/// replay path enforce the same invariant.
76#[repr(transparent)]
77#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
78#[serde(try_from = "f64", into = "f64")]
79pub struct IsotropicScale(f64);
80
81impl IsotropicScale {
82    pub const ONE: Self = Self(1.0);
83
84    pub fn new(value: f64) -> Result<Self, IsotropicScaleError> {
85        if value.is_finite() && value > 0.0 && value.recip().is_finite() {
86            Ok(Self(value))
87        } else {
88            Err(IsotropicScaleError { value })
89        }
90    }
91
92    pub fn get(self) -> f64 {
93        self.0
94    }
95
96    pub fn reciprocal(self) -> f64 {
97        self.0.recip()
98    }
99
100    pub fn to_bits(self) -> u64 {
101        self.0.to_bits()
102    }
103
104    /// Convert a coordinate-valued length from original to standardized units.
105    ///
106    /// This is the ONLY conversion between the two frames.  Because it
107    /// consumes an [`OriginalUnits`] and produces a [`StandardizedUnits`],
108    /// applying it to a value that is already standardized — the
109    /// double-divide that silently halves or doubles a kernel bandwidth — is
110    /// a type error rather than a comment for the next reader to remember.
111    pub fn to_standardized_units(self, value: OriginalUnits) -> StandardizedUnits {
112        StandardizedUnits(value.original_value() * self.reciprocal())
113    }
114
115    /// The inverse of [`Self::to_standardized_units`]: express a standardized
116    /// length back in the user's original coordinate units.
117    pub fn to_original_units(self, value: StandardizedUnits) -> OriginalUnits {
118        OriginalUnits(value.standardized_value() * self.0)
119    }
120
121    /// Apply the uniform coordinate pullback in place.
122    pub fn standardize(self, coordinates: &mut Array2<f64>) {
123        let reciprocal = self.reciprocal();
124        coordinates.mapv_inplace(|value| value * reciprocal);
125    }
126}
127
128impl TryFrom<f64> for IsotropicScale {
129    type Error = IsotropicScaleError;
130
131    fn try_from(value: f64) -> Result<Self, Self::Error> {
132        Self::new(value)
133    }
134}
135
136impl From<IsotropicScale> for f64 {
137    fn from(value: IsotropicScale) -> Self {
138        value.get()
139    }
140}
141
142#[derive(Clone, Copy, Debug, PartialEq)]
143pub struct IsotropicScaleError {
144    value: f64,
145}
146
147impl fmt::Display for IsotropicScaleError {
148    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
149        write!(
150            formatter,
151            "isotropic scale must be positive and finite with a finite reciprocal, got {}",
152            self.value
153        )
154    }
155}
156
157impl std::error::Error for IsotropicScaleError {}
158
159#[cfg(test)]
160mod tests {
161    use super::*;
162
163    #[test]
164    fn construction_enforces_the_operational_invariant() {
165        assert_eq!(IsotropicScale::new(1.0), Ok(IsotropicScale::ONE));
166        assert!(IsotropicScale::new(f64::MIN_POSITIVE).is_ok());
167        assert!(IsotropicScale::new(0.0).is_err());
168        assert!(IsotropicScale::new(-1.0).is_err());
169        assert!(IsotropicScale::new(f64::NAN).is_err());
170        assert!(IsotropicScale::new(f64::INFINITY).is_err());
171        assert!(IsotropicScale::new(f64::from_bits(1)).is_err());
172    }
173
174    #[test]
175    fn the_two_frames_round_trip_through_the_only_conversion() {
176        let scale = IsotropicScale::new(4.0).unwrap();
177        let original = OriginalUnits::new(500.0);
178        let standardized = scale.to_standardized_units(original);
179        assert_eq!(standardized, StandardizedUnits::new(125.0));
180        assert_eq!(scale.to_original_units(standardized), original);
181        // The frame tag is what distinguishes the two, not the magnitude:
182        // a bare `125.0` carries no evidence of which side it came from,
183        // which is precisely the hazard the tags remove.
184        assert_eq!(
185            scale.to_standardized_units(OriginalUnits::new(125.0)),
186            StandardizedUnits::new(31.25)
187        );
188    }
189
190    #[test]
191    fn wire_representation_is_a_checked_scalar() {
192        let encoded = serde_json::to_string(&IsotropicScale::new(2.5).unwrap()).unwrap();
193        assert_eq!(encoded, "2.5");
194        assert_eq!(
195            serde_json::from_str::<IsotropicScale>(&encoded).unwrap(),
196            IsotropicScale::new(2.5).unwrap()
197        );
198        assert!(serde_json::from_str::<IsotropicScale>("[2.5,2.5]").is_err());
199        assert!(serde_json::from_str::<IsotropicScale>("0.0").is_err());
200        assert!(serde_json::from_str::<IsotropicScale>("-1.0").is_err());
201    }
202}
gam_terms/isotropic_scale.rs

gam_terms/
isotropic_scale.rs
