spectrafit_models/
asym_ir.rs

1use crate::Model;
2
3/// Asymmetric IR band: a Gaussian multiplied by a logistic sigmoid.
4///
5/// `f(x) = A·exp(−(x − c)²/(2σ²)) · 1/(1 + exp(−k·(x − c)))`.
6///
7/// Parameters (in order): `[amplitude, center, sigma, k]`. `amplitude` is the Gaussian scale
8/// (peak ≈ A/2 at center because of the sigmoid); `k` is the asymmetry. The sigmoid exponent
9/// is clamped to ≤ 50 to avoid overflow — the numpy oracle clamps identically
10/// (`np.clip(-k*(x-c), None, 50.0)`), so numpy↔Rust parity is exact.
11pub struct AsymIr;
12
13impl Model for AsymIr {
14    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
15        let (a, c, sigma, k) = (params[0], params[1], params[2], params[3]);
16        let dx = x[0] - c;
17        let g = a * (-(dx * dx) / (2.0 * sigma * sigma)).exp();
18        let arg = (-(k * dx)).min(50.0);
19        g / (1.0 + arg.exp())
20    }
21
22    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
23        let mut p = params.to_vec();
24        (0..params.len())
25            .map(|i| {
26                let h = 1e-7_f64 * params[i].abs().max(1e-7);
27                p[i] = params[i] + h;
28                let f_plus = self.eval(x, &p);
29                p[i] = params[i] - h;
30                let f_minus = self.eval(x, &p);
31                p[i] = params[i];
32                (f_plus - f_minus) / (2.0 * h)
33            })
34            .collect()
35    }
36
37    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
38        vec![
39            "amplitude".into(),
40            "center".into(),
41            "sigma".into(),
42            "k".into(),
43        ]
44    }
45}
46
47#[cfg(test)]
48mod tests {
49    use super::*;
50    use approx::assert_relative_eq;
51
52    #[test]
53    fn eval_at_center_is_half_amplitude() {
54        // At x == center, dx = 0 ⇒ Gaussian = A, sigmoid = 1/(1+1) = 0.5.
55        let v = AsymIr.eval(&[1.5], &[4.0, 1.5, 0.8, 1.0]);
56        assert_relative_eq!(v, 2.0, epsilon = 1e-12);
57    }
58
59    #[test]
60    fn param_names_and_jacobian() {
61        assert_eq!(
62            AsymIr
63                .param_names()
64                .iter()
65                .map(|c| c.as_ref())
66                .collect::<Vec<_>>(),
67            &["amplitude", "center", "sigma", "k"]
68        );
69        assert_eq!(AsymIr.jacobian(&[1.0], &[4.0, 1.5, 0.8, 1.0]).len(), 4);
70    }
71
72    // ----- Limiting-case asymptotic (ground-truth verification) -----
73    //
74    // At k = 0 the sigmoid collapses to a constant 1/2 (1/(1+exp(0))):
75    //
76    //     AsymIr(A, c, σ, k=0)  ≡  Gaussian(A/2, c, σ)   at every x
77    //
78    // Catches: missing/wrong constant factor in the no-asymmetry case,
79    // sigmoid-evaluation bugs that break the k=0 symmetric reduction.
80
81    #[test]
82    fn k_zero_equals_half_gaussian() {
83        use crate::gaussian::Gaussian;
84        let asym = AsymIr;
85        let gauss = Gaussian;
86        let a = 5.0_f64;
87        let c = 0.3_f64;
88        let sigma = 0.8_f64;
89        let p_asym = [a, c, sigma, 0.0]; // k = 0 → symmetric
90        let p_gauss = [a / 2.0, c, sigma];
91        for &xi in &[-2.0_f64, -0.5, c, 0.5, 2.0] {
92            assert_relative_eq!(
93                asym.eval(&[xi], &p_asym),
94                gauss.eval(&[xi], &p_gauss),
95                epsilon = 1e-12
96            );
97        }
98    }
99}
spectrafit_models/asym_ir.rs

spectrafit_models/
asym_ir.rs
