spectrafit_models/
kww.rs

1use crate::Model;
2
3/// Kohlrausch–Williams–Watts (KWW) stretched exponential: `A · exp(−(x/τ)^β)`.
4///
5/// Parameters (in order): `[amplitude, tau, beta]`
6///
7/// - `amplitude` (`A`) is the value at `x == 0` (where `(0/τ)^β = 0` ⇒ `exp(0) = 1`).
8/// - `tau` (`τ > 0`) is the characteristic relaxation time.
9/// - `beta` (`0 < β ≤ 1`) is the stretching exponent: `β = 1` recovers a plain
10///   exponential, `β < 1` gives a stretched (multi-timescale) relaxation.
11///
12/// Defined for `x ≥ 0`. For `x < 0` the base `x/τ` is negative and a fractional `β`
13/// would yield a NaN, so the kernel returns `0.0` there. The numpy benchmark formula is
14/// identical — `np.where(x >= 0, A·exp(−(x/τ)^β), 0)` — so numpy↔Rust parity is exact.
15pub struct Kww;
16
17impl Model for Kww {
18    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
19        let (a, tau, beta) = (params[0], params[1], params[2]);
20        if x[0] >= 0.0 {
21            a * (-(x[0] / tau).powf(beta)).exp()
22        } else {
23            0.0
24        }
25    }
26
27    /// Central finite-difference Jacobian. `∂/∂β` involves `ln(x/τ)` which diverges as
28    /// `x → 0⁺`, so a numerical Jacobian is used (matching `log_normal`).
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!["amplitude".into(), "tau".into(), "beta".into()]
46    }
47}
48
49#[cfg(test)]
50mod tests {
51    use super::*;
52    use approx::assert_relative_eq;
53
54    #[test]
55    fn value_at_zero_equals_amplitude() {
56        // (0/τ)^β = 0 ⇒ exp(0) = 1 ⇒ value == amplitude.
57        let m = Kww;
58        assert_relative_eq!(m.eval(&[0.0], &[3.0, 2.0, 0.7]), 3.0, epsilon = 1e-12);
59    }
60
61    #[test]
62    fn beta_one_is_plain_exponential() {
63        // β = 1 ⇒ A·exp(−x/τ). At x = τ that is A·exp(−1).
64        let m = Kww;
65        assert_relative_eq!(
66            m.eval(&[2.0], &[3.0, 2.0, 1.0]),
67            3.0 * (-1.0_f64).exp(),
68            epsilon = 1e-12
69        );
70    }
71
72    #[test]
73    fn negative_x_is_zero() {
74        let m = Kww;
75        assert_eq!(m.eval(&[-1.0], &[3.0, 2.0, 0.7]), 0.0);
76    }
77
78    #[test]
79    fn param_names_are_canonical() {
80        assert_eq!(
81            Kww.param_names()
82                .iter()
83                .map(|c| c.as_ref())
84                .collect::<Vec<_>>(),
85            &["amplitude", "tau", "beta"]
86        );
87    }
88
89    #[test]
90    fn jacobian_shape_and_amplitude() {
91        let m = Kww;
92        let j = m.jacobian(&[2.0], &[3.0, 2.0, 0.7]);
93        assert_eq!(j.len(), 3);
94        // ∂/∂amplitude = exp(−(x/τ)^β).
95        let expected = (-(2.0_f64 / 2.0).powf(0.7)).exp();
96        assert_relative_eq!(j[0], expected, epsilon = 1e-6);
97    }
98
99    // ----- Limiting-case asymptotic (ground-truth verification) -----
100    //
101    // KWW reduces to a plain single exponential when β = 1:
102    //
103    //     A · exp(−(x/τ)^β)  with β = 1  ≡  A · exp(−x/τ)
104    //
105    // Test at several x to pin the identity, not just the happy-path point.
106
107    #[test]
108    fn beta_one_equals_single_exponential() {
109        let m = Kww;
110        let a = 2.5_f64;
111        let tau = 1.7_f64;
112        let p = [a, tau, 1.0]; // β = 1 collapse
113        for &xi in &[0.0_f64, 0.5, 1.0, 2.0, 5.0] {
114            let kww_val = m.eval(&[xi], &p);
115            let single_exp = a * (-xi / tau).exp();
116            assert_relative_eq!(kww_val, single_exp, epsilon = 1e-12);
117        }
118    }
119}
spectrafit_models/kww.rs

spectrafit_models/
kww.rs
