geop_core_math/solvers/
subdivide.rs

1//! Subdivision search: the global method that finds *every* solution in a
2//! region by discarding the pieces that provably hold none and splitting
3//! the rest (see "Subdivide to isolate, then Newton to refine" in
4//! `AGENTS.md`).
5//!
6//! Which pieces hold no solution, when a piece is small enough to report,
7//! and how to split one are the caller's, in `visit`; this is the search
8//! around them. It promises completeness: every piece `visit` keeps is
9//! visited, in a fixed order, unless the node budget runs out — and then
10//! it says so, with how many pieces were left, rather than returning as if
11//! the search were done. A search cut short is incomplete, and neither
12//! "nothing found" nor anything else may be read into it.
13
14use std::collections::VecDeque;
15
16/// What a visit made of its piece.
17pub enum Split<T> {
18    /// Nothing more to search in it: it holds no solution, or the visit
19    /// has recorded the one it holds.
20    Done,
21    /// The piece narrowed to where its solutions can be (a clip).
22    Restrict(T),
23    /// The piece cut in two, both to be searched.
24    Halves(T, T),
25}
26
27/// In which order pieces are visited.
28#[derive(Clone, Copy, Debug, PartialEq, Eq)]
29pub enum Order {
30    /// Level by level, each piece's halves after every piece already
31    /// waiting: what keeps a budget spent evenly over the region.
32    BreadthFirst,
33    /// Each piece's first half, and all of it, before its second.
34    DepthFirst,
35}
36
37/// How a search ended.
38#[derive(Clone, Copy, Debug, PartialEq, Eq)]
39#[must_use]
40pub enum Searched {
41    /// Every piece was visited.
42    Complete,
43    /// The budget ran out with `pending` pieces left unvisited: the search
44    /// is incomplete.
45    Exhausted { pending: usize },
46}
47
48/// Searches `roots` — in their order — and every piece a visit restricts
49/// or splits them into, visiting at most `max_nodes` pieces. `visit` gets a
50/// piece and the budget left after it (for a search of its own) and says
51/// what to do with it.
52pub fn subdivide<T, E>(
53    roots: impl IntoIterator<Item = T>,
54    order: Order,
55    max_nodes: usize,
56    mut visit: impl FnMut(T, usize) -> Result<Split<T>, E>,
57) -> Result<Searched, E> {
58    let mut pending: VecDeque<T> = VecDeque::new();
59    match order {
60        Order::BreadthFirst => pending.extend(roots),
61        Order::DepthFirst => roots.into_iter().for_each(|r| pending.push_front(r)),
62    }
63    let mut explored = 0usize;
64    // Next at the front, either way: depth first keeps a stack there.
65    while let Some(piece) = pending.pop_front() {
66        if explored >= max_nodes {
67            return Ok(Searched::Exhausted {
68                pending: pending.len() + 1,
69            });
70        }
71        explored += 1;
72        match visit(piece, max_nodes - explored)? {
73            Split::Done => {}
74            Split::Restrict(restricted) => match order {
75                Order::BreadthFirst => pending.push_back(restricted),
76                Order::DepthFirst => pending.push_front(restricted),
77            },
78            Split::Halves(first, second) => match order {
79                Order::BreadthFirst => pending.extend([first, second]),
80                Order::DepthFirst => {
81                    pending.push_front(second);
82                    pending.push_front(first);
83                }
84            },
85        }
86    }
87    Ok(Searched::Complete)
88}
89
90#[cfg(test)]
91mod tests {
92    use super::*;
93
94    /// The pieces of `[0, 1)` that a bisection to `1/16` keeps around
95    /// `0.3` and `0.7`, in either order, and how many it visits.
96    fn search(order: Order, max_nodes: usize) -> (Searched, Vec<(f64, f64)>) {
97        let mut found = Vec::new();
98        let searched = subdivide::<_, ()>([(0.0, 1.0)], order, max_nodes, |(lo, hi), _| {
99            if ![0.3, 0.7].iter().any(|&r| lo <= r && r < hi) {
100                return Ok(Split::Done);
101            }
102            if hi - lo <= 1.0 / 16.0 {
103                found.push((lo, hi));
104                return Ok(Split::Done);
105            }
106            let mid = 0.5 * (lo + hi);
107            Ok(Split::Halves((lo, mid), (mid, hi)))
108        })
109        .unwrap();
110        (searched, found)
111    }
112
113    #[test]
114    fn every_solution_is_found_in_either_order() {
115        let both = [(0.25, 0.3125), (0.6875, 0.75)];
116        for order in [Order::BreadthFirst, Order::DepthFirst] {
117            let (searched, found) = search(order, 1000);
118            assert_eq!(searched, Searched::Complete);
119            assert_eq!(found, both);
120        }
121    }
122
123    /// Breadth first visits by level; depth first finishes the left half
124    /// before the right.
125    #[test]
126    fn the_orders_differ_only_in_order() {
127        let mut seen = [Vec::new(), Vec::new()];
128        for (k, order) in [Order::BreadthFirst, Order::DepthFirst]
129            .into_iter()
130            .enumerate()
131        {
132            let _ = subdivide::<_, ()>([(0u32, 4u32)], order, 100, |(lo, hi), _| {
133                seen[k].push((lo, hi));
134                Ok(if hi - lo == 1 {
135                    Split::Done
136                } else {
137                    let mid = (lo + hi) / 2;
138                    Split::Halves((lo, mid), (mid, hi))
139                })
140            })
141            .unwrap();
142        }
143        assert_eq!(
144            seen[0],
145            [(0, 4), (0, 2), (2, 4), (0, 1), (1, 2), (2, 3), (3, 4)]
146        );
147        assert_eq!(
148            seen[1],
149            [(0, 4), (0, 2), (0, 1), (1, 2), (2, 4), (2, 3), (3, 4)]
150        );
151    }
152
153    /// A budget too small for the search reports it incomplete, never a
154    /// result.
155    #[test]
156    fn running_out_is_reported() {
157        let (searched, _) = search(Order::BreadthFirst, 5);
158        assert_eq!(searched, Searched::Exhausted { pending: 4 });
159    }
160}
geop_core_math/solvers/subdivide.rs

geop_core_math/solvers/
subdivide.rs
