spectrafit_models/
lorentzian.rs

1use crate::Model;
2
3/// Lorentzian kernel: `A / (1 + ((x₀ - c) / σ)²)`
4///
5/// Parameters (in order): `[amplitude, center, sigma]`
6pub struct Lorentzian;
7
8impl Model for Lorentzian {
9    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
10        let (a, c, sigma) = (params[0], params[1], params[2]);
11        let d = (x[0] - c) / sigma;
12        a / (1.0 + d * d)
13    }
14
15    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
16        let (a, c, sigma) = (params[0], params[1], params[2]);
17        let dx = x[0] - c;
18        let d = dx / sigma;
19        let big_d = 1.0 + d * d;
20
21        // ∂/∂amplitude = 1/D
22        let da = 1.0 / big_d;
23        // ∂/∂center    = 2*A*(x₀-c) / (σ²*D²)
24        let dc = 2.0 * a * dx / (sigma * sigma * big_d * big_d);
25        // ∂/∂sigma     = 2*A*(x₀-c)² / (σ³*D²)
26        let ds = 2.0 * a * dx * dx / (sigma * sigma * sigma * big_d * big_d);
27
28        vec![da, dc, ds]
29    }
30
31    #[inline]
32    fn jacobian_into(&self, x: &[f64], params: &[f64], out: &mut [f64]) {
33        let (a, c, sigma) = (params[0], params[1], params[2]);
34        let dx = x[0] - c;
35        let d = dx / sigma;
36        let big_d = 1.0 + d * d;
37        let big_d2 = big_d * big_d;
38        let s2 = sigma * sigma;
39        out[0] = 1.0 / big_d;
40        out[1] = 2.0 * a * dx / (s2 * big_d2);
41        out[2] = 2.0 * a * dx * dx / (s2 * sigma * big_d2);
42    }
43
44    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
45        vec!["amplitude".into(), "center".into(), "sigma".into()]
46    }
47}
48
49#[cfg(test)]
50mod tests {
51    use super::*;
52    use approx::assert_relative_eq;
53
54    #[test]
55    fn eval_at_center() {
56        // At x=center the Lorentzian equals amplitude.
57        let l = Lorentzian;
58        let v = l.eval(&[0.0], &[3.0, 0.0, 1.0]);
59        assert_relative_eq!(v, 3.0, epsilon = 1e-12);
60    }
61
62    #[test]
63    fn eval_at_sigma_offset() {
64        // At x = center + σ: A / (1 + 1) = A/2.
65        let l = Lorentzian;
66        let v = l.eval(&[1.0], &[4.0, 0.0, 1.0]);
67        assert_relative_eq!(v, 2.0, epsilon = 1e-12);
68    }
69
70    #[test]
71    fn jacobian_shape() {
72        let l = Lorentzian;
73        let j = l.jacobian(&[0.5], &[1.0, 0.0, 1.0]);
74        assert_eq!(j.len(), 3);
75    }
76
77    #[test]
78    fn jacobian_at_center_dc_zero() {
79        // ∂/∂center at x == center is 0.
80        let l = Lorentzian;
81        let j = l.jacobian(&[0.0], &[1.0, 0.0, 1.0]);
82        assert_relative_eq!(j[1], 0.0, epsilon = 1e-12);
83    }
84
85    // Central-difference (O(h^2)) sweep over three parameter regimes and
86    // several x-points relative to the centre, replacing the old
87    // single-point forward differences (h=1e-6, epsilon=1e-5). Each test
88    // keeps its name (traceability) but now checks only its own parameter
89    // index across the full regime/x-point grid.
90    const LORENTZIAN_REGIMES: [[f64; 3]; 3] = [
91        [2.0, 0.0, 1.5],   // nominal
92        [1e-3, 0.0, 0.05], // small amplitude/width
93        [5.0, -2.0, 3.0],  // offset centre, wide
94    ];
95
96    fn check_lorentzian_param_central_diff(idx: usize) {
97        let l = Lorentzian;
98        for p in LORENTZIAN_REGIMES {
99            let (c, sigma) = (p[1], p[2]);
100            for &mult in &[-3.0_f64, -1.0, -0.1, 0.0, 0.1, 1.0, 3.0] {
101                let x = [c + mult * sigma];
102                let j = l.jacobian(&x, &p);
103                let h = 1e-6 * p[idx].abs().max(1.0);
104                let (mut a, mut b) = (p, p);
105                a[idx] += h;
106                b[idx] -= h;
107                let fd = (l.eval(&x, &a) - l.eval(&x, &b)) / (2.0 * h);
108                assert_relative_eq!(j[idx], fd, epsilon = 1e-7, max_relative = 1e-6);
109            }
110        }
111    }
112
113    #[test]
114    fn jacobian_numerical_check_amplitude() {
115        check_lorentzian_param_central_diff(0);
116    }
117
118    #[test]
119    fn jacobian_numerical_check_center() {
120        check_lorentzian_param_central_diff(1);
121    }
122
123    #[test]
124    fn jacobian_numerical_check_sigma() {
125        check_lorentzian_param_central_diff(2);
126    }
127}
spectrafit_models/lorentzian.rs

spectrafit_models/
lorentzian.rs
