gam_problem/
log_strength.rs

1//! Exact supported domain for logarithmic penalty strengths.
2//!
3//! Every smoothing precision has the form `lambda = exp(rho)`.  The value,
4//! gradient, and Hessian with respect to `rho` agree only while that
5//! exponentiation is evaluated exactly: clamping `rho` or flooring/ceilinging
6//! `lambda` creates a constant tail with a fictitious nonzero derivative.
7//! The inclusive `[-700, 700]` interval deliberately stays inside binary64's
8//! finite, normal exponential range: its lower face avoids subnormal-strength
9//! arithmetic and its upper face leaves overflow guard margin.  It is a solver
10//! policy domain, not a claim about the widest representable binary64 input.
11//! This module owns the single domain used by all penalty implementations.
12
13/// `ln √ε`: the log of the relative resolution of a criterion gradient carried
14/// through an inverse whose conditioning is the strength ratio itself. Per
15/// direction the ρ-gradient is the effective degrees of freedom `γ/(γ+λ)`
16/// through `(H+λS)⁻¹`, whose condition in that direction is `λ/γ` once the
17/// penalty dominates; a quantity through an inverse of condition `κ` holds
18/// relative error `εκ`, and value and error cross at `λ/γ = 1/√ε`. Every
19/// derived ρ-domain edge (#2812) is this many e-folds from the spectrum.
20pub fn log_gradient_resolution() -> f64 {
21    0.5 * f64::EPSILON.ln()
22}
23
24/// The precision box `[ln √ε, ln(1/√ε)]` around unit strength: the domain of
25/// a coordinate whose penalty geometry cannot be projected, and the envelope a
26/// seed is placed in before the domain is derived.
27pub fn precision_box() -> (f64, f64) {
28    (log_gradient_resolution(), -log_gradient_resolution())
29}
30
31/// Smallest supported logarithmic strength (inclusive).
32pub const LOG_STRENGTH_MIN: f64 = -700.0;
33
34/// Largest supported logarithmic strength (inclusive).
35pub const LOG_STRENGTH_MAX: f64 = 700.0;
36
37/// A logarithmic strength is outside the exact supported solver contract.
38#[derive(Debug, Clone, Copy, PartialEq)]
39pub struct LogStrengthDomainError {
40    pub value: f64,
41}
42
43impl std::fmt::Display for LogStrengthDomainError {
44    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
45        write!(
46            f,
47            "log strength must be finite and in [{LOG_STRENGTH_MIN}, {LOG_STRENGTH_MAX}]; got {}",
48            self.value
49        )
50    }
51}
52
53impl std::error::Error for LogStrengthDomainError {}
54
55#[derive(Debug, Clone, Copy, PartialEq)]
56pub struct PhysicalStrengthDomainError {
57    pub value: f64,
58}
59
60impl std::fmt::Display for PhysicalStrengthDomainError {
61    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
62        write!(
63            f,
64            "physical strength must be positive and finite with its logarithm in [{LOG_STRENGTH_MIN}, {LOG_STRENGTH_MAX}]; got {}",
65            self.value
66        )
67    }
68}
69
70impl std::error::Error for PhysicalStrengthDomainError {}
71
72/// Coordinate-aware failure returned when validating a vector of logarithmic
73/// strengths.
74#[derive(Debug, Clone, Copy, PartialEq)]
75pub struct IndexedLogStrengthDomainError {
76    pub coordinate: usize,
77    pub value: f64,
78}
79
80impl std::fmt::Display for IndexedLogStrengthDomainError {
81    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
82        write!(
83            f,
84            "log strength coordinate {} must be finite and in [{LOG_STRENGTH_MIN}, {LOG_STRENGTH_MAX}]; got {}",
85            self.coordinate, self.value
86        )
87    }
88}
89
90impl std::error::Error for IndexedLogStrengthDomainError {}
91
92impl From<IndexedLogStrengthDomainError> for crate::EstimationError {
93    fn from(error: IndexedLogStrengthDomainError) -> Self {
94        Self::LogStrengthDomainViolation {
95            coordinate: error.coordinate,
96            value: error.value,
97            lower: LOG_STRENGTH_MIN,
98            upper: LOG_STRENGTH_MAX,
99        }
100    }
101}
102
103/// Validate a logarithmic strength without changing it.
104#[inline]
105pub fn validate_log_strength(log_strength: f64) -> Result<(), LogStrengthDomainError> {
106    if log_strength.is_finite() && (LOG_STRENGTH_MIN..=LOG_STRENGTH_MAX).contains(&log_strength) {
107        Ok(())
108    } else {
109        Err(LogStrengthDomainError {
110            value: log_strength,
111        })
112    }
113}
114
115/// Validate a complete vector, reporting the deterministic smallest invalid
116/// coordinate before any caller-visible computation begins.
117pub fn validate_log_strengths(
118    values: impl IntoIterator<Item = f64>,
119) -> Result<(), IndexedLogStrengthDomainError> {
120    for (coordinate, value) in values.into_iter().enumerate() {
121        validate_log_strength(value)
122            .map_err(|_| IndexedLogStrengthDomainError { coordinate, value })?;
123    }
124    Ok(())
125}
126
127/// Convert a complete vector atomically on the exact supported domain.
128pub fn checked_exp_log_strengths(
129    values: impl IntoIterator<Item = f64>,
130) -> Result<Vec<f64>, IndexedLogStrengthDomainError> {
131    let mut strengths = Vec::new();
132    for (coordinate, value) in values.into_iter().enumerate() {
133        strengths.push(
134            checked_exp_log_strength(value)
135                .map_err(|_| IndexedLogStrengthDomainError { coordinate, value })?,
136        );
137    }
138    Ok(strengths)
139}
140
141/// Exponentiate a logarithmic strength on the exact closed solver domain.
142///
143/// No input is clamped and no output is floored or capped.  Thus the returned
144/// value is exactly the one whose first and second `rho` derivatives are both
145/// `exp(rho)`.
146#[inline]
147pub fn checked_exp_log_strength(log_strength: f64) -> Result<f64, LogStrengthDomainError> {
148    validate_log_strength(log_strength)?;
149    Ok(log_strength.exp())
150}
151
152/// Recover a canonical logarithmic coordinate without flooring or ceilinging
153/// a physical strength.
154pub fn checked_log_strength(strength: f64) -> Result<f64, PhysicalStrengthDomainError> {
155    if !(strength.is_finite() && strength > 0.0) {
156        return Err(PhysicalStrengthDomainError { value: strength });
157    }
158    let log_strength = strength.ln();
159    validate_log_strength(log_strength)
160        .map_err(|_| PhysicalStrengthDomainError { value: strength })?;
161    Ok(log_strength)
162}
163
164#[cfg(test)]
165mod tests {
166    use super::*;
167
168    #[test]
169    fn exact_closed_domain_accepts_both_endpoints_without_saturation() {
170        for endpoint in [LOG_STRENGTH_MIN, LOG_STRENGTH_MAX] {
171            let strength = checked_exp_log_strength(endpoint).expect("closed endpoint");
172            assert_eq!(strength.to_bits(), endpoint.exp().to_bits());
173            assert!(strength.is_finite() && strength > 0.0);
174        }
175    }
176
177    #[test]
178    fn exact_closed_domain_rejects_unsupported_and_nonfinite_values() {
179        for value in [
180            LOG_STRENGTH_MIN - 1.0,
181            LOG_STRENGTH_MAX + 1.0,
182            f64::NEG_INFINITY,
183            f64::INFINITY,
184            f64::NAN,
185        ] {
186            assert_eq!(
187                checked_exp_log_strength(value).unwrap_err().value.to_bits(),
188                value.to_bits()
189            );
190        }
191    }
192
193    #[test]
194    fn vector_validation_reports_the_smallest_bad_coordinate_atomically() {
195        let values = [0.0, LOG_STRENGTH_MAX + 1.0, f64::NAN];
196        let error = checked_exp_log_strengths(values).unwrap_err();
197        assert_eq!(error.coordinate, 1);
198        assert_eq!(error.value, LOG_STRENGTH_MAX + 1.0);
199    }
200
201    #[test]
202    fn physical_strength_conversion_refuses_floor_and_ceiling_cases() {
203        for value in [0.0, -1.0, f64::INFINITY, f64::NAN] {
204            assert!(checked_log_strength(value).is_err());
205        }
206        for endpoint in [LOG_STRENGTH_MIN, LOG_STRENGTH_MAX] {
207            assert!(checked_log_strength(endpoint.exp()).is_ok());
208        }
209    }
210}
gam_problem/log_strength.rs

gam_problem/
log_strength.rs
