sim_lib_numbers_optimize/
scalar.rs

1//! Bounded scalar minimization and finite-difference gradients.
2
3use super::*;
4
5/// Bounded Brent minimization. This proves local bracket convergence only.
6pub fn minimize_scalar<F>(
7    mut f: F,
8    mut a: f64,
9    mut b: f64,
10    tol: f64,
11    limits: Limits,
12) -> Result<ScalarResult, Error>
13where
14    F: FnMut(f64) -> f64,
15{
16    if !a.is_finite() || !b.is_finite() || a >= b || !tol.is_finite() || tol <= 0.0 {
17        return Err(Error::InvalidPlan(
18            "scalar interval and tolerance must be finite and ordered",
19        ));
20    }
21    let golden = 0.3819660112501051;
22    let mut x = a + golden * (b - a);
23    let (mut w, mut v) = (x, x);
24    let mut fx = f(x);
25    let (mut fw, mut fv) = (fx, fx);
26    let (mut d, mut e) = (0.0_f64, 0.0_f64);
27    let mut evals = 1;
28    if !fx.is_finite() {
29        return Ok(ScalarResult {
30            minimizer: x,
31            value: fx,
32            final_bracket: (a, b),
33            termination: Termination::NonFinite,
34            work: Work {
35                evaluations: evals,
36                iterations: 0,
37                memory_bytes: 0,
38            },
39        });
40    }
41    for iter in 0..limits.iterations {
42        let m = 0.5 * (a + b);
43        let t = tol * x.abs() + f64::EPSILON.sqrt();
44        if (x - m).abs() <= 2.0 * t - 0.5 * (b - a) {
45            return Ok(ScalarResult {
46                minimizer: x,
47                value: fx,
48                final_bracket: (a, b),
49                termination: Termination::Converged,
50                work: Work {
51                    evaluations: evals,
52                    iterations: iter,
53                    memory_bytes: 0,
54                },
55            });
56        }
57        let old = e;
58        e = d;
59        if old.abs() > t {
60            let r = (x - w) * (fx - fv);
61            let mut q = (x - v) * (fx - fw);
62            let mut p = (x - v) * q - (x - w) * r;
63            q = 2.0 * (q - r);
64            if q > 0.0 {
65                p = -p
66            } else {
67                q = -q
68            };
69            if p.abs() >= 0.5 * q * old.abs() || p <= q * (a - x) || p >= q * (b - x) {
70                e = if x < m { b - x } else { a - x };
71                d = golden * e
72            } else {
73                d = p / q;
74            }
75        } else {
76            e = if x < m { b - x } else { a - x };
77            d = golden * e
78        }
79        let u = x + if d.abs() >= t { d } else { t.copysign(d) };
80        if evals >= limits.evaluations {
81            return Ok(ScalarResult {
82                minimizer: x,
83                value: fx,
84                final_bracket: (a, b),
85                termination: Termination::WorkLimit,
86                work: Work {
87                    evaluations: evals,
88                    iterations: iter,
89                    memory_bytes: 0,
90                },
91            });
92        }
93        let fu = f(u);
94        evals += 1;
95        if !fu.is_finite() {
96            return Ok(ScalarResult {
97                minimizer: x,
98                value: fx,
99                final_bracket: (a, b),
100                termination: Termination::NonFinite,
101                work: Work {
102                    evaluations: evals,
103                    iterations: iter,
104                    memory_bytes: 0,
105                },
106            });
107        }
108        if fu <= fx {
109            if u < x {
110                b = x
111            } else {
112                a = x
113            };
114            v = w;
115            fv = fw;
116            w = x;
117            fw = fx;
118            x = u;
119            fx = fu
120        } else {
121            if u < x {
122                a = u
123            } else {
124                b = u
125            };
126            if fu <= fw || w == x {
127                v = w;
128                fv = fw;
129                w = u;
130                fw = fu
131            } else if fu <= fv || v == x || v == w {
132                v = u;
133                fv = fu
134            }
135        }
136    }
137    Ok(ScalarResult {
138        minimizer: x,
139        value: fx,
140        final_bracket: (a, b),
141        termination: Termination::WorkLimit,
142        work: Work {
143            evaluations: evals,
144            iterations: limits.iterations,
145            memory_bytes: 0,
146        },
147    })
148}
149
150pub(crate) fn numerical_gradient<F: FnMut(&[f64]) -> f64>(
151    f: &mut F,
152    x: &[f64],
153    fx: f64,
154    scale: &[f64],
155    evals: &mut usize,
156    limit: usize,
157) -> Option<Vec<f64>> {
158    let mut g = vec![0.0; x.len()];
159    for i in 0..x.len() {
160        if *evals >= limit {
161            return None;
162        }
163        let mut y = x.to_vec();
164        let h = f64::EPSILON.sqrt() * (x[i].abs() + scale[i]);
165        y[i] += h;
166        let fy = f(&y);
167        *evals += 1;
168        if !fy.is_finite() {
169            return None;
170        }
171        g[i] = (fy - fx) / h;
172    }
173    Some(g)
174}
sim_lib_numbers_optimize/scalar.rs

sim_lib_numbers_optimize/
scalar.rs
