spectrafit_models/
exponential.rs

1use crate::math_backend::batch_exp;
2use crate::Model;
3
4/// Double-exponential decay: `A1·exp(-λ1·x) + A2·exp(-λ2·x)`
5///
6/// Parameters (in order): `[A1, lam1, A2, lam2]`
7pub struct DoubleExponential;
8
9impl Model for DoubleExponential {
10    fn eval(&self, x: &[f64], p: &[f64]) -> f64 {
11        let xi = x[0];
12        p[0] * (-p[1] * xi).exp() + p[2] * (-p[3] * xi).exp()
13    }
14
15    fn jacobian_into(&self, x: &[f64], p: &[f64], out: &mut [f64]) {
16        let xi = x[0];
17        let e1 = (-p[1] * xi).exp();
18        let e2 = (-p[3] * xi).exp();
19        out[0] = e1; // ∂/∂A1
20        out[1] = -p[0] * xi * e1; // ∂/∂lam1
21        out[2] = e2; // ∂/∂A2
22        out[3] = -p[2] * xi * e2; // ∂/∂lam2
23    }
24
25    fn jacobian(&self, x: &[f64], p: &[f64]) -> Vec<f64> {
26        let xi = x[0];
27        let e1 = (-p[1] * xi).exp();
28        let e2 = (-p[3] * xi).exp();
29        vec![e1, -p[0] * xi * e1, e2, -p[2] * xi * e2]
30    }
31
32    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
33        vec!["A1".into(), "lam1".into(), "A2".into(), "lam2".into()]
34    }
35
36    fn eval_slice_into(&self, xs: &[f64], params: &[f64], out: &mut [f64]) {
37        debug_assert_eq!(out.len(), xs.len());
38        let (a1, lam1, a2, lam2) = (params[0], params[1], params[2], params[3]);
39        // Compute exp arguments for each component separately then sum.
40        let args1: Vec<f64> = xs.iter().map(|xi| -lam1 * xi).collect();
41        let args2: Vec<f64> = xs.iter().map(|xi| -lam2 * xi).collect();
42        let mut e1 = vec![0.0_f64; xs.len()];
43        let mut e2 = vec![0.0_f64; xs.len()];
44        batch_exp(&mut e1, &args1);
45        batch_exp(&mut e2, &args2);
46        for ((slot, v1), v2) in out.iter_mut().zip(e1.iter()).zip(e2.iter()) {
47            *slot = a1 * v1 + a2 * v2;
48        }
49    }
50
51    fn jac_slice_into(&self, xs: &[f64], params: &[f64], out: &mut [f64]) {
52        debug_assert_eq!(out.len(), xs.len() * 4);
53        let (a1, lam1, a2, lam2) = (params[0], params[1], params[2], params[3]);
54        let n = xs.len();
55        let args1: Vec<f64> = xs.iter().map(|xi| -lam1 * xi).collect();
56        let args2: Vec<f64> = xs.iter().map(|xi| -lam2 * xi).collect();
57        let mut e1 = vec![0.0_f64; n];
58        let mut e2 = vec![0.0_f64; n];
59        batch_exp(&mut e1, &args1);
60        batch_exp(&mut e2, &args2);
61        for (i, xi) in xs.iter().enumerate() {
62            out[i * 4] = e1[i]; // ∂/∂A1
63            out[i * 4 + 1] = -a1 * xi * e1[i]; // ∂/∂lam1
64            out[i * 4 + 2] = e2[i]; // ∂/∂A2
65            out[i * 4 + 3] = -a2 * xi * e2[i]; // ∂/∂lam2
66        }
67    }
68}
69
70#[cfg(test)]
71mod tests {
72    use super::*;
73    use approx::assert_relative_eq;
74
75    fn model() -> DoubleExponential {
76        DoubleExponential
77    }
78
79    // A1=2, lam1=0.5, A2=1, lam2=2 at x=0: 2*1 + 1*1 = 3
80    #[test]
81    fn eval_at_zero() {
82        let v = model().eval(&[0.0], &[2.0, 0.5, 1.0, 2.0]);
83        assert_relative_eq!(v, 3.0, epsilon = 1e-12);
84    }
85
86    #[test]
87    fn jacobian_shape() {
88        let j = model().jacobian(&[1.0], &[2.0, 0.5, 1.0, 2.0]);
89        assert_eq!(j.len(), 4);
90    }
91
92    // Central-difference (O(h^2)) sweep over three parameter regimes and
93    // several x-points, replacing the old single-point forward differences
94    // (h=1e-6, epsilon=1e-5). Each test keeps its name (traceability) but
95    // now checks only its own parameter index across the full grid.
96    const DOUBLE_EXP_REGIMES: [[f64; 4]; 3] = [
97        [2.0, 0.5, 1.0, 2.0],     // nominal
98        [1e-3, 1e-2, 1e-3, 1e-2], // small amplitudes/rates
99        [5.0, -1.0, -3.0, 4.0],   // negative rates (growth), large amplitudes
100    ];
101    const DOUBLE_EXP_XS: [f64; 5] = [-2.0, -0.5, 0.0, 1.0, 3.0];
102
103    fn check_double_exp_param_central_diff(idx: usize) {
104        for p in DOUBLE_EXP_REGIMES {
105            for &x in &DOUBLE_EXP_XS {
106                let j = model().jacobian(&[x], &p);
107                let h = 1e-6 * p[idx].abs().max(1.0);
108                let mut pp = p;
109                pp[idx] += h;
110                let mut pm = p;
111                pm[idx] -= h;
112                let fd = (model().eval(&[x], &pp) - model().eval(&[x], &pm)) / (2.0 * h);
113                assert_relative_eq!(j[idx], fd, epsilon = 1e-7, max_relative = 1e-6);
114            }
115        }
116    }
117
118    #[test]
119    fn jacobian_numerical_a1() {
120        check_double_exp_param_central_diff(0);
121    }
122
123    #[test]
124    fn jacobian_numerical_lam1() {
125        check_double_exp_param_central_diff(1);
126    }
127
128    #[test]
129    fn jacobian_numerical_a2() {
130        check_double_exp_param_central_diff(2);
131    }
132
133    #[test]
134    fn jacobian_numerical_lam2() {
135        check_double_exp_param_central_diff(3);
136    }
137
138    #[test]
139    fn eval_slice_matches_scalar() {
140        let xs: Vec<f64> = (0..10).map(|i| i as f64 * 0.5).collect();
141        let p = [2.0, 0.5, 1.0, 2.0];
142        let mut out = vec![0.0_f64; xs.len()];
143        model().eval_slice_into(&xs, &p, &mut out);
144        for (xi, &bi) in xs.iter().zip(out.iter()) {
145            assert_relative_eq!(bi, model().eval(&[*xi], &p), epsilon = 1e-10);
146        }
147    }
148
149    #[test]
150    fn jac_slice_matches_scalar() {
151        let xs: Vec<f64> = (0..8).map(|i| i as f64 * 0.5).collect();
152        let p = [2.0, 0.5, 1.0, 2.0];
153        let n = xs.len();
154        let mut out = vec![0.0_f64; n * 4];
155        model().jac_slice_into(&xs, &p, &mut out);
156        for (i, xi) in xs.iter().enumerate() {
157            let j = model().jacobian(&[*xi], &p);
158            for k in 0..4 {
159                assert_relative_eq!(out[i * 4 + k], j[k], epsilon = 1e-10);
160            }
161        }
162    }
163}
spectrafit_models/exponential.rs

spectrafit_models/
exponential.rs
