spectrafit_models/
cauchy_dispersion.rs

1use crate::Model;
2
3/// Cauchy refractive-index dispersion: `n(x) = a + b/x² + c/x⁴`.
4///
5/// Parameters (in order): `[a, b, c]`
6///
7/// - `a` is the high-frequency (constant) refractive-index offset.
8/// - `b` is the first dispersion coefficient (units of x²).
9/// - `c` is the second dispersion coefficient (units of x⁴).
10///
11/// The independent variable `x` is a wavelength and must be `> 0`; the kernel is smooth
12/// and analytic there. At `x == 0` the `1/x²` / `1/x⁴` terms are undefined, so the kernel
13/// returns `0.0` (matching the numpy oracle `np.where(x > 0, a + b/x² + c/x⁴, 0)`), keeping
14/// numpy↔Rust parity exact. The analytical Jacobian is `[1, 1/x², 1/x⁴]`.
15pub struct CauchyDispersion;
16
17impl Model for CauchyDispersion {
18    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
19        let (a, b, c) = (params[0], params[1], params[2]);
20        if x[0] > 0.0 {
21            let x2 = x[0] * x[0];
22            let x4 = x2 * x2;
23            a + b / x2 + c / x4
24        } else {
25            0.0
26        }
27    }
28
29    /// Analytical Jacobian: `∂n/∂a = 1`, `∂n/∂b = 1/x²`, `∂n/∂c = 1/x⁴` (for `x > 0`);
30    /// all zero at `x ≤ 0` where the kernel is clamped to `0`.
31    fn jacobian(&self, x: &[f64], _params: &[f64]) -> Vec<f64> {
32        if x[0] > 0.0 {
33            let x2 = x[0] * x[0];
34            let x4 = x2 * x2;
35            vec![1.0, 1.0 / x2, 1.0 / x4]
36        } else {
37            vec![0.0, 0.0, 0.0]
38        }
39    }
40
41    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
42        vec!["a".into(), "b".into(), "c".into()]
43    }
44}
45
46#[cfg(test)]
47mod tests {
48    use super::*;
49    use approx::assert_relative_eq;
50
51    #[test]
52    fn dispersion_value() {
53        // a + b/x² + c/x⁴ at x = 2: 1.5 + 0.4/4 + 0.2/16 = 1.5 + 0.1 + 0.0125.
54        let m = CauchyDispersion;
55        assert_relative_eq!(
56            m.eval(&[2.0], &[1.5, 0.4, 0.2]),
57            1.5 + 0.1 + 0.0125,
58            epsilon = 1e-12
59        );
60    }
61
62    #[test]
63    fn non_positive_x_is_zero() {
64        let m = CauchyDispersion;
65        assert_eq!(m.eval(&[0.0], &[1.5, 0.4, 0.2]), 0.0);
66        assert_eq!(m.eval(&[-1.0], &[1.5, 0.4, 0.2]), 0.0);
67    }
68
69    #[test]
70    fn analytical_jacobian_matches_closed_form() {
71        let m = CauchyDispersion;
72        let j = m.jacobian(&[2.0], &[1.5, 0.4, 0.2]);
73        assert_eq!(j.len(), 3);
74        assert_relative_eq!(j[0], 1.0, epsilon = 1e-12);
75        assert_relative_eq!(j[1], 1.0 / 4.0, epsilon = 1e-12);
76        assert_relative_eq!(j[2], 1.0 / 16.0, epsilon = 1e-12);
77    }
78
79    #[test]
80    fn param_names_are_canonical() {
81        assert_eq!(
82            CauchyDispersion
83                .param_names()
84                .iter()
85                .map(|c| c.as_ref())
86                .collect::<Vec<_>>(),
87            &["a", "b", "c"]
88        );
89    }
90
91    #[test]
92    fn jacobian_matches_central_difference_across_regimes() {
93        // Central difference is O(h^2); the old forward-difference pattern
94        // (see fano.rs::jacobian_finite_diff_check) was O(h), i.e. its own
95        // truncation error was the size of the 1e-4 tolerance it asserted
96        // against.
97        let m = CauchyDispersion;
98        let param_sets = [
99            [1.5, 0.4, 0.2],    // nominal
100            [1e-3, 1e-3, 1e-3], // small coefficients
101            [5.0, -2.0, 8.0],   // large / negative coefficients
102        ];
103        for p in param_sets {
104            // x <= 0 clamps eval (and jacobian) to 0.0 regardless of params, so
105            // the central difference agrees trivially there (0 == 0), proving
106            // the clamp itself doesn't leak a spurious derivative. x = 0.3 is
107            // the smallest positive point kept: below that, c/x^4 for the
108            // large-coefficient regime dwarfs the ~1e-6-scaled perturbation of
109            // `a`, and the resulting catastrophic cancellation in the finite
110            // difference (not the analytic Jacobian, which is exactly
111            // param-independent here) pushes its own error above 1e-7.
112            for &x in &[0.3_f64, 0.5, 1.0, 2.0, 7.0, 0.0, -1.0] {
113                let j = m.jacobian(&[x], &p);
114                for i in 0..p.len() {
115                    let h = 1e-6 * p[i].abs().max(1.0);
116                    let (mut a, mut b) = (p, p);
117                    a[i] += h;
118                    b[i] -= h;
119                    let fd = (m.eval(&[x], &a) - m.eval(&[x], &b)) / (2.0 * h);
120                    assert_relative_eq!(j[i], fd, epsilon = 1e-7, max_relative = 1e-6);
121                }
122            }
123        }
124    }
125}
spectrafit_models/cauchy_dispersion.rs

spectrafit_models/
cauchy_dispersion.rs
