slain-num 0.4.0

SLAIN's own numeric crate: dense matrices, least squares, polynomials, FFT, convolution, statistics and 1-D/N-D minimisation, with AVX2 kernels.
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

pub fn golden(mut f: impl FnMut(f64) -> f64, mut a: f64, mut b: f64, tol: f64) -> (f64, f64) {
    let g = (5f64.sqrt() - 1.0) / 2.0;
    let mut c = b - g * (b - a);
    let mut d = a + g * (b - a);
    let (mut fc, mut fd) = (f(c), f(d));
    while (b - a).abs() > tol {
        if fc < fd {
            b = d;
            (d, fd) = (c, fc);
            c = b - g * (b - a);
            fc = f(c);
        } else {
            a = c;
            (c, fc) = (d, fd);
            d = a + g * (b - a);
            fd = f(d);
        }
    }
    let x = (a + b) / 2.0;
    (x, f(x))
}

pub fn nelder_mead(mut f: impl FnMut(&[f64]) -> f64, start: &[f64], step: f64, iters: usize) -> (Vec<f64>, f64) {
    let n = start.len();
    let mut pts: Vec<Vec<f64>> = vec![start.to_vec()];
    for i in 0..n {
        let mut p = start.to_vec();
        p[i] += step;
        pts.push(p);
    }
    let mut vals: Vec<f64> = pts.iter().map(|p| f(p)).collect();
    let lerp = |a: &[f64], b: &[f64], t: f64| -> Vec<f64> { a.iter().zip(b).map(|(x, y)| x + t * (y - x)).collect() };
    for _ in 0..iters {
        let mut order: Vec<usize> = (0..=n).collect();
        order.sort_by(|&i, &j| vals[i].total_cmp(&vals[j]));
        pts = order.iter().map(|&i| pts[i].clone()).collect();
        vals = order.iter().map(|&i| vals[i]).collect();
        let centroid: Vec<f64> = (0..n).map(|k| pts[..n].iter().map(|p| p[k]).sum::<f64>() / n as f64).collect();
        let worst = pts[n].clone();
        let refl = lerp(&worst, &centroid, 2.0);
        let fr = f(&refl);
        if fr < vals[0] {
            let exp = lerp(&worst, &centroid, 3.0);
            let fe = f(&exp);
            (pts[n], vals[n]) = if fe < fr { (exp, fe) } else { (refl, fr) };
        } else if fr < vals[n - 1] {
            (pts[n], vals[n]) = (refl, fr);
        } else {
            let con = lerp(&worst, &centroid, 0.5);
            let fcn = f(&con);
            if fcn < vals[n] {
                (pts[n], vals[n]) = (con, fcn);
            } else {
                for i in 1..=n {
                    pts[i] = lerp(&pts[0], &pts[i], 0.5);
                    vals[i] = f(&pts[i]);
                }
            }
        }
    }
    let best = (0..=n).min_by(|&i, &j| vals[i].total_cmp(&vals[j])).unwrap_or(0);
    (pts[best].clone(), vals[best])
}

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

    #[test]
    fn finds_minima() {
        let (x, _) = golden(|x| (x - 1.7).powi(2) + 3.0, -5.0, 5.0, 1e-8);
        assert!((x - 1.7).abs() < 1e-6);
        let (p, v) = nelder_mead(|p| (p[0] - 1.0).powi(2) + 10.0 * (p[1] + 2.0).powi(2), &[0.0, 0.0], 0.5, 400);
        assert!((p[0] - 1.0).abs() < 1e-4 && (p[1] + 2.0).abs() < 1e-4 && v < 1e-7);
    }
}