carapace_core/analysis/
line_search.rs

1/// Post-processes a raw Newton correction to actually reduce the residual,
2/// rather than accepting it outright. Closed enum. Ported from
3/// Xara's `SRC/analysis/algorithm/equiSolnAlgo/search/` (`docs/
4/// algorithms.md` §4): both variants bracket a sign change in the 1-D
5/// "energy" `s(eta) = du . residual(eta)` by geometric expansion of `eta`,
6/// then refine within the bracket — bisection, or false-position
7/// (`RegulaFalsi`, usually fewer iterations than bisection once bracketed).
8/// Only these two for now — `Secant`/`InitialInterpolatedLineSearch` can
9/// diverge without a bracket and aren't needed until a specific model
10/// demands the extra speed (`docs/algorithms.md` §4).
11#[derive(Debug, Clone, Copy, PartialEq)]
12pub enum LineSearch {
13    Bisection {
14        tol: f64,
15        max_iter: usize,
16        max_eta: f64,
17    },
18    RegulaFalsi {
19        tol: f64,
20        max_iter: usize,
21        max_eta: f64,
22    },
23}
24
25impl LineSearch {
26    fn params(&self) -> (f64, usize, f64) {
27        match *self {
28            LineSearch::Bisection {
29                tol,
30                max_iter,
31                max_eta,
32            }
33            | LineSearch::RegulaFalsi {
34                tol,
35                max_iter,
36                max_eta,
37            } => (tol, max_iter, max_eta),
38        }
39    }
40
41    /// Given `s0 = du . residual` at `eta = 0` and `s1 = du . residual` at
42    /// `eta = 1` (the raw, un-rescaled correction), and a closure to
43    /// evaluate `s(eta)` at any other trial `eta`, returns the `eta` this
44    /// strategy converges to. `s0` and `s1` must have opposite signs
45    /// (checked by the caller — `Analysis::try_step` only calls this once
46    /// the ordinary convergence test at `eta = 1` has already failed, per
47    /// `docs/algorithms.md` §4 step 3; a same-sign `s0`/`s1` means the raw
48    /// step already reduced the energy monotonically, so bracketing
49    /// expansion below is what finds the far side of the sign change).
50    pub(crate) fn resolve(&self, s0: f64, s1: f64, mut eval: impl FnMut(f64) -> f64) -> f64 {
51        let (tol, max_iter, max_eta) = self.params();
52
53        // Bracket a sign change by geometric expansion, starting from the
54        // raw step's own [0, 1] interval (Xara's `BisectionLineSearch`/
55        // `RegulaFalsiLineSearch` both start this way).
56        let (mut eta_lo, mut s_lo) = (0.0, s0);
57        let (mut eta_hi, mut s_hi) = (1.0, s1);
58        let mut expand_iter = 0;
59        while s_lo.signum() == s_hi.signum() && eta_hi < max_eta && expand_iter < max_iter {
60            eta_lo = eta_hi;
61            s_lo = s_hi;
62            eta_hi = (eta_hi * 4.0).min(max_eta);
63            s_hi = eval(eta_hi);
64            expand_iter += 1;
65        }
66        if s_lo.signum() == s_hi.signum() {
67            // Never bracketed (e.g. `s` monotonically one-signed out to
68            // `max_eta`) — the best available point is whichever endpoint
69            // has the smaller `|s|`, not a fabricated root.
70            return if s_lo.abs() < s_hi.abs() {
71                eta_lo
72            } else {
73                eta_hi
74            };
75        }
76
77        for _ in 0..max_iter {
78            let eta_trial = match self {
79                LineSearch::Bisection { .. } => 0.5 * (eta_lo + eta_hi),
80                // False position: linear-interpolate the root between the
81                // bracket's two endpoints instead of always bisecting.
82                LineSearch::RegulaFalsi { .. } => {
83                    eta_lo + (eta_hi - eta_lo) * (-s_lo) / (s_hi - s_lo)
84                }
85            };
86            let s_trial = eval(eta_trial);
87            if s_trial.abs() < tol {
88                return eta_trial;
89            }
90            if s_trial.signum() == s_lo.signum() {
91                eta_lo = eta_trial;
92                s_lo = s_trial;
93            } else {
94                eta_hi = eta_trial;
95                s_hi = s_trial;
96            }
97        }
98        0.5 * (eta_lo + eta_hi)
99    }
100}
carapace_core/analysis/line_search.rs

carapace_core/analysis/
line_search.rs
