sidereon_core/astro/math/
robust.rs

1//! Robust M-estimation primitives for iteratively reweighted least squares.
2//!
3//! These are the per-outer-iteration reweighting pieces of a Huber IRLS loop:
4//! a median-absolute-deviation scale estimate and the Huber weight function.
5//! They are deliberately pure `f64` arithmetic (abs, compare, divide, sort by
6//! [`f64::total_cmp`]) with no fused-multiply-add and no contraction, so the
7//! per-iteration weight vector is bit-reproducible against an explicit
8//! outer-loop reference recipe. The trust-region linear-algebra step that
9//! consumes the weights is BLAS-bound and is NOT a 0-ULP target.
10
11/// The default Huber tuning constant. Residuals scaled below this (in units of
12/// the robust scale) keep full weight; larger ones are down-weighted as
13/// `k / |u|`. `1.345` gives ~95% efficiency at the Gaussian model.
14pub const HUBER_K: f64 = 1.345;
15
16/// The MAD-to-sigma consistency constant for a normal distribution,
17/// `1 / Phi^-1(3/4)`. Multiplying the median absolute deviation by this makes
18/// it a consistent estimator of the standard deviation under normality.
19pub const MAD_NORMAL_CONST: f64 = 1.4826;
20
21#[derive(Debug, Clone, Copy, PartialEq, Eq, thiserror::Error)]
22pub enum RobustError {
23    #[error("invalid robust statistic {field}: {reason}")]
24    InvalidInput {
25        field: &'static str,
26        reason: &'static str,
27    },
28}
29
30impl RobustError {
31    pub const fn field(&self) -> &'static str {
32        match self {
33            Self::InvalidInput { field, .. } => field,
34        }
35    }
36
37    pub const fn reason(&self) -> &'static str {
38        match self {
39            Self::InvalidInput { reason, .. } => reason,
40        }
41    }
42}
43
44/// The median of `values`, computed on a `total_cmp` sort so the order (and
45/// thus the result for an even count, which averages the two central values)
46/// is deterministic. An empty slice yields `0.0`. The two central elements use
47/// `(a + b) / 2.0` when its addition is finite, with `a / 2.0 + b / 2.0` only
48/// as the overflow-safe fallback. Neither path uses FMA.
49pub fn median(values: &[f64]) -> Result<f64, RobustError> {
50    validate_finite_slice(values, "values")?;
51    if values.is_empty() {
52        return Ok(0.0);
53    }
54    let mut v: Vec<f64> = values.to_vec();
55    Ok(median_sorting_in_place(&mut v).unwrap_or(0.0))
56}
57
58/// The shared median kernel: sorts `values` in place by `total_cmp` and
59/// returns the middle element (odd count) or an overflow-safe average of the
60/// two central elements (even count, no FMA). `None` for an empty slice. No
61/// finiteness validation; callers own their input contracts.
62pub(crate) fn median_sorting_in_place(values: &mut [f64]) -> Option<f64> {
63    if values.is_empty() {
64        return None;
65    }
66    values.sort_by(|a, b| a.total_cmp(b));
67    let n = values.len();
68    if n % 2 == 1 {
69        Some(values[n / 2])
70    } else {
71        let lower = values[n / 2 - 1];
72        let upper = values[n / 2];
73        let sum = lower + upper;
74        Some(if sum.is_finite() {
75            sum / 2.0
76        } else {
77            lower / 2.0 + upper / 2.0
78        })
79    }
80}
81
82/// The median-absolute-deviation scale of `residuals`, scaled to a normal-sigma
83/// estimate and floored at `scale_floor`.
84///
85/// `s = max(scale_floor, MAD_NORMAL_CONST * median(|r_i - median(r)|))`. The
86/// floor prevents a near-perfect fit (MAD approaching zero) from blowing up the
87/// scaled residuals `u_i = r_i / s` and spuriously down-weighting every
88/// observation. Both medians use [`median`]'s `total_cmp` sort.
89pub fn mad_scale(residuals: &[f64], scale_floor: f64) -> Result<f64, RobustError> {
90    validate_finite_positive(scale_floor, "scale_floor")?;
91    let med = median(residuals)?;
92    let abs_dev: Vec<f64> = residuals.iter().map(|r| (r - med).abs()).collect();
93    let mad = median(&abs_dev)?;
94    let scaled = MAD_NORMAL_CONST * mad;
95    if scaled > scale_floor {
96        Ok(scaled)
97    } else {
98        Ok(scale_floor)
99    }
100}
101
102/// The Huber weight for a scaled residual `u = r / s`.
103///
104/// `w(u) = 1` for `|u| <= k` and `w(u) = k / |u|` otherwise (the Huber
105/// `psi(u) / u` form, always in `(0, 1]`). At `u == 0` the weight is `1`. This
106/// is the multiplier applied on top of any base (elevation) weight to obtain the
107/// effective per-observation weight of the current outer iteration.
108pub fn huber_weight(u: f64, k: f64) -> f64 {
109    let a = u.abs();
110    if a <= k {
111        1.0
112    } else {
113        k / a
114    }
115}
116
117fn validate_finite_slice(values: &[f64], field: &'static str) -> Result<(), RobustError> {
118    if values.iter().all(|value| value.is_finite()) {
119        Ok(())
120    } else {
121        Err(invalid_input(field, "not finite"))
122    }
123}
124
125fn validate_finite_positive(value: f64, field: &'static str) -> Result<(), RobustError> {
126    if !value.is_finite() {
127        Err(invalid_input(field, "not finite"))
128    } else if value <= 0.0 {
129        Err(invalid_input(field, "not positive"))
130    } else {
131        Ok(())
132    }
133}
134
135fn invalid_input(field: &'static str, reason: &'static str) -> RobustError {
136    RobustError::InvalidInput { field, reason }
137}
138
139#[cfg(test)]
140mod tests {
141    use super::*;
142
143    #[test]
144    fn median_odd_even() {
145        assert_eq!(median(&[3.0, 1.0, 2.0]).unwrap(), 2.0);
146        assert_eq!(median(&[1.0, 2.0, 3.0, 4.0]).unwrap(), 2.5);
147        assert_eq!(median(&[]).unwrap(), 0.0);
148    }
149
150    #[test]
151    fn median_does_not_overflow_for_finite_central_values() {
152        assert_eq!(median(&[f64::MAX, f64::MAX]).unwrap(), f64::MAX);
153        assert_eq!(median(&[-f64::MAX, -f64::MAX]).unwrap(), -f64::MAX);
154    }
155
156    #[test]
157    fn median_rejects_nonfinite_sample() {
158        assert_eq!(
159            median(&[1.0, f64::NAN]),
160            Err(RobustError::InvalidInput {
161                field: "values",
162                reason: "not finite"
163            })
164        );
165    }
166
167    #[test]
168    fn huber_weight_breaks_at_k() {
169        assert_eq!(huber_weight(0.0, HUBER_K), 1.0);
170        assert_eq!(huber_weight(HUBER_K, HUBER_K), 1.0);
171        let w = huber_weight(2.0 * HUBER_K, HUBER_K);
172        assert!((w - 0.5).abs() < 1e-15);
173    }
174
175    #[test]
176    fn mad_scale_floored() {
177        // All-equal residuals give MAD 0, so the floor governs.
178        assert_eq!(mad_scale(&[5.0, 5.0, 5.0], 0.25).unwrap(), 0.25);
179    }
180
181    #[test]
182    fn mad_scale_rejects_nonfinite_sample() {
183        assert_eq!(
184            mad_scale(&[5.0, f64::INFINITY], 0.25),
185            Err(RobustError::InvalidInput {
186                field: "values",
187                reason: "not finite"
188            })
189        );
190    }
191}
sidereon_core/astro/math/robust.rs

sidereon_core/astro/math/
robust.rs
