geop_core_math/solvers/system/
driver.rs

1//! The driver every solve of a [`System`] shares: one group solved
2//! exactly, and among where its constraints hold preferring what a caller
3//! pulls towards and, after that, where things were — a sketch's numbers
4//! and a program's state alike.
5//!
6//! A solve is one minimization by constrained Levenberg–Marquardt
7//! ([`least_squares::minimize`]): the constraints held exactly, the pulls
8//! and the damping in the sum (see [`DAMPING`]). Where the constraints
9//! cannot all hold — or the minimization could not tell how to meet them —
10//! a second one takes them into the sum too, as near as they come, least
11//! squares, still damped so that what they leave free stays put. Each
12//! minimization is reported as a [`Phase`]. Each block moves by steps along
13//! it ([`Variable::retract`]).
14
15use crate::{
16    geop_error::GeopResult,
17    scalars::Scalar,
18    solvers::least_squares::{self, Evaluation, Rows, Space, Stop},
19};
20
21use super::{Bound, Equations, Solving, Term, Variable};
22
23/// A residual holds once it is within this fraction of the system's size:
24/// one part in `RELATIVE_TOLERANCE`.
25pub const RELATIVE_TOLERANCE: i64 = 1_000_000_000;
26
27/// How near zero a residual of a system of size `scale` must be to hold:
28/// one part in [`RELATIVE_TOLERANCE`] of it.
29pub fn tolerance(scale: Scalar) -> GeopResult<Scalar> {
30    scale.div(Scalar::from_i64(RELATIVE_TOLERANCE))
31}
32
33/// How strongly every value a solve may move and nothing pulls is held
34/// where it was, as a ratio to a pull, which weighs one: the damping that
35/// makes a solve prefer, among the many configurations the constraints
36/// allow, the one nearest where things were. A value the constraints leave
37/// free stays put; one they need moves as far as they need, and no further.
38///
39/// Ranked ten thousand to one below a pull: far enough apart that the
40/// damping yields to a pull (to a hundred-millionth, as they enter
41/// squared), and no further — the minimizer resolves the weaker only as
42/// well as the stronger lets it, and weights ten orders of magnitude apart
43/// in the sum leave the weaker unresolved.
44pub const DAMPING: (i64, i64) = (1, 10_000);
45
46/// The most steps one minimization takes: how hard a solve tries, never
47/// what its answer means.
48const MAX_ITERATIONS: usize = 200;
49
50/// What a solve may do with a value.
51#[derive(Clone, Copy, Debug, PartialEq, Eq)]
52pub enum Mobility {
53    /// Nothing: it stays exactly as it is.
54    Fixed,
55    /// Move it as far as the constraints need, and no further: it is held
56    /// where it was, by [`DAMPING`].
57    Held,
58    /// Move it as the constraints have it, and hold it nowhere: it follows
59    /// the others — a joint's coordinate, measured from its parts. Held as
60    /// well, it would pull against them.
61    Follows,
62}
63
64impl Mobility {
65    /// Whether a solve may move it.
66    pub fn free(self) -> bool {
67        self != Mobility::Fixed
68    }
69}
70
71/// One minimization of a solve: why it stopped, after how many steps, and
72/// the largest residual it left — an upper bound.
73#[derive(Clone, Copy, Debug, PartialEq, serde::Serialize, serde::Deserialize)]
74pub struct Phase {
75    pub stop: Stop,
76    pub iterations: usize,
77    /// Steps proposed and turned down on the way (see
78    /// [`least_squares::Outcome::turned_down`]).
79    pub turned_down: usize,
80    #[serde(with = "crate::scalars::as_f64")]
81    pub max_residual: Scalar,
82}
83
84/// The largest size of `residuals`: an upper bound.
85pub(super) fn largest(residuals: impl IntoIterator<Item = Scalar>) -> Scalar {
86    residuals.into_iter().fold(Scalar::ZERO, |largest, r| {
87        let size = r.abs().upper();
88        if size.could_be_greater(largest) {
89            size
90        } else {
91            largest
92        }
93    })
94}
95
96/// Where a solve of a group ended, and the phases that put it there.
97pub(super) type Phased<E> = (Vec<<E as Equations>::Value>, Vec<Phase>);
98
99/// Where one minimization ended, how, and whether the constraints hold
100/// there.
101type Minimized<E> = (Vec<<E as Equations>::Value>, Phase, bool);
102
103/// The values of a group of blocks, as the points a minimization moves:
104/// each stepped along its own steps ([`Variable::retract`]).
105pub(super) struct Values<V>(pub(super) std::marker::PhantomData<V>);
106
107impl<V: Variable> Space for Values<V> {
108    type Point = Vec<V>;
109
110    fn dimension(&self, at: &Vec<V>) -> usize {
111        at.iter().map(Variable::dimension).sum()
112    }
113
114    fn step(&self, at: &Vec<V>, step: &[Scalar]) -> Option<Vec<V>> {
115        let mut rest = step;
116        let mut next = Vec::with_capacity(at.len());
117        for value in at {
118            let (own, after) = rest.split_at_checked(value.dimension())?;
119            next.push(value.retract(own).ok()?);
120            rest = after;
121        }
122        Some(next)
123    }
124
125    fn same(&self, a: &Vec<V>, b: &Vec<V>) -> bool {
126        a.iter().zip(b).all(|(a, b)| a.same(b))
127    }
128}
129
130impl<E: Equations> Solving<'_, '_, E> {
131    /// The group solved from `start` with the limits `active` held at
132    /// their bounds and the constraints `left_out` left out (see
133    /// [`Asked::solve`]).
134    pub(super) async fn solve(
135        &self,
136        start: Vec<E::Value>,
137        active: &[usize],
138        left_out: &[usize],
139    ) -> GeopResult<Phased<E>> {
140        let hard: Vec<usize> = (self.group.constraints.iter())
141            .filter(|c| !left_out.contains(c))
142            .filter(|&&c| self.system.equations.bound(c) == Bound::Zero || active.contains(&c))
143            .copied()
144            .collect();
145        Asked {
146            solving: self,
147            hard,
148        }
149        .solve(start)
150        .await
151    }
152}
153
154/// A group as one solve has it: the constraints `hard` to hold.
155struct Asked<'a, 's, 'e, E: Equations> {
156    solving: &'a Solving<'s, 'e, E>,
157    hard: Vec<usize>,
158}
159
160/// The residuals of a group at a point, each a length, and their Jacobians
161/// by steps from it: the constraints' — zero where they hold — the pulls'
162/// — each the distance from where a value is pulled to — and those of every
163/// value [`Mobility::Held`] and not pulled, from where it was.
164struct Terms {
165    constraints: Rows,
166    pulls: Rows,
167    stays: Rows,
168}
169
170impl<E: Equations> Asked<'_, '_, '_, E> {
171    /// The residuals at `values` — the pulls only if `pulling` (else
172    /// nothing is pulled, and every value held stays), every value held
173    /// staying where it is at `anchor`. `None` where one is undecidable,
174    /// which makes `values` infeasible.
175    async fn terms(
176        &self,
177        values: &[E::Value],
178        anchor: &[E::Value],
179        pulling: bool,
180    ) -> Option<Terms> {
181        let solving = self.solving;
182        let identity = |_: usize| None;
183        let mut which: Vec<Term> = self.hard.iter().map(|&c| Term::Constraint(c)).collect();
184        if pulling {
185            which.extend(
186                solving
187                    .group
188                    .preferences
189                    .iter()
190                    .map(|&p| Term::Preference(p)),
191            );
192        }
193        let made = solving.made(values, &which).await;
194        let mut constraints = Rows::default();
195        let mut pulls = Rows::default();
196        let mut pulled = Vec::new();
197        for (term, made) in which.iter().zip(made) {
198            let made = made?;
199            match term {
200                Term::Constraint(_) => constraints.append(solving.rows(&made, &identity)),
201                Term::Preference(_) => {
202                    pulled.extend(made.jacobian.iter().map(|(block, _)| *block));
203                    pulls.append(solving.rows(&made, &identity));
204                }
205            }
206        }
207        let mut stays = Rows::default();
208        for (k, ((&block, value), at)) in (solving.group.blocks.iter())
209            .zip(values)
210            .zip(anchor)
211            .enumerate()
212        {
213            if pulled.contains(&block) || solving.system.mobility[block] != Mobility::Held {
214                continue;
215            }
216            let stay = value.stay(at, solving.system.scale).ok()?;
217            let column = solving.columns[k];
218            for (value, own) in stay.values.into_iter().zip(stay.jacobian) {
219                let mut row = vec![Scalar::ZERO; solving.n];
220                for (slot, v) in row.iter_mut().skip(column).zip(own) {
221                    *slot = v;
222                }
223                stays.push(value, row);
224            }
225        }
226        Some(Terms {
227            constraints,
228            pulls,
229            stays,
230        })
231    }
232
233    /// The largest residual of the constraints at `values` — an upper
234    /// bound, infinite where one is undecidable (see [`largest`]) — and
235    /// whether they all hold there.
236    async fn check(&self, values: &[E::Value]) -> GeopResult<(Scalar, bool)> {
237        let tolerance = tolerance(self.solving.system.scale)?;
238        Ok(self.solving.check(values, &self.hard, tolerance).await)
239    }
240
241    /// Where the group has to be, from `start`, for its constraints to
242    /// hold, as near where its pulls want it and, among those, as near
243    /// `start` as the constraints allow — or, where they cannot all hold,
244    /// where they come closest (see the module docs). The point is given
245    /// whether or not they hold, with the phases that put it there.
246    async fn solve(&self, start: Vec<E::Value>) -> GeopResult<Phased<E>> {
247        let (solved, phase, holds) = self.minimize(start, self.solving.start(), true).await?;
248        let mut phases = vec![phase];
249        if holds {
250            return Ok((solved, phases));
251        }
252        let (solved, phase, _) = self.minimize(solved.clone(), solved, false).await?;
253        phases.push(phase);
254        Ok((solved, phases))
255    }
256
257    /// One minimization from `from`: the constraints held exactly and the
258    /// group pulled if `hard`, else the constraints in the sum; every value
259    /// held where it is at `anchor` — where the group was, for a pendant
260    /// started where what it hangs off carried it. Where it ended, how, and
261    /// whether the constraints hold there.
262    async fn minimize(
263        &self,
264        from: Vec<E::Value>,
265        anchor: Vec<E::Value>,
266        hard: bool,
267    ) -> GeopResult<Minimized<E>> {
268        let space = &self.solving.space;
269        if space.dimension(&from) == 0 {
270            let (max_residual, holds) = self.check(&from).await?;
271            let phase = Phase {
272                stop: Stop::NoChange { rejected: None },
273                iterations: 0,
274                turned_down: 0,
275                max_residual,
276            };
277            return Ok((from, phase, holds));
278        }
279        let damping = Scalar::from_ratio(DAMPING.0, DAMPING.1)?;
280        let evaluate = async |at: &Vec<E::Value>| {
281            let Terms {
282                constraints,
283                pulls,
284                stays,
285            } = self.terms(at, &anchor, hard).await?;
286            let mut evaluation = Evaluation::default();
287            if hard {
288                evaluation.constraints = constraints;
289            } else {
290                evaluation.sum = constraints;
291            }
292            let sum = &mut evaluation.sum;
293            sum.append(pulls);
294            for (value, row) in stays.values.into_iter().zip(stays.jacobian) {
295                sum.push(
296                    value.mul(damping),
297                    row.into_iter().map(|d| d.mul(damping)).collect(),
298                );
299            }
300            // An unbounded residual or slope is no number either.
301            let finite = |rows: &Rows| {
302                (rows.values.iter())
303                    .chain(rows.jacobian.iter().flatten())
304                    .all(|v| v.is_finite())
305            };
306            (finite(&evaluation.constraints) && finite(&evaluation.sum)).then_some(evaluation)
307        };
308        let outcome = least_squares::minimize(
309            space,
310            evaluate,
311            from,
312            least_squares::Options {
313                max_iterations: MAX_ITERATIONS,
314                max_step: self.solving.system.scale,
315            },
316        )
317        .await?;
318        let (max_residual, holds) = self.check(&outcome.x).await?;
319        let phase = Phase {
320            stop: outcome.stop,
321            iterations: outcome.iterations,
322            turned_down: outcome.turned_down,
323            max_residual,
324        };
325        Ok((outcome.x, phase, holds))
326    }
327}
geop_core_math/solvers/system/driver.rs

geop_core_math/solvers/system/
driver.rs
