spectrafit_models/
log_normal.rs

1use crate::Model;
2
3/// Log-normal peak: `A · exp(−(ln(x/c))² / (2σ²))` for `x > 0`, else `0`.
4///
5/// Parameters (in order): `[amplitude, center, sigma]`
6///
7/// - `amplitude` is the peak height attained at `x == center`.
8/// - `center > 0` is the peak location (log-space mode).
9/// - `sigma` is the log-space width.
10///
11/// The kernel is defined only for `x > 0`; at `x <= 0` it returns `0.0` (the
12/// log argument is undefined there). The numpy benchmark formula is identical —
13/// `np.where(x > 0, A·exp(−(ln(x/c))²/(2σ²)), 0)` — so numpy↔Rust parity is exact.
14pub struct LogNormal;
15
16impl Model for LogNormal {
17    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
18        let (a, c, sigma) = (params[0], params[1], params[2]);
19        if x[0] > 0.0 {
20            let l = (x[0] / c).ln();
21            let z = -(l * l) / (2.0 * sigma * sigma);
22            a * z.exp()
23        } else {
24            0.0
25        }
26    }
27
28    /// Central finite-difference Jacobian.
29    ///
30    /// The closed form has a logarithmic singularity at `x == 0`, so a numerical
31    /// (central-difference) Jacobian is used rather than an analytical one. Step
32    /// `h = 1e-7 · |p[i]|.max(1e-7)` (relative + absolute floor), matching the
33    /// trait's default magnitude.
34    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
35        let mut p = params.to_vec();
36        (0..params.len())
37            .map(|i| {
38                let h = 1e-7_f64 * params[i].abs().max(1e-7);
39                p[i] = params[i] + h;
40                let f_plus = self.eval(x, &p);
41                p[i] = params[i] - h;
42                let f_minus = self.eval(x, &p);
43                p[i] = params[i];
44                (f_plus - f_minus) / (2.0 * h)
45            })
46            .collect()
47    }
48
49    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
50        vec!["amplitude".into(), "center".into(), "sigma".into()]
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        // At x == center, ln(x/c) = 0 ⇒ exp(0) = 1, so value == amplitude.
62        let m = LogNormal;
63        let v = m.eval(&[2.5], &[3.0, 2.5, 0.4]);
64        assert_relative_eq!(v, 3.0, epsilon = 1e-12);
65    }
66
67    #[test]
68    fn eval_non_positive_is_zero() {
69        // x <= 0 (where ln is undefined) returns exactly 0.0.
70        let m = LogNormal;
71        assert_eq!(m.eval(&[0.0], &[1.0, 2.0, 0.5]), 0.0);
72        assert_eq!(m.eval(&[-1.0], &[1.0, 2.0, 0.5]), 0.0);
73    }
74
75    #[test]
76    fn param_names_are_canonical() {
77        assert_eq!(
78            LogNormal
79                .param_names()
80                .iter()
81                .map(|c| c.as_ref())
82                .collect::<Vec<_>>(),
83            &["amplitude", "center", "sigma"]
84        );
85    }
86
87    #[test]
88    fn jacobian_shape() {
89        let j = LogNormal.jacobian(&[1.5], &[2.0, 2.0, 0.6]);
90        assert_eq!(j.len(), 3);
91    }
92
93    #[test]
94    fn jacobian_amplitude_numerical_check() {
95        // ∂/∂amplitude at x == center is exp(0) = 1.
96        let m = LogNormal;
97        let j = m.jacobian(&[2.0], &[3.0, 2.0, 0.5]);
98        assert_relative_eq!(j[0], 1.0, epsilon = 1e-6);
99    }
100
101    #[test]
102    fn jacobian_center_zero_at_peak() {
103        // ∂/∂center at x == center is 0 (peak is stationary in log space).
104        let m = LogNormal;
105        let j = m.jacobian(&[2.0], &[3.0, 2.0, 0.5]);
106        assert_relative_eq!(j[1], 0.0, epsilon = 1e-6);
107    }
108}
spectrafit_models/log_normal.rs

spectrafit_models/
log_normal.rs
