rustamath_mnmz/
golden_section.rs

1//! Golden section search for a minimum.
2//!
3//! (c) 2023 Igor Lesik
4//! MIT license
5//!
6//! References:
7//!
8//! 1. William H. Press - Numerical recipes, the art of scientific computing.
9//!   Cambridge University Press (2007).
10//!
11use super::bracket::{find_bracket, shft3, shft2};
12
13/// Smallest tolerance.
14///
15/// See book "Numerical recipes, the art of scientific computing."
16/// sqrt(f64 precision 10^16), by Tailor series for `f(x+eps)`
17const MIN_TOLERANCE: f64 = 3.0e-8_f64;
18
19/// Golden section search for a minimum.
20///
21/// - William H. Press - Numerical recipes, the art of scientific computing.
22///   Cambridge University Press (2007).
23///
24/// Given a function `f`, and given a bracketing triplet of abscissas a, b, c
25/// (such that b is between a and c, and f(x) is less than both f(a) and f(c)),
26/// this routine performs a golden section search for the minimum,
27/// isolating it to a fractional precision of about `tolerance`.
28///
29/// # Example
30///
31/// ```
32/// use rustamath_mnmz::golden_section_search;
33/// use assert_float_eq::*;
34/// // Roots 1.0 and 2.0, minimum at 1.5.
35/// let poly2 = |x: f64| (x-1.0)*(x-2.0);
36/// let ranges = vec![(10.0, 20.0), (20.0, 10.0), (-10.0, 0.0)];
37/// for range in ranges {
38///     let (xmin, f, nr_iterations) = golden_section_search(poly2, range.0, range.1, 0.0, 0);
39///     println!("MIN: {:.8} f(xmin): {:6.2} iterations:{}",
40///         xmin, f, nr_iterations
41///     );
42///     assert_float_relative_eq!(xmin, 1.5, 1.0e-8);
43/// }
44///
45pub fn golden_section_search<F: Fn (f64) -> f64>(
46    fun: F,
47    a: f64,
48    b: f64,
49    tol: f64,
50    max_iterations: usize
51) -> (f64, f64, usize)
52{
53    let tol = tol.max(MIN_TOLERANCE);
54    let max_iterations = if max_iterations < 1 { 500 } else { max_iterations.min(1000) };
55    const R: f64 = 0.61803399_f64;
56    const C: f64 = 1.0 - R; // The golden ratios.
57
58    let bracket = find_bracket(&fun, a, b);
59    let a = bracket.a;
60    let b = bracket.b;
61    let c = bracket.c;
62
63    // At any given time we will keep track of four points, x0,x1,x2,x3.
64    let mut x1: f64;
65    let mut x2: f64;
66    let mut x0 = a;
67    let mut x3 = c;
68
69    // Make x0 to x1 the smaller segment, and fill in the new point to be tried.
70    if (c-b).abs() > (b-a).abs() {
71        x1 = b;
72        x2 = b + C*(c-b);
73    } else {
74        x2 = b;
75        x1 = b - C*(b-a);
76    }
77
78    // The initial function evaluations. Note that we never need to evaluate
79    // the function at the original endpoints.
80    let mut f1 = fun(x1);
81    let mut f2 = fun(x2);
82    let mut nr_iterations: usize = 0;
83
84    while (x3-x0).abs() > tol*(x1.abs() + x2.abs()) {
85        if f2 < f1 {
86            let d = R*x2 + C*x3;
87            shft3(&mut x0, &mut x1, &mut x2, d);
88            shft2(&mut f1, &mut f2, fun(x2));
89        }
90        else {
91            let d = R*x1 + C*x0;
92            shft3(&mut x3, &mut x2, &mut x1, d);
93            shft2(&mut f2, &mut f1, fun(x1));
94        }
95        nr_iterations += 1;
96
97        if nr_iterations >= max_iterations { break; }
98
99        // @igor: saw/non-smooth functions demostrate that `tol*(x1.abs() + x2.abs())`
100        // gets smaller faster than `(x3-x0).abs()` preventing the conversion;
101        // here if we see that x3 is close to x0 and f1 to f2 we force the exit.
102        if nr_iterations > 10 && (x3-x0).abs() < tol && (f1 - f2).abs() < tol { break; }
103    }
104
105    // Output the best of the two current values.
106    if f1 < f2 {
107        (x1, f1, nr_iterations)
108    }
109    else {
110        (x2, f2, nr_iterations)
111    }
112}
113
114#[cfg(test)]
115#[test]
116fn test_poly2() {
117    // Roots 1.0 and 2.0, minimum at 1.5.
118    let poly2 = |x: f64| (x-1.0)*(x-2.0);
119
120    let ranges = vec![(10.0, 20.0), (20.0, 10.0), (-10.0, 0.0),
121        (-2000.0, -1000.0), (-10_000.0, 30_000.0), (0.0001, 0.0002), (-0.00001, 1.4999)];
122
123    for range in ranges {
124        let (xmin, f, nr_iterations) = golden_section_search(poly2, range.0, range.1, 0.0, 0);
125
126        println!("MIN: {:.8} f(xmin): {:6.2} iterations:{}",
127            xmin, f, nr_iterations
128        );
129
130        assert_float_relative_eq!(xmin, 1.5, 1.0e-8);
131    }
132}
133
134#[cfg(test)]
135#[test]
136fn test_cosine() {
137    // Minimum at Pi on [0, 2*Pi].
138    let cosine = |x: f64| x.cos();
139
140    let ranges = vec![(0.01, 1.0)];
141
142    for range in ranges {
143        let (xmin, f, nr_iterations) = golden_section_search(cosine, range.0, range.1, 0.0, 0);
144
145        println!("MIN: {:.8} f(xmin): {:6.2} iterations:{}",
146            xmin, f, nr_iterations
147        );
148
149        assert_float_relative_eq!(xmin, std::f64::consts::PI, 1.0e-8);
150    }
151}
152
153#[cfg(test)]
154#[test]
155fn test_saw() {
156    let saw = |x: f64| if  x >= 0.0 { x*x*x } else { -x / 1000.0 } ;
157
158    let ranges = vec![(10.0, 20.0), (20.0, 10.0), (-10.0, 0.0),
159        (-2000.0, -1000.0), (-10_000.0, 30_000.0), (0.0001, 0.0002), (-0.00001, 1.4999)];
160
161        for range in ranges {
162            let (xmin, f, nr_iterations) = golden_section_search(saw, range.0, range.1, 1.0e-5, 0);
163
164            println!("MIN: {:.8} f(xmin): {:6.2} iterations:{}",
165                xmin, f, nr_iterations
166            );
167
168            assert_float_absolute_eq!(xmin, 0.0, 1.0e-5);
169        }
170}
rustamath_mnmz/golden_section.rs

rustamath_mnmz/
golden_section.rs
