spectrafit_models/
exp_over_linear.rs

1use crate::Model;
2
3/// Exponential decay divided by a line.
4///
5/// ```text
6/// y = exp(−rate·x) / (lin_const + lin_slope·x)
7/// ```
8///
9/// Parameters (in order): `[rate, lin_const, lin_slope]`
10///
11/// This is the NIST StRD Chwirut model, shared by two datasets: Chwirut1 (214
12/// observations) and Chwirut2 (54), which differ only in the data.
13///
14/// **Domain guard:** the denominator must be non-zero. Unlike a quadratic
15/// denominator this one has a single root at `x = −lin_const/lin_slope`, which a
16/// solver can walk onto during search even when the certified parameters keep it
17/// clear of the data, so `D = 0` returns `f64::NAN`.
18///
19/// **Analytic Jacobian** (let `E = exp(−rate·x)`, `D = lin_const + lin_slope·x`):
20/// - ∂y/∂rate      = −x·E / D
21/// - ∂y/∂lin_const = −E / D²
22/// - ∂y/∂lin_slope   = −x·E / D²
23pub struct ExpOverLinear;
24
25#[inline]
26fn ed(xi: f64, p: &[f64]) -> (f64, f64) {
27    ((-p[0] * xi).exp(), p[1] + p[2] * xi)
28}
29
30#[inline]
31fn jac_at(xi: f64, p: &[f64], out: &mut [f64]) {
32    let (e, d) = ed(xi, p);
33    if d == 0.0 {
34        out[..3].fill(f64::NAN);
35        return;
36    }
37    let inv_d = 1.0 / d;
38    let e_over_d2 = e * inv_d * inv_d;
39    out[0] = -xi * e * inv_d;
40    out[1] = -e_over_d2;
41    out[2] = -xi * e_over_d2;
42}
43
44impl Model for ExpOverLinear {
45    fn eval(&self, x: &[f64], p: &[f64]) -> f64 {
46        let (e, d) = ed(x[0], p);
47        if d == 0.0 {
48            return f64::NAN;
49        }
50        e / d
51    }
52
53    fn jacobian_into(&self, x: &[f64], p: &[f64], out: &mut [f64]) {
54        jac_at(x[0], p, out);
55    }
56
57    fn jacobian(&self, x: &[f64], p: &[f64]) -> Vec<f64> {
58        let mut out = vec![0.0_f64; 3];
59        jac_at(x[0], p, &mut out);
60        out
61    }
62
63    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
64        vec!["rate".into(), "lin_const".into(), "lin_slope".into()]
65    }
66
67    fn eval_slice_into(&self, xs: &[f64], params: &[f64], out: &mut [f64]) {
68        debug_assert_eq!(out.len(), xs.len());
69        for (slot, &xi) in out.iter_mut().zip(xs.iter()) {
70            let (e, d) = ed(xi, params);
71            *slot = if d == 0.0 { f64::NAN } else { e / d };
72        }
73    }
74
75    fn jac_slice_into(&self, xs: &[f64], params: &[f64], out: &mut [f64]) {
76        debug_assert_eq!(out.len(), xs.len() * 3);
77        for (i, &xi) in xs.iter().enumerate() {
78            jac_at(xi, params, &mut out[i * 3..i * 3 + 3]);
79        }
80    }
81}
82
83#[cfg(test)]
84mod tests {
85    use super::*;
86    use approx::assert_relative_eq;
87
88    fn model() -> ExpOverLinear {
89        ExpOverLinear
90    }
91
92    // Chwirut2's certified values.
93    const CHWIRUT: [f64; 3] = [1.6657666537e-01, 5.1653291286e-03, 1.2150007096e-02];
94
95    // Anchored on the certified fit, not the observations: at these parameters the
96    // model reproduces Chwirut2's certified residual sum of squares (5.130480e2)
97    // exactly. The data scatters widely about it (92.9 observed against 81.86
98    // fitted at x=0.5), so an observation-anchored test measures the scatter.
99    #[test]
100    fn eval_chwirut_reproduces_the_certified_fit() {
101        assert_relative_eq!(
102            model().eval(&[0.5], &CHWIRUT),
103            81.855_746_144_55,
104            max_relative = 1e-10
105        );
106        assert_relative_eq!(model().eval(&[6.0], &CHWIRUT), 4.715_0, max_relative = 1e-3);
107    }
108
109    #[test]
110    fn zero_denominator_returns_nan() {
111        // D = lin_const + lin_slope·x = 0 at x = 1 when lin_const = -lin_slope.
112        let p = [0.1_f64, -1.0, 1.0];
113        assert!(model().eval(&[1.0], &p).is_nan());
114        assert!(model().jacobian(&[1.0], &p).iter().all(|v| v.is_nan()));
115    }
116
117    #[test]
118    fn jacobian_matches_finite_difference() {
119        // Already a central difference; widened from one regime/x-point to
120        // three regimes and several x-points (avoiding D = lin_const +
121        // lin_slope·x = 0, the domain-guard boundary), and tightened epsilon
122        // from 1e-6 toward 1e-7.
123        let param_sets = [
124            CHWIRUT,                // certified Chwirut2 fit (nominal)
125            [1e-3_f64, 1e-2, 1e-3], // small rate/slope
126            [5.0_f64, -2.0, 8.0],   // large rate, offset/steep line
127        ];
128        for p in param_sets {
129            for &x in &[0.1_f64, 0.5, 2.5, 6.0, 12.0] {
130                let j = model().jacobian(&[x], &p);
131                for k in 0..3 {
132                    let h = 1e-6 * p[k].abs().max(1.0);
133                    let mut pp = p;
134                    pp[k] += h;
135                    let mut pm = p;
136                    pm[k] -= h;
137                    let fd = (model().eval(&[x], &pp) - model().eval(&[x], &pm)) / (2.0 * h);
138                    assert_relative_eq!(j[k], fd, max_relative = 1e-6, epsilon = 1e-7);
139                }
140            }
141        }
142    }
143
144    #[test]
145    fn slices_match_scalar() {
146        let xs = [0.5_f64, 1.0, 3.0, 6.0];
147        let mut ys = vec![0.0_f64; xs.len()];
148        model().eval_slice_into(&xs, &CHWIRUT, &mut ys);
149        let mut js = vec![0.0_f64; xs.len() * 3];
150        model().jac_slice_into(&xs, &CHWIRUT, &mut js);
151        for (i, &xi) in xs.iter().enumerate() {
152            assert_relative_eq!(ys[i], model().eval(&[xi], &CHWIRUT), epsilon = 1e-12);
153            let j = model().jacobian(&[xi], &CHWIRUT);
154            for k in 0..3 {
155                assert_relative_eq!(js[i * 3 + k], j[k], epsilon = 1e-12);
156            }
157        }
158    }
159}
spectrafit_models/exp_over_linear.rs

spectrafit_models/
exp_over_linear.rs
