spectrafit_models/
split_gaussian.rs

1use crate::Model;
2
3/// Split (asymmetric) Gaussian: a Gaussian with a different width on each side of
4/// the center. Covers both "asymmetric split-σ Gaussian" and "bi-Gaussian".
5///
6/// `f(x) = A · exp(−(x−c)² / (2σ²))` with `σ = σ_L` for `x < c`, else `σ_R`.
7///
8/// Parameters (in order): `[amplitude, center, sigma_l, sigma_r]`
9///
10/// - `amplitude` is the peak height at `x == center` (both branches equal `A` there,
11///   so the curve is continuous).
12/// - `center` is the peak location.
13/// - `sigma_l` / `sigma_r` are the left/right Gaussian widths.
14///
15/// The numpy oracle is identical —
16/// `np.where(x < c, A·exp(−½((x−c)/σ_L)²), A·exp(−½((x−c)/σ_R)²))` — so parity is exact.
17pub struct SplitGaussian;
18
19impl Model for SplitGaussian {
20    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
21        let (a, c, sl, sr) = (params[0], params[1], params[2], params[3]);
22        let dx = x[0] - c;
23        let sigma = if x[0] < c { sl } else { sr };
24        a * (-(dx * dx) / (2.0 * sigma * sigma)).exp()
25    }
26
27    /// Central finite-difference Jacobian (the side-selecting branch makes the
28    /// analytic derivative piecewise; numerical is simpler and matches `log_normal`).
29    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
30        let mut p = params.to_vec();
31        (0..params.len())
32            .map(|i| {
33                let h = 1e-7_f64 * params[i].abs().max(1e-7);
34                p[i] = params[i] + h;
35                let f_plus = self.eval(x, &p);
36                p[i] = params[i] - h;
37                let f_minus = self.eval(x, &p);
38                p[i] = params[i];
39                (f_plus - f_minus) / (2.0 * h)
40            })
41            .collect()
42    }
43
44    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
45        vec![
46            "amplitude".into(),
47            "center".into(),
48            "sigma_l".into(),
49            "sigma_r".into(),
50        ]
51    }
52}
53
54#[cfg(test)]
55mod tests {
56    use super::*;
57    use approx::assert_relative_eq;
58
59    #[test]
60    fn eval_at_center_equals_amplitude() {
61        let m = SplitGaussian;
62        assert_relative_eq!(m.eval(&[1.5], &[3.0, 1.5, 0.6, 1.2]), 3.0, epsilon = 1e-12);
63    }
64
65    #[test]
66    fn left_and_right_use_their_own_width() {
67        // At ±σ from center the value is A·exp(−0.5) on each side, using its own σ.
68        let m = SplitGaussian;
69        let p = [4.0, 0.0, 0.6, 1.2];
70        assert_relative_eq!(m.eval(&[-0.6], &p), 4.0 * (-0.5_f64).exp(), epsilon = 1e-12);
71        assert_relative_eq!(m.eval(&[1.2], &p), 4.0 * (-0.5_f64).exp(), epsilon = 1e-12);
72    }
73
74    #[test]
75    fn param_names_are_canonical() {
76        assert_eq!(
77            SplitGaussian
78                .param_names()
79                .iter()
80                .map(|c| c.as_ref())
81                .collect::<Vec<_>>(),
82            &["amplitude", "center", "sigma_l", "sigma_r"]
83        );
84    }
85
86    #[test]
87    fn jacobian_shape() {
88        let j = SplitGaussian.jacobian(&[0.5], &[3.0, 0.0, 0.6, 1.2]);
89        assert_eq!(j.len(), 4);
90        // x = 0.5 > center ⇒ right width 1.2; ∂/∂A = exp(−0.5²/(2·1.2²)).
91        assert_relative_eq!(j[0], (-0.25_f64 / (2.0 * 1.44)).exp(), epsilon = 1e-6);
92    }
93
94    // ----- Limiting-case asymptotic (ground-truth verification) -----
95    //
96    // SplitGaussian reduces to a plain symmetric Gaussian when σ_l == σ_r:
97    //
98    //     A · exp(−(x−c)²/(2σ²))   with σ_l = σ_r = σ
99    //
100    // Both branches give the same value at any x, so the curve is identical to
101    // a plain Gaussian everywhere — pin this across left/right/center points.
102
103    #[test]
104    fn symmetric_widths_equal_plain_gaussian() {
105        use crate::gaussian::Gaussian;
106        let split = SplitGaussian;
107        let gauss = Gaussian;
108        let sigma = 0.8_f64;
109        let a = 3.5_f64;
110        let c = 0.4_f64;
111        let p_split = [a, c, sigma, sigma]; // σ_l = σ_r
112        let p_gauss = [a, c, sigma];
113        for &xi in &[-2.0_f64, -0.6, c, 0.6, 2.0] {
114            assert_relative_eq!(
115                split.eval(&[xi], &p_split),
116                gauss.eval(&[xi], &p_gauss),
117                epsilon = 1e-12
118            );
119        }
120    }
121}
spectrafit_models/split_gaussian.rs

spectrafit_models/
split_gaussian.rs
