spectrafit_models/
pearson7.rs

1use crate::Model;
2
3/// Pearson VII peak: `A / [1 + ((x − c)/σ)² · (2^{1/m} − 1)]^m`.
4///
5/// Parameters (in order): `[amplitude, center, sigma, m]`
6///
7/// - `amplitude` is the peak height attained at `x == center`.
8/// - `center` is the peak location.
9/// - `sigma` is the half-width at half-maximum (the `2^{1/m}−1` factor normalizes
10///   so `σ` is the HWHM for any shape exponent).
11/// - `m` is the shape exponent: `m → 1` gives a Lorentzian, `m → ∞` a Gaussian.
12///
13/// The numpy benchmark oracle is identical —
14/// `A / (1 + ((x−c)/σ)² · (2**(1/m) − 1))**m` — so numpy↔Rust parity is exact.
15pub struct Pearson7;
16
17impl Model for Pearson7 {
18    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
19        let (a, c, sigma, m) = (params[0], params[1], params[2], params[3]);
20        let z = (x[0] - c) / sigma;
21        let base = 1.0 + z * z * (2.0_f64.powf(1.0 / m) - 1.0);
22        a / base.powf(m)
23    }
24
25    /// Central finite-difference Jacobian (the `2^{1/m}` term makes the analytic
26    /// `m`-derivative awkward, so a numerical Jacobian is used, matching `log_normal`).
27    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
28        let mut p = params.to_vec();
29        (0..params.len())
30            .map(|i| {
31                let h = 1e-7_f64 * params[i].abs().max(1e-7);
32                p[i] = params[i] + h;
33                let f_plus = self.eval(x, &p);
34                p[i] = params[i] - h;
35                let f_minus = self.eval(x, &p);
36                p[i] = params[i];
37                (f_plus - f_minus) / (2.0 * h)
38            })
39            .collect()
40    }
41
42    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
43        vec![
44            "amplitude".into(),
45            "center".into(),
46            "sigma".into(),
47            "m".into(),
48        ]
49    }
50}
51
52#[cfg(test)]
53mod tests {
54    use super::*;
55    use approx::assert_relative_eq;
56
57    #[test]
58    fn eval_at_center_equals_amplitude() {
59        // At x == center, z = 0 ⇒ base = 1, so value == amplitude for any m.
60        let m = Pearson7;
61        assert_relative_eq!(m.eval(&[1.5], &[3.0, 1.5, 0.8, 2.0]), 3.0, epsilon = 1e-12);
62    }
63
64    #[test]
65    fn sigma_is_hwhm() {
66        // At |x − c| == sigma the normalization gives exactly half-maximum.
67        let m = Pearson7;
68        let v = m.eval(&[2.3], &[4.0, 1.5, 0.8, 2.5]); // x = c + sigma
69        assert_relative_eq!(v, 2.0, epsilon = 1e-12); // A/2
70    }
71
72    #[test]
73    fn param_names_are_canonical() {
74        assert_eq!(
75            Pearson7
76                .param_names()
77                .iter()
78                .map(|c| c.as_ref())
79                .collect::<Vec<_>>(),
80            &["amplitude", "center", "sigma", "m"]
81        );
82    }
83
84    #[test]
85    fn jacobian_shape_and_amplitude() {
86        let j = Pearson7.jacobian(&[1.5], &[3.0, 1.5, 0.8, 2.0]);
87        assert_eq!(j.len(), 4);
88        // ∂/∂amplitude at x == center is 1/base = 1.
89        assert_relative_eq!(j[0], 1.0, epsilon = 1e-6);
90        // ∂/∂center at the peak is 0 (stationary).
91        assert_relative_eq!(j[1], 0.0, epsilon = 1e-6);
92    }
93}
spectrafit_models/pearson7.rs

spectrafit_models/
pearson7.rs
