ogeom_algo/
cutoff.rs

1//! Which points a fit keeps: the rule that drops outliers between rounds of
2//! a fit, shared by every fit that iterates over a point set.
3
4use ogeom_core::{OgeomResult, ogeom_bail};
5
6/// Which points a fit keeps.
7///
8/// A fit keeps the points whose residual is no larger than a threshold, fits
9/// again on those, and repeats until the kept set stops changing. The
10/// threshold is measured over the points the last round kept.
11#[derive(Debug, Clone, Copy, PartialEq, Default)]
12pub enum Cutoff {
13    /// Every point, in one round.
14    #[default]
15    All,
16    /// The points within `k` times the RMS residual.
17    Sigma(f64),
18    /// The points within `k` times the median residual size.
19    Median(f64),
20    /// The points within a fixed distance.
21    Distance(f64),
22}
23
24impl Cutoff {
25    /// The rule's number checked: finite and positive.
26    ///
27    /// # Errors
28    ///
29    /// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if the factor or
30    /// the distance is not finite and positive.
31    pub fn check(self) -> OgeomResult<()> {
32        let (name, value) = match self {
33            Self::All => return Ok(()),
34            Self::Sigma(k) => ("sigma factor", k),
35            Self::Median(k) => ("median factor", k),
36            Self::Distance(d) => ("cutoff distance", d),
37        };
38        if !(value.is_finite() && value > 0.0) {
39            ogeom_bail!(Domain, "the {name} {value} must be finite and positive");
40        }
41        Ok(())
42    }
43
44    /// The largest residual size kept, measured over the residuals of the
45    /// points the last round kept; infinite for [`Cutoff::All`] and where
46    /// there are none.
47    #[must_use]
48    pub fn threshold(self, kept: &[f64]) -> f64 {
49        if kept.is_empty() {
50            return f64::INFINITY;
51        }
52        #[allow(
53            clippy::cast_precision_loss,
54            reason = "point counts are far below 2^52"
55        )]
56        let count = kept.len() as f64;
57        match self {
58            Self::All => f64::INFINITY,
59            Self::Sigma(k) => k * (kept.iter().map(|r| r * r).sum::<f64>() / count).sqrt(),
60            Self::Median(k) => {
61                let mut sizes: Vec<f64> = kept.iter().map(|r| r.abs()).collect();
62                sizes.sort_by(f64::total_cmp);
63                let mid = sizes.len() / 2;
64                let median = if sizes.len().is_multiple_of(2) {
65                    f64::midpoint(sizes[mid - 1], sizes[mid])
66                } else {
67                    sizes[mid]
68                };
69                k * median
70            }
71            Self::Distance(d) => d,
72        }
73    }
74
75    /// As [`Cutoff::threshold`], each residual counting by its weight: the
76    /// RMS and the median of a distribution of weight rather than of
77    /// count, so a fit over triangles weighs them by area. Equal weights
78    /// give the unweighted threshold.
79    pub(crate) fn threshold_weighted(self, kept: &[(f64, f64)]) -> f64 {
80        let total: f64 = kept.iter().map(|(_, w)| w).sum();
81        if kept.is_empty() || total.is_nan() || total <= 0.0 {
82            return f64::INFINITY;
83        }
84        match self {
85            Self::All => f64::INFINITY,
86            Self::Sigma(k) => k * (kept.iter().map(|(r, w)| w * r * r).sum::<f64>() / total).sqrt(),
87            Self::Median(k) => {
88                let mut sizes: Vec<(f64, f64)> = kept.iter().map(|(r, w)| (r.abs(), *w)).collect();
89                sizes.sort_by(|p, q| p.0.total_cmp(&q.0));
90                let half = total * 0.5;
91                let mut below = 0.0;
92                for (i, (size, w)) in sizes.iter().enumerate() {
93                    below += w;
94                    if below > half {
95                        return k * size;
96                    }
97                    if below == half {
98                        return k * sizes
99                            .get(i + 1)
100                            .map_or(*size, |next| f64::midpoint(*size, next.0));
101                    }
102                }
103                k * sizes.last().map_or(0.0, |s| s.0)
104            }
105            Self::Distance(d) => d,
106        }
107    }
108}
109
110#[cfg(test)]
111mod tests {
112    use super::*;
113
114    #[test]
115    fn each_rule_measures_its_threshold_over_the_kept_residuals() {
116        let kept = [3.0, -4.0, 1.0, -2.0];
117        assert!(Cutoff::All.threshold(&kept).is_infinite());
118        let rms = (30.0_f64 / 4.0).sqrt();
119        assert!((Cutoff::Sigma(2.0).threshold(&kept) - 2.0 * rms).abs() < 1e-15);
120        assert!((Cutoff::Median(3.0).threshold(&kept) - 7.5).abs() < 1e-15);
121        assert!((Cutoff::Distance(0.5).threshold(&kept) - 0.5).abs() < 1e-15);
122    }
123
124    #[test]
125    fn equal_weights_give_the_unweighted_threshold_and_weights_move_it() {
126        let kept = [3.0, -4.0, 1.0, -2.0];
127        let even: Vec<(f64, f64)> = kept.iter().map(|r| (*r, 0.5)).collect();
128        for rule in [
129            Cutoff::Sigma(2.0),
130            Cutoff::Median(3.0),
131            Cutoff::Distance(0.5),
132        ] {
133            let (a, b) = (rule.threshold(&kept), rule.threshold_weighted(&even));
134            assert!((a - b).abs() <= 1e-15 * a, "{rule:?}: {a} and {b}");
135        }
136        let odd = [(1.0, 1.0), (2.0, 1.0), (9.0, 1.0)];
137        assert!((Cutoff::Median(1.0).threshold_weighted(&odd) - 2.0).abs() < 1e-15);
138        // The large residual weighs as much as the others together and
139        // more: the median moves to it, and the RMS follows.
140        let heavy = [(1.0, 1.0), (2.0, 1.0), (9.0, 3.0)];
141        assert!((Cutoff::Median(1.0).threshold_weighted(&heavy) - 9.0).abs() < 1e-15);
142        let rms = ((1.0 + 4.0 + 3.0 * 81.0) / 5.0_f64).sqrt();
143        assert!((Cutoff::Sigma(1.0).threshold_weighted(&heavy) - rms).abs() < 1e-14);
144        assert!(Cutoff::All.threshold_weighted(&heavy).is_infinite());
145    }
146
147    #[test]
148    fn a_rule_without_a_positive_number_is_refused() {
149        for bad in [0.0, -1.0, f64::NAN, f64::INFINITY] {
150            assert!(Cutoff::Sigma(bad).check().is_err());
151            assert!(Cutoff::Median(bad).check().is_err());
152            assert!(Cutoff::Distance(bad).check().is_err());
153        }
154        assert!(Cutoff::All.check().is_ok());
155        assert!(Cutoff::Sigma(3.0).check().is_ok());
156    }
157}
ogeom_algo/cutoff.rs

ogeom_algo/
cutoff.rs
