scivex-optim 0.1.0

Scivex — Optimization, root finding, and numerical integration
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

//! Golden-section search for 1-D minimization.

use scivex_core::Float;

use crate::error::Result;

use super::Minimize1dResult;

/// Find the minimum of a unimodal function `f` on `[a, b]` using
/// golden-section search.
///
/// Linear convergence, reducing the bracket by a factor of the golden
/// ratio (~0.618) each iteration.
///
/// # Examples
///
/// ```
/// # use scivex_optim::minimize_1d::golden_section;
/// // Minimum of (x - 3)^2 on [0, 10]
/// let result = golden_section(|x: f64| (x - 3.0) * (x - 3.0), 0.0, 10.0, 1e-8, 100).unwrap();
/// assert!((result.x_min - 3.0).abs() < 1e-6);
/// ```
pub fn golden_section<T: Float, F: Fn(T) -> T>(
    f: F,
    a: T,
    b: T,
    xtol: T,
    max_iter: usize,
) -> Result<Minimize1dResult<T>> {
    let gr = T::from_f64(0.618_033_988_749_895); // (sqrt(5) - 1) / 2
    let one_minus_gr = T::one() - gr;

    let mut a = a;
    let mut b = b;

    let mut x1 = a + one_minus_gr * (b - a);
    let mut x2 = a + gr * (b - a);
    let mut f1 = f(x1);
    let mut f2 = f(x2);

    for i in 0..max_iter {
        if (b - a).abs() < xtol {
            let (x_min, f_min) = if f1 < f2 { (x1, f1) } else { (x2, f2) };
            return Ok(Minimize1dResult {
                x_min,
                f_min,
                iterations: i + 1,
                converged: true,
            });
        }

        if f1 < f2 {
            b = x2;
            x2 = x1;
            f2 = f1;
            x1 = a + one_minus_gr * (b - a);
            f1 = f(x1);
        } else {
            a = x1;
            x1 = x2;
            f1 = f2;
            x2 = a + gr * (b - a);
            f2 = f(x2);
        }
    }

    let (x_min, f_min) = if f1 < f2 { (x1, f1) } else { (x2, f2) };
    Ok(Minimize1dResult {
        x_min,
        f_min,
        iterations: max_iter,
        converged: false,
    })
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_golden_section_quadratic() {
        // (x - 3)^2 on [0, 5] => minimum at x = 3
        let result = golden_section(|x: f64| (x - 3.0) * (x - 3.0), 0.0, 5.0, 1e-10, 200).unwrap();
        assert!(result.converged);
        assert!((result.x_min - 3.0).abs() < 1e-8);
        assert!(result.f_min < 1e-14);
    }

    #[test]
    fn test_golden_section_cubic_min() {
        // x^4 - 2x^2 has local minima at ±1. On [0, 3] the minimum is at 1.
        let result =
            golden_section(|x: f64| x.powi(4) - 2.0 * x * x, 0.0, 3.0, 1e-10, 200).unwrap();
        assert!(result.converged);
        assert!((result.x_min - 1.0).abs() < 1e-6);
    }
}