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/// Smallest supported logarithmic strength (inclusive).
14pub const LOG_STRENGTH_MIN: f64 = -700.0;
15
16/// Largest supported logarithmic strength (inclusive).
17pub const LOG_STRENGTH_MAX: f64 = 700.0;
18
19/// A logarithmic strength is outside the exact supported solver contract.
20#[derive(Debug, Clone, Copy, PartialEq)]
21pub struct LogStrengthDomainError {
22    pub value: f64,
23}
24
25impl std::fmt::Display for LogStrengthDomainError {
26    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
27        write!(
28            f,
29            "log strength must be finite and in [{LOG_STRENGTH_MIN}, {LOG_STRENGTH_MAX}]; got {}",
30            self.value
31        )
32    }
33}
34
35impl std::error::Error for LogStrengthDomainError {}
36
37#[derive(Debug, Clone, Copy, PartialEq)]
38pub struct PhysicalStrengthDomainError {
39    pub value: f64,
40}
41
42impl std::fmt::Display for PhysicalStrengthDomainError {
43    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
44        write!(
45            f,
46            "physical strength must be positive and finite with its logarithm in [{LOG_STRENGTH_MIN}, {LOG_STRENGTH_MAX}]; got {}",
47            self.value
48        )
49    }
50}
51
52impl std::error::Error for PhysicalStrengthDomainError {}
53
54/// Coordinate-aware failure returned when validating a vector of logarithmic
55/// strengths.
56#[derive(Debug, Clone, Copy, PartialEq)]
57pub struct IndexedLogStrengthDomainError {
58    pub coordinate: usize,
59    pub value: f64,
60}
61
62impl std::fmt::Display for IndexedLogStrengthDomainError {
63    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
64        write!(
65            f,
66            "log strength coordinate {} must be finite and in [{LOG_STRENGTH_MIN}, {LOG_STRENGTH_MAX}]; got {}",
67            self.coordinate, self.value
68        )
69    }
70}
71
72impl std::error::Error for IndexedLogStrengthDomainError {}
73
74impl From<IndexedLogStrengthDomainError> for crate::EstimationError {
75    fn from(error: IndexedLogStrengthDomainError) -> Self {
76        Self::LogStrengthDomainViolation {
77            coordinate: error.coordinate,
78            value: error.value,
79            lower: LOG_STRENGTH_MIN,
80            upper: LOG_STRENGTH_MAX,
81        }
82    }
83}
84
85/// Validate a logarithmic strength without changing it.
86#[inline]
87pub fn validate_log_strength(log_strength: f64) -> Result<(), LogStrengthDomainError> {
88    if log_strength.is_finite() && (LOG_STRENGTH_MIN..=LOG_STRENGTH_MAX).contains(&log_strength) {
89        Ok(())
90    } else {
91        Err(LogStrengthDomainError {
92            value: log_strength,
93        })
94    }
95}
96
97/// Validate a complete vector, reporting the deterministic smallest invalid
98/// coordinate before any caller-visible computation begins.
99pub fn validate_log_strengths(
100    values: impl IntoIterator<Item = f64>,
101) -> Result<(), IndexedLogStrengthDomainError> {
102    for (coordinate, value) in values.into_iter().enumerate() {
103        validate_log_strength(value)
104            .map_err(|_| IndexedLogStrengthDomainError { coordinate, value })?;
105    }
106    Ok(())
107}
108
109/// Convert a complete vector atomically on the exact supported domain.
110pub fn checked_exp_log_strengths(
111    values: impl IntoIterator<Item = f64>,
112) -> Result<Vec<f64>, IndexedLogStrengthDomainError> {
113    let mut strengths = Vec::new();
114    for (coordinate, value) in values.into_iter().enumerate() {
115        strengths.push(
116            checked_exp_log_strength(value)
117                .map_err(|_| IndexedLogStrengthDomainError { coordinate, value })?,
118        );
119    }
120    Ok(strengths)
121}
122
123/// Exponentiate a logarithmic strength on the exact closed solver domain.
124///
125/// No input is clamped and no output is floored or capped.  Thus the returned
126/// value is exactly the one whose first and second `rho` derivatives are both
127/// `exp(rho)`.
128#[inline]
129pub fn checked_exp_log_strength(log_strength: f64) -> Result<f64, LogStrengthDomainError> {
130    validate_log_strength(log_strength)?;
131    Ok(log_strength.exp())
132}
133
134/// Recover a canonical logarithmic coordinate without flooring or ceilinging
135/// a physical strength.
136pub fn checked_log_strength(strength: f64) -> Result<f64, PhysicalStrengthDomainError> {
137    if !(strength.is_finite() && strength > 0.0) {
138        return Err(PhysicalStrengthDomainError { value: strength });
139    }
140    let log_strength = strength.ln();
141    validate_log_strength(log_strength)
142        .map_err(|_| PhysicalStrengthDomainError { value: strength })?;
143    Ok(log_strength)
144}
145
146#[cfg(test)]
147mod tests {
148    use super::*;
149
150    #[test]
151    fn exact_closed_domain_accepts_both_endpoints_without_saturation() {
152        for endpoint in [LOG_STRENGTH_MIN, LOG_STRENGTH_MAX] {
153            let strength = checked_exp_log_strength(endpoint).expect("closed endpoint");
154            assert_eq!(strength.to_bits(), endpoint.exp().to_bits());
155            assert!(strength.is_finite() && strength > 0.0);
156        }
157    }
158
159    #[test]
160    fn exact_closed_domain_rejects_unsupported_and_nonfinite_values() {
161        for value in [
162            LOG_STRENGTH_MIN - 1.0,
163            LOG_STRENGTH_MAX + 1.0,
164            f64::NEG_INFINITY,
165            f64::INFINITY,
166            f64::NAN,
167        ] {
168            assert_eq!(
169                checked_exp_log_strength(value).unwrap_err().value.to_bits(),
170                value.to_bits()
171            );
172        }
173    }
174
175    #[test]
176    fn vector_validation_reports_the_smallest_bad_coordinate_atomically() {
177        let values = [0.0, LOG_STRENGTH_MAX + 1.0, f64::NAN];
178        let error = checked_exp_log_strengths(values).unwrap_err();
179        assert_eq!(error.coordinate, 1);
180        assert_eq!(error.value, LOG_STRENGTH_MAX + 1.0);
181    }
182
183    #[test]
184    fn physical_strength_conversion_refuses_floor_and_ceiling_cases() {
185        for value in [0.0, -1.0, f64::INFINITY, f64::NAN] {
186            assert!(checked_log_strength(value).is_err());
187        }
188        for endpoint in [LOG_STRENGTH_MIN, LOG_STRENGTH_MAX] {
189            assert!(checked_log_strength(endpoint.exp()).is_ok());
190        }
191    }
192}
gam_problem/log_strength.rs

gam_problem/
log_strength.rs
