spectrafit_models/
split_pearson7.rs

1use crate::Model;
2
3/// Split Pearson VII: a Pearson VII with a different width AND exponent on each side.
4///
5/// `f(x) = A / [1 + ((x − c)/σ_i)²·(2^{1/m_i} − 1)]^{m_i}`, with `(σ_i, m_i) = (σ_L, m_L)`
6/// for `x < c`, else `(σ_R, m_R)`.
7///
8/// Parameters (in order): `[amplitude, center, sigma_l, sigma_r, m_l, m_r]`. `amplitude` is
9/// the peak height at `x == center` (continuous: both branches give `A` there). numpy oracle
10/// identical via `np.where(x < c, left, right)`.
11pub struct SplitPearson7;
12
13fn p7(a: f64, z: f64, m: f64) -> f64 {
14    a / (1.0 + z * z * (2.0_f64.powf(1.0 / m) - 1.0)).powf(m)
15}
16
17impl Model for SplitPearson7 {
18    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
19        let (a, c, sl, sr, ml, mr) = (
20            params[0], params[1], params[2], params[3], params[4], params[5],
21        );
22        if x[0] < c {
23            p7(a, (x[0] - c) / sl, ml)
24        } else {
25            p7(a, (x[0] - c) / sr, mr)
26        }
27    }
28
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            "m_l".into(),
51            "m_r".into(),
52        ]
53    }
54}
55
56#[cfg(test)]
57mod tests {
58    use super::*;
59    use approx::assert_relative_eq;
60
61    #[test]
62    fn eval_at_center_equals_amplitude() {
63        let v = SplitPearson7.eval(&[1.5], &[3.0, 1.5, 0.6, 1.2, 2.0, 3.0]);
64        assert_relative_eq!(v, 3.0, epsilon = 1e-12);
65    }
66
67    #[test]
68    fn param_names_and_jacobian() {
69        assert_eq!(
70            SplitPearson7
71                .param_names()
72                .iter()
73                .map(|c| c.as_ref())
74                .collect::<Vec<_>>(),
75            &["amplitude", "center", "sigma_l", "sigma_r", "m_l", "m_r"]
76        );
77        let j = SplitPearson7.jacobian(&[1.0], &[3.0, 1.5, 0.6, 1.2, 2.0, 3.0]);
78        assert_eq!(j.len(), 6);
79    }
80
81    // ----- Limiting-case asymptotic (ground-truth verification) -----
82    //
83    // SplitPearson7 reduces to a plain (symmetric) Pearson VII when
84    // σ_l = σ_r and m_l = m_r:
85    //
86    //     A / [1 + ε²·(2^(1/m) − 1)]^m   uniformly for x < c and x ≥ c
87    //
88    // Verify against the plain Pearson7 kernel across symmetric x points.
89
90    #[test]
91    fn symmetric_widths_and_exponents_equal_plain_pearson7() {
92        use crate::pearson7::Pearson7;
93        let split = SplitPearson7;
94        let plain = Pearson7;
95        let sigma = 0.7_f64;
96        let m = 2.3_f64;
97        let a = 4.2_f64;
98        let c = 0.3_f64;
99        let p_split = [a, c, sigma, sigma, m, m]; // symmetric collapse
100        let p_plain = [a, c, sigma, m];
101        for &xi in &[-1.5_f64, -0.3, c, 0.5, 1.5] {
102            assert_relative_eq!(
103                split.eval(&[xi], &p_split),
104                plain.eval(&[xi], &p_plain),
105                epsilon = 1e-12
106            );
107        }
108    }
109}
spectrafit_models/split_pearson7.rs

spectrafit_models/
split_pearson7.rs
