spectrafit_models/
gaussian2d.rs

1//! 2-D Gaussian kernel scaffold.
2//!
3//! This is **scaffolding** for 2-D / N-D fitting.  The `eval` and analytical
4//! Jacobian bodies are implemented for the axis-aligned (no rotation) MVP so
5//! the kernel compiles and can be exercised by `#[ignore]` round-trip /
6//! finite-difference tests.  Rotation (`theta`) is intentionally omitted from
7//! the MVP — see the parameter note below.
8
9use crate::Model;
10
11/// Axis-aligned 2-D Gaussian peak.
12///
13/// Formula (no rotation):
14///
15/// ```text
16/// f(x, y) = A · exp( −(x − cₓ)² / (2·σₓ²) − (y − c_y)² / (2·σ_y²) )
17/// ```
18///
19/// Parameters (in order):
20/// `[amplitude, center_x, center_y, sigma_x, sigma_y]`
21///
22/// - `amplitude` — peak value at `(center_x, center_y)` (not area).
23/// - `center_x`, `center_y` — peak location.
24/// - `sigma_x`, `sigma_y` — standard deviations along each axis
25///   (FWHM ≈ 2.355·σ).
26///
27/// MVP note: a `theta` rotation parameter is intentionally **omitted** to keep
28/// the first 2-D kernel minimal; rotated Gaussians are a follow-up unit.
29pub struct Gaussian2D;
30
31impl Model for Gaussian2D {
32    fn eval(&self, x: &[f64], params: &[f64]) -> f64 {
33        let (a, cx, cy, sx, sy) = (params[0], params[1], params[2], params[3], params[4]);
34        let dx = x[0] - cx;
35        let dy = x[1] - cy;
36        let z = -(dx * dx) / (2.0 * sx * sx) - (dy * dy) / (2.0 * sy * sy);
37        a * z.exp()
38    }
39
40    fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64> {
41        let mut out = vec![0.0_f64; self.param_names().len()];
42        self.jacobian_into(x, params, &mut out);
43        out
44    }
45
46    #[inline]
47    fn jacobian_into(&self, x: &[f64], params: &[f64], out: &mut [f64]) {
48        let (a, cx, cy, sx, sy) = (params[0], params[1], params[2], params[3], params[4]);
49        let dx = x[0] - cx;
50        let dy = x[1] - cy;
51        let sx2 = sx * sx;
52        let sy2 = sy * sy;
53        let g = (-(dx * dx) / (2.0 * sx2) - (dy * dy) / (2.0 * sy2)).exp();
54
55        // ∂/∂amplitude = g
56        out[0] = g;
57        // ∂/∂center_x  = A · g · (x − cₓ) / σₓ²
58        out[1] = a * g * dx / sx2;
59        // ∂/∂center_y  = A · g · (y − c_y) / σ_y²
60        out[2] = a * g * dy / sy2;
61        // ∂/∂sigma_x   = A · g · (x − cₓ)² / σₓ³
62        out[3] = a * g * dx * dx / (sx2 * sx);
63        // ∂/∂sigma_y   = A · g · (y − c_y)² / σ_y³
64        out[4] = a * g * dy * dy / (sy2 * sy);
65    }
66
67    fn param_names(&self) -> Vec<std::borrow::Cow<'static, str>> {
68        vec![
69            "amplitude".into(),
70            "center_x".into(),
71            "center_y".into(),
72            "sigma_x".into(),
73            "sigma_y".into(),
74        ]
75    }
76
77    fn n_dims(&self) -> usize {
78        2
79    }
80}
81
82#[cfg(test)]
83mod tests {
84    use super::*;
85    use approx::assert_relative_eq;
86
87    #[test]
88    fn n_dims_is_two() {
89        assert_eq!(Gaussian2D.n_dims(), 2);
90    }
91
92    #[test]
93    fn param_names_match_convention() {
94        assert_eq!(
95            Gaussian2D
96                .param_names()
97                .iter()
98                .map(|c| c.as_ref())
99                .collect::<Vec<_>>(),
100            &["amplitude", "center_x", "center_y", "sigma_x", "sigma_y"]
101        );
102    }
103
104    #[test]
105    fn eval_at_center_equals_amplitude() {
106        // At (cx, cy) the 2-D Gaussian equals amplitude.
107        let v = Gaussian2D.eval(&[0.5, -1.0], &[3.0, 0.5, -1.0, 1.0, 2.0]);
108        assert_relative_eq!(v, 3.0, epsilon = 1e-12);
109    }
110
111    #[test]
112    fn jacobian_into_matches_finite_difference() {
113        // Upgraded from a single-point forward difference (h=1e-6,
114        // epsilon=1e-5) to a central difference (O(h^2)) swept over three
115        // parameter regimes and several (x, y) offsets from each regime's
116        // own centre.
117        let m = Gaussian2D;
118        let param_sets = [
119            [2.0, 0.5, -1.0, 1.0, 1.5],   // nominal
120            [1e-3, 0.0, 0.0, 0.05, 0.05], // small amplitude/width
121            [5.0, -2.0, 3.0, 2.5, 1.5],   // offset centre, wide
122        ];
123        for params in param_sets {
124            let (cx, cy) = (params[1], params[2]);
125            for &(ox, oy) in &[
126                (0.0_f64, 0.0),
127                (0.7, -0.3),
128                (-2.5, 2.0),
129                (3.0, -1.5),
130                (0.05, 0.05),
131            ] {
132                let x = [cx + ox, cy + oy];
133                let mut analytic = vec![0.0_f64; params.len()];
134                m.jacobian_into(&x, &params, &mut analytic);
135                for i in 0..params.len() {
136                    let h = 1e-6 * params[i].abs().max(1.0);
137                    let (mut a, mut b) = (params, params);
138                    a[i] += h;
139                    b[i] -= h;
140                    let fd = (m.eval(&x, &a) - m.eval(&x, &b)) / (2.0 * h);
141                    assert_relative_eq!(analytic[i], fd, epsilon = 1e-7, max_relative = 1e-6);
142                }
143            }
144        }
145    }
146}
spectrafit_models/gaussian2d.rs

spectrafit_models/
gaussian2d.rs
