geop_core_math/solvers/
bisect.rs

1//! Bisection of a bracket on which side of a crossing its middle lies —
2//! in plain numbers, where where-exactly is a free choice found as far as
3//! `f64` goes, or in sharp scalars, where the side is told by an interval
4//! test that may not decide — and the scan of samples that finds the
5//! brackets: every pair of neighbouring samples on different sides.
6
7use crate::{
8    geop_error::{GeopError, GeopResult},
9    scalars::Scalar,
10};
11
12/// A value a bracket can be bisected in.
13pub trait Bisect: Copy {
14    /// The middle of `lo` and `hi`: where the next side is asked.
15    fn middle(lo: Self, hi: Self) -> GeopResult<Self>;
16    /// Whether `mid` lies strictly between `lo` and `hi` — once it does not,
17    /// the bracket cannot shrink any further.
18    fn strictly_between(lo: Self, mid: Self, hi: Self) -> bool;
19}
20
21impl Bisect for f64 {
22    fn middle(lo: Self, hi: Self) -> GeopResult<Self> {
23        Ok(0.5 * (lo + hi))
24    }
25    fn strictly_between(lo: Self, mid: Self, hi: Self) -> bool {
26        mid > lo.min(hi) && mid < lo.max(hi)
27    }
28}
29
30/// Sharp: which point is asked is the search's own free choice, and a
31/// sharp middle keeps the bracket from widening as it shrinks.
32impl Bisect for Scalar {
33    fn middle(lo: Self, hi: Self) -> GeopResult<Self> {
34        Ok(lo.add(hi).div(Scalar::TWO)?.sharpen())
35    }
36    /// For `lo` below `hi`.
37    fn strictly_between(lo: Self, mid: Self, hi: Self) -> bool {
38        lo.definitely_less(mid) && mid.definitely_less(hi)
39    }
40}
41
42/// Where a bisection ended: its bracket, and the middle whose side could
43/// not be told, if that is what stopped it.
44#[derive(Clone, Copy, Debug)]
45pub struct Bracket<T> {
46    pub lo: T,
47    pub hi: T,
48    pub undecided: Option<T>,
49}
50
51/// Bisects `(lo, hi)` at most `steps` times. `side(mid)` says whether the
52/// middle lies on `lo`'s side (`Some(true)`, it becomes `lo`), on `hi`'s
53/// (`Some(false)`), or cannot be told (`None`, which stops the bisection).
54/// With `until_collapsed`, it also stops once the middle no longer lies
55/// strictly inside the bracket. `side` may keep state of its own — a seed
56/// carried from one middle to the next — since every step asks it once.
57pub fn bisect<T: Bisect, E: From<GeopError>>(
58    lo: T,
59    hi: T,
60    steps: usize,
61    until_collapsed: bool,
62    mut side: impl FnMut(T) -> Result<Option<bool>, E>,
63) -> Result<Bracket<T>, E> {
64    let (mut lo, mut hi) = (lo, hi);
65    for _ in 0..steps {
66        let mid = T::middle(lo, hi)?;
67        if until_collapsed && !T::strictly_between(lo, mid, hi) {
68            break;
69        }
70        match side(mid)? {
71            Some(true) => lo = mid,
72            Some(false) => hi = mid,
73            None => {
74                return Ok(Bracket {
75                    lo,
76                    hi,
77                    undecided: Some(mid),
78                });
79            }
80        }
81    }
82    Ok(Bracket {
83        lo,
84        hi,
85        undecided: None,
86    })
87}
88
89/// Neighbouring samples on different sides: see [`sign_changes`].
90pub struct SignChanges<I: Iterator, F> {
91    samples: I,
92    side: F,
93    /// The last sample whose side was told, and its side.
94    last: Option<(I::Item, bool)>,
95}
96
97/// The pairs of neighbouring samples of `samples` that `side` puts on
98/// different sides, in order, each with the side of its first: brackets
99/// of a crossing. `side` says which side a sample lies on (`< cut` against
100/// `<= 0`, a sign of a number or of an interval: the caller's choice), or
101/// `None` where it cannot tell; such a sample is passed over, its
102/// neighbours paired across it. Lazy: a scan for the first crossing asks
103/// no sample past it. [`SignChanges::bisected`] bisects each bracket.
104pub fn sign_changes<I: IntoIterator, F>(samples: I, side: F) -> SignChanges<I::IntoIter, F> {
105    SignChanges {
106        samples: samples.into_iter(),
107        side,
108        last: None,
109    }
110}
111
112impl<T: Copy, E, I: Iterator<Item = T>, F: FnMut(T) -> Result<Option<bool>, E>> Iterator
113    for SignChanges<I, F>
114{
115    type Item = Result<(T, T, bool), E>;
116
117    fn next(&mut self) -> Option<Self::Item> {
118        for t in self.samples.by_ref() {
119            let s = match (self.side)(t) {
120                Ok(Some(s)) => s,
121                Ok(None) => continue,
122                Err(e) => return Some(Err(e)),
123            };
124            if let Some((a, side_a)) = self.last.replace((t, s))
125                && side_a != s
126            {
127                return Some(Ok((a, t, side_a)));
128            }
129        }
130        None
131    }
132}
133
134impl<I: Iterator, F> SignChanges<I, F> {
135    /// Each bracket bisected at most `steps` times on the same `side`
136    /// (see [`bisect`] for `until_collapsed`).
137    pub fn bisected(self, steps: usize, until_collapsed: bool) -> Crossings<I, F> {
138        Crossings {
139            scan: self,
140            steps,
141            until_collapsed,
142        }
143    }
144}
145
146/// The brackets of a [`SignChanges`] scan, bisected: see
147/// [`SignChanges::bisected`].
148pub struct Crossings<I: Iterator, F> {
149    scan: SignChanges<I, F>,
150    steps: usize,
151    until_collapsed: bool,
152}
153
154impl<T: Bisect, E: From<GeopError>, I: Iterator<Item = T>, F: FnMut(T) -> Result<Option<bool>, E>>
155    Iterator for Crossings<I, F>
156{
157    type Item = Result<Bracket<T>, E>;
158
159    fn next(&mut self) -> Option<Self::Item> {
160        let (a, b, side_a) = match self.scan.next()? {
161            Ok(bracket) => bracket,
162            Err(e) => return Some(Err(e)),
163        };
164        let side = &mut self.scan.side;
165        Some(bisect(a, b, self.steps, self.until_collapsed, |mid| {
166            Ok(side(mid)?.map(|s| s == side_a))
167        }))
168    }
169}
170
171#[cfg(test)]
172mod tests {
173    use super::*;
174
175    #[test]
176    fn bisects_to_the_last_bit() {
177        let b = bisect::<f64, GeopError>(1.0, 2.0, 200, true, |x| Ok(Some(x * x < 2.0))).unwrap();
178        assert_eq!(b.hi, f64::from_bits(b.lo.to_bits() + 1));
179        assert!(b.lo * b.lo < 2.0 && b.hi * b.hi >= 2.0);
180    }
181
182    /// A middle whose side cannot be told stops the bisection and is
183    /// reported, the bracket left as it was told.
184    #[test]
185    fn an_undecided_middle_is_reported() {
186        let b =
187            bisect::<f64, GeopError>(0.0, 1.0, 50, false, |mid| Ok((mid < 0.6).then_some(true)))
188                .unwrap();
189        assert_eq!((b.lo, b.hi, b.undecided), (0.5, 1.0, Some(0.75)));
190    }
191
192    /// Every sign change between samples is bracketed and bisected, an
193    /// undecided sample passed over.
194    #[test]
195    fn scans_and_bisects_every_sign_change() {
196        let f = |x: f64| (x - 0.3) * (x - 2.7);
197        let side = |x: f64| Ok::<_, GeopError>((x != 1.5).then_some(f(x) > 0.0));
198        let roots: Vec<Bracket<f64>> = sign_changes((0..=6).map(|k| 0.5 * k as f64), side)
199            .bisected(200, true)
200            .collect::<Result<_, _>>()
201            .unwrap();
202        assert_eq!(roots.len(), 2);
203        for (b, below) in roots.iter().zip([0.0, 2.5]) {
204            assert!(b.lo > below && b.hi < below + 0.5);
205            assert_eq!(b.hi, f64::from_bits(b.lo.to_bits() + 1));
206            assert!(f(b.lo) * f(b.hi) <= 0.0);
207        }
208    }
209}
geop_core_math/solvers/bisect.rs

geop_core_math/solvers/
bisect.rs
