spectrafit_models/
skewed_gaussian.rs

1use crate::Model;
2
3/// Skewed Gaussian: a Gaussian modulated by an error-function skew factor.
4///
5/// `A · exp(−½((x−c)/σ)²) · (1 + erf(γ·(x−c)/(σ·√2)))`
6///
7/// Parameters (in order): `[amplitude, center, sigma, gamma]`
8///
9/// - `gamma` is the skewness: `0` ⇒ symmetric Gaussian, `γ>0` ⇒ tail toward high
10///   `x`, `γ<0` ⇒ tail toward low `x`. Uses `libm::erf`.
11pub struct SkewedGaussian;
12
13impl Model for SkewedGaussian {
14    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
15        let (a, c, sigma, gamma) = (params[0], params[1], params[2], params[3]);
16        let dx = x[0] - c;
17        let g = (-0.5 * (dx / sigma) * (dx / sigma)).exp();
18        let beta = gamma / (sigma * std::f64::consts::SQRT_2);
19        a * g * (1.0 + libm::erf(beta * dx))
20    }
21
22    /// Analytical Jacobian of the skewed Gaussian.
23    ///
24    /// Let `g = exp(−½(dx/σ)²)`, `β = γ/(σ√2)`, `skew = 1 + erf(β·dx)`,
25    /// `erf_d = (2/√π)·exp(−(β·dx)²)` (the erf derivative factor).
26    ///
27    /// ∂f/∂A      = g · skew
28    /// ∂f/∂center = A·g · [  dx/σ² · skew − β · erf_d ]   (∂dx/∂c = −1)
29    /// ∂f/∂sigma  = A·g · [ dx²/σ³ · skew − β·dx/σ · erf_d ]
30    /// ∂f/∂gamma  = A·g · [ dx/(σ√2) · erf_d ]
31    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
32        let (a, c, sigma, gamma) = (params[0], params[1], params[2], params[3]);
33        let sqrt2 = std::f64::consts::SQRT_2;
34        let dx = x[0] - c;
35        let inv_sigma = 1.0 / sigma;
36        let u = dx * inv_sigma; // dx/σ
37        let g = (-0.5 * u * u).exp();
38        let beta = gamma / (sigma * sqrt2); // γ/(σ√2)
39        let beta_dx = beta * dx;
40        let skew = 1.0 + libm::erf(beta_dx);
41        // erf_d = (2/√π)·exp(−(β·dx)²)
42        let erf_d = (2.0 / std::f64::consts::PI.sqrt()) * (-beta_dx * beta_dx).exp();
43
44        let da = g * skew;
45        // ∂f/∂c = A·g·[dx/σ²·skew − β·erf_d]
46        // Note: ∂g/∂c = g·u/σ (positive for dx>0, ∂(−½u²)/∂c = u·(−∂u/∂c) = u/σ);
47        //       ∂skew/∂c = erf_d·β·(∂dx/∂c) = erf_d·β·(−1)
48        let dc = a * g * (u * inv_sigma * skew - beta * erf_d);
49        // ∂f/∂σ: ∂g/∂σ = g·u²/σ; ∂(β·dx)/∂σ = −β·u
50        let ds = a * g * (u * u * inv_sigma * skew - beta * dx * inv_sigma * erf_d);
51        // ∂f/∂γ: ∂(β·dx)/∂γ = dx/(σ√2)
52        let dg = a * g * (dx / (sigma * sqrt2)) * erf_d;
53
54        vec![da, dc, ds, dg]
55    }
56
57    #[inline]
58    fn jacobian_into(&self, x: &[f64], params: &[f64], out: &mut [f64]) {
59        let jac = self.jacobian(x, params);
60        out[..4].copy_from_slice(&jac);
61    }
62
63    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
64        vec![
65            "amplitude".into(),
66            "center".into(),
67            "sigma".into(),
68            "gamma".into(),
69        ]
70    }
71}
72
73#[cfg(test)]
74mod tests {
75    use super::*;
76    use approx::assert_relative_eq;
77
78    #[test]
79    fn gamma_zero_reduces_to_gaussian() {
80        // erf(0)=0 ⇒ factor (1+0)=1 ⇒ pure Gaussian (value A at center).
81        let m = SkewedGaussian;
82        assert_relative_eq!(m.eval(&[0.0], &[3.0, 0.0, 1.0, 0.0]), 3.0, epsilon = 1e-12);
83    }
84
85    #[test]
86    fn positive_skew_lifts_high_side() {
87        // γ>0: the high-x side is enhanced relative to the symmetric mirror point.
88        let m = SkewedGaussian;
89        let hi = m.eval(&[1.0], &[1.0, 0.0, 1.0, 1.5]);
90        let lo = m.eval(&[-1.0], &[1.0, 0.0, 1.0, 1.5]);
91        assert!(hi > lo);
92    }
93
94    #[test]
95    fn jacobian_matches_central_difference_across_regimes() {
96        // Central difference is O(h^2); the old forward-difference pattern
97        // (see fano.rs::jacobian_finite_diff_check) was O(h), i.e. its own
98        // truncation error was the size of the 1e-4 tolerance it asserted
99        // against.
100        let m = SkewedGaussian;
101        let param_sets = [
102            [3.0, 0.0, 1.0, 1.5],   // nominal, positive skew
103            [1e-3, 0.0, 0.05, 0.5], // small amplitude and width
104            [5.0, -2.0, 3.0, -2.0], // offset centre, wide, negative skew
105        ];
106        for p in param_sets {
107            for &x in &[-3.0_f64, -0.5, 0.0, p[1], 0.5, 1.0, 7.0] {
108                let j = m.jacobian(&[x], &p);
109                for i in 0..p.len() {
110                    let h = 1e-6 * p[i].abs().max(1.0);
111                    let (mut a, mut b) = (p, p);
112                    a[i] += h;
113                    b[i] -= h;
114                    let fd = (m.eval(&[x], &a) - m.eval(&[x], &b)) / (2.0 * h);
115                    assert_relative_eq!(j[i], fd, epsilon = 1e-7, max_relative = 1e-6);
116                }
117            }
118        }
119    }
120}
spectrafit_models/skewed_gaussian.rs

spectrafit_models/
skewed_gaussian.rs
