lox_core/math/
optim.rs

1// SPDX-FileCopyrightText: 2026 Helge Eichhorn <git@helgeeichhorn.de>
2//
3// SPDX-License-Identifier: MPL-2.0
4
5//! Bracketed optimization algorithms.
6
7use lox_approx::approx_eq;
8use thiserror::Error;
9
10use crate::error::LoxError;
11use crate::math::callback::Callback;
12use crate::math::float::abs;
13
14/// Error returned by bracketed minimization algorithms.
15#[derive(Debug, Error)]
16pub enum MinimizerError {
17    /// The algorithm did not converge within the maximum number of iterations.
18    #[error(
19        "minimization did not converge after {iterations} iterations at x = {x}, f(x) = {value}"
20    )]
21    NotConverged {
22        /// Number of iterations performed before giving up.
23        iterations: u32,
24        /// The best minimizer estimate reached.
25        x: f64,
26        /// The objective value `f(x)` at the best estimate.
27        value: f64,
28    },
29    /// The objective function returned an error.
30    #[error(transparent)]
31    Callback(#[from] LoxError),
32}
33
34/// Finds the minimum of a function within a bracket.
35pub trait FindBracketedMinimum<F>
36where
37    F: Callback,
38{
39    /// Finds the x value that minimizes `f` within the given `bracket`.
40    fn find_minimum_in_bracket(&self, f: F, bracket: (f64, f64)) -> Result<f64, MinimizerError>;
41}
42
43/// Brent's method for finding the minimum of a unimodal function in a bracket.
44///
45/// Combines golden section search with parabolic interpolation.
46#[derive(Debug, Copy, Clone, PartialEq)]
47pub struct BrentMinimizer {
48    /// Maximum number of iterations.
49    pub max_iter: u32,
50    /// Absolute tolerance for convergence.
51    pub abs_tol: f64,
52}
53
54impl Default for BrentMinimizer {
55    fn default() -> Self {
56        Self {
57            max_iter: 500,
58            abs_tol: 1e-10,
59        }
60    }
61}
62
63/// Golden ratio constant used in Brent minimization.
64const GOLDEN: f64 = 0.381_966_011_250_105_1; // (3 - sqrt(5)) / 2
65
66impl<F> FindBracketedMinimum<F> for BrentMinimizer
67where
68    F: Callback,
69{
70    fn find_minimum_in_bracket(&self, f: F, bracket: (f64, f64)) -> Result<f64, MinimizerError> {
71        let (mut a, mut b) = bracket;
72        if a > b {
73            core::mem::swap(&mut a, &mut b);
74        }
75
76        // x is the point with the least function value found so far.
77        // w is the point with the second least value.
78        // v is the previous value of w.
79        let mut x = a + GOLDEN * (b - a);
80        let mut w = x;
81        let mut v = x;
82        let mut fx = f.call(x)?;
83        let mut fw = fx;
84        let mut fv = fx;
85
86        // e is the distance moved on the step before last.
87        // d is the distance moved on the last step.
88        let mut e = 0.0_f64;
89        let mut d = 0.0_f64;
90
91        for _ in 0..self.max_iter {
92            let midpoint = 0.5 * (a + b);
93            let tol1 = self.abs_tol * abs(x) + 1e-10;
94            let tol2 = 2.0 * tol1;
95
96            // Check convergence.
97            if abs(x - midpoint) <= tol2 - 0.5 * (b - a) {
98                return Ok(x);
99            }
100
101            // Try parabolic interpolation.
102            let mut use_golden = true;
103            if abs(e) > tol1 {
104                // Fit parabola through x, v, w.
105                let r = (x - w) * (fx - fv);
106                let q = (x - v) * (fx - fw);
107                let p = (x - v) * q - (x - w) * r;
108                let q = 2.0 * (q - r);
109                let (p, q) = if q > 0.0 { (-p, q) } else { (p, -q) };
110
111                // Is the parabola acceptable?
112                if abs(p) < abs(0.5 * q * e) && p > q * (a - x) && p < q * (b - x) {
113                    e = d;
114                    d = p / q;
115                    let u = x + d;
116
117                    // f must not be evaluated too close to a or b.
118                    if (u - a) < tol2 || (b - u) < tol2 {
119                        d = if x < midpoint { tol1 } else { -tol1 };
120                    }
121                    use_golden = false;
122                }
123            }
124
125            if use_golden {
126                // Golden section step.
127                e = if x < midpoint { b - x } else { a - x };
128                d = GOLDEN * e;
129            }
130
131            // f must not be evaluated too close to x.
132            let u = if abs(d) >= tol1 {
133                x + d
134            } else if d > 0.0 {
135                x + tol1
136            } else {
137                x - tol1
138            };
139
140            let fu = f.call(u)?;
141
142            // Update a, b, v, w, x.
143            if fu <= fx {
144                if u < x {
145                    b = x;
146                } else {
147                    a = x;
148                }
149                v = w;
150                fv = fw;
151                w = x;
152                fw = fx;
153                x = u;
154                fx = fu;
155            } else {
156                if u < x {
157                    a = u;
158                } else {
159                    b = u;
160                }
161                if fu <= fw || approx_eq!(w, x, atol <= 1e-15) {
162                    v = w;
163                    fv = fw;
164                    w = u;
165                    fw = fu;
166                } else if fu <= fv
167                    || approx_eq!(v, x, atol <= 1e-15)
168                    || approx_eq!(v, w, atol <= 1e-15)
169                {
170                    v = u;
171                    fv = fu;
172                }
173            }
174        }
175
176        Err(MinimizerError::NotConverged {
177            iterations: self.max_iter,
178            x,
179            value: fx,
180        })
181    }
182}
183
184#[cfg(test)]
185mod tests {
186    use core::f64::consts::PI;
187    use lox_approx::assert_approx_eq;
188
189    use super::*;
190    use crate::math::float::{cos, powi};
191
192    #[test]
193    fn test_brent_minimizer_quadratic() {
194        let minimizer = BrentMinimizer::default();
195        let x = minimizer
196            .find_minimum_in_bracket(|x: f64| powi(x - 3.0, 2), (0.0, 5.0))
197            .expect("should converge");
198        assert_approx_eq!(x, 3.0, atol <= 1e-8);
199    }
200
201    #[test]
202    fn test_brent_minimizer_cosine() {
203        // cos(x) has a minimum at PI in [PI/2, 3*PI/2]
204        let minimizer = BrentMinimizer::default();
205        let x = minimizer
206            .find_minimum_in_bracket(|x: f64| cos(x), (PI / 2.0, 3.0 * PI / 2.0))
207            .expect("should converge");
208        assert_approx_eq!(x, PI, atol <= 1e-8);
209    }
210
211    #[test]
212    fn test_brent_minimizer_reversed_bracket() {
213        let minimizer = BrentMinimizer::default();
214        let x = minimizer
215            .find_minimum_in_bracket(|x: f64| powi(x - 2.0, 2), (5.0, 0.0))
216            .expect("should converge");
217        assert_approx_eq!(x, 2.0, atol <= 1e-8);
218    }
219
220    #[test]
221    fn test_brent_minimizer_custom_tolerance() {
222        let minimizer = BrentMinimizer {
223            max_iter: 100,
224            abs_tol: 1e-4,
225        };
226        let x = minimizer
227            .find_minimum_in_bracket(|x: f64| powi(x - 1.0, 2), (-2.0, 5.0))
228            .expect("should converge");
229        assert_approx_eq!(x, 1.0, atol <= 1e-3);
230    }
231
232    #[test]
233    fn test_brent_minimizer_not_converged() {
234        let minimizer = BrentMinimizer {
235            max_iter: 0,
236            abs_tol: 1e-15,
237        };
238        let result = minimizer.find_minimum_in_bracket(|x: f64| powi(x - 1.0, 2), (0.0, 5.0));
239        assert!(result.is_err());
240    }
241}
lox_core/math/optim.rs

lox_core/math/
optim.rs
