scialg 0.3.0

A collection of scientific algorithms
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
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

//! Minimization of functions

use num::Zero;

use crate::vector::Vector;

/// Find a local minimum of *f* in *[a, c]* using the Golden-section search
///
/// # Example
/// ```
/// use scialg::optimize::golden_section;
///
/// let min = golden_section(f32::sin, 3.0, 6.0, 1e-4);
///
/// assert!((min - 1.5 * std::f32::consts::PI).abs() < 5e-4);
/// ```
///
/// # References
///  - [Wikipedia: Golden-section search](https://en.wikipedia.org/wiki/Golden-section_search)
pub fn golden_section(f: fn(f32) -> f32, a: f32, c: f32, tol: f32) -> f32 {
    let phi = 1.618033988749894848204586834365638118_f32;

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

    while (b - a).abs() > tol {
        let c = b - (b - a) / phi;
        let d = a + (b - a) / phi;
        if f(c) < f(d) {
            b = d;
        } else {
            a = c;
        }
    }

    (b + a) / 2.0
}

fn co_sort<O: PartialOrd, T>(arr: &mut [O], co_arr: &mut [T]) {
    for i in 1..arr.len() {
        let mut j = i;
        while j > 0 && arr[j] < arr[j - 1] {
            arr.swap(j, j - 1);
            co_arr.swap(j, j - 1);
            j -= 1;
        }
    }
}

/// Find a local minima of *f* using the Nelder-Mead algorithm starting at *p0*
///
/// # Example
/// ```
/// use num::Zero;
///
/// use scialg::vector::Vector;
/// use scialg::optimize::nelder_mead;
///
/// let func = |v: Vector<2>| v.coeff[0] * v.coeff[0] + v.coeff[1] * v.coeff[1];
/// let minima = nelder_mead(func, &Vector::new(&[3.0, 5.0]), 100);
///
/// assert!((minima - Vector::zero()).length() < 1e-5);
/// ```
///
/// # References
///  - [Wikipedia: Nelder-Mead method](https://en.wikipedia.org/wiki/Nelder%E2%80%93Mead_method)
///  - [Wikipedia: Centroid](https://en.wikipedia.org/wiki/Centroid)
pub fn nelder_mead<const N: usize>(
    f: fn(Vector<N>) -> f64,
    p0: &Vector<N>,
    max_iter: usize,
) -> Vector<N> {
    let alpha = 1.0;
    let gamma = 2.0;
    let rho = 0.5;
    let sigma = 0.5;

    // starting simplex
    let mut xs: Vec<Vector<N>> = vec![*p0];
    for i in 0..N {
        let mut coeffs = p0.coeff;
        coeffs[i] += 1.0;
        xs.push(Vector { coeff: coeffs });
    }

    for _ in 0..max_iter {
        let mut fs: Vec<f64> = xs.iter().map(|x| f(*x)).collect();
        co_sort(&mut fs, &mut xs);

        // compute centoid
        let mut x0 = Vector::zero();
        for i in 0..N {
            x0 = x0 + xs[i];
        }
        x0 = x0 / N as f64;

        // reflect
        let xr = x0 + (x0 - xs[N]) * alpha;
        let fr = f(xr);
        if fs[0] <= fr && fr <= fs[N - 1] {
            xs[N] = xr;
            continue;
        }

        // expand
        if fr < fs[0] {
            let xe = x0 + (xr - x0) * gamma;
            let fe = f(xe);
            if fe < fr {
                xs[N] = xe;
            } else {
                xs[N] = xr;
            }
            continue;
        }

        // contract
        if fr < fs[N] {
            // contract outside
            let xc = x0 + (xr - x0) * rho;
            let fc = f(xc);
            if fc < fr {
                xs[N] = xc;
                continue;
            }
        } else {
            // contract inside
            let xc = x0 + (xs[N] - x0) * rho;
            let fc = f(xc);
            if fc < fs[N] {
                xs[N] = xc;
                continue;
            }
        }

        // shrink
        for i in 0..xs.len() {
            xs[i] = xs[0] + (xs[i] - xs[0]) * sigma;
        }
    }

    // compute centoid
    let mut x0 = Vector::zero();
    for i in 0..N {
        x0 = x0 + xs[i];
    }
    x0 / N as f64
}