rustyqlib/core/optimization/
differential_evolution.rs

1//! Differential evolution (DE/rand/1/bin): seeded, bounded global
2//! search. The tool for multimodal calibration landscapes — find the
3//! basin globally, then polish with BFGS or Levenberg-Marquardt.
4
5use super::{OptimConfig, OptimResult};
6
7const DIFFERENTIAL_WEIGHT: f64 = 0.8; // F
8const CROSSOVER_RATE: f64 = 0.9; // CR
9
10/// Minimize `f` inside the per-parameter `bounds` box with
11/// DE/rand/1/bin. Deterministic for a given `seed`; the population size
12/// is `max(15, 10 * dim)`. Converged when the population's value spread
13/// falls below `tol`; `iterations` counts generations.
14pub fn differential_evolution(
15    cfg: &OptimConfig,
16    f: &dyn Fn(&[f64]) -> f64,
17    bounds: &[(f64, f64)],
18    seed: u64,
19) -> OptimResult {
20    let dim = bounds.len();
21    assert!(dim > 0, "bounds must give at least one parameter");
22    for &(lo, hi) in bounds {
23        assert!(lo < hi, "each bound needs lo < hi");
24    }
25    let np = (10 * dim).max(15);
26    let mut rng = Xorshift64Star::new(seed);
27
28    // random initial population in the box
29    let mut pop: Vec<Vec<f64>> = (0..np)
30        .map(|_| bounds.iter().map(|&(lo, hi)| lo + (hi - lo) * rng.uniform()).collect())
31        .collect();
32    let mut values: Vec<f64> = pop.iter().map(|x| f(x)).collect();
33
34    for gen in 1..=cfg.max_iter {
35        for i in 0..np {
36            // three distinct partners, none equal to i
37            let mut pick = || loop {
38                let j = (rng.next() % np as u64) as usize;
39                if j != i {
40                    return j;
41                }
42            };
43            let (a, b, c) = {
44                let a = pick();
45                let b = loop {
46                    let b = pick();
47                    if b != a {
48                        break b;
49                    }
50                };
51                let c = loop {
52                    let c = pick();
53                    if c != a && c != b {
54                        break c;
55                    }
56                };
57                (a, b, c)
58            };
59            // mutate + binomial crossover (j_rand guarantees one gene moves)
60            let j_rand = (rng.next() % dim as u64) as usize;
61            let mut trial = pop[i].clone();
62            for j in 0..dim {
63                if j == j_rand || rng.uniform() < CROSSOVER_RATE {
64                    let v = pop[a][j] + DIFFERENTIAL_WEIGHT * (pop[b][j] - pop[c][j]);
65                    trial[j] = v.clamp(bounds[j].0, bounds[j].1);
66                }
67            }
68            let f_trial = f(&trial);
69            if f_trial <= values[i] {
70                pop[i] = trial;
71                values[i] = f_trial;
72            }
73        }
74        let best = values.iter().cloned().fold(f64::INFINITY, f64::min);
75        let worst = values.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
76        if worst - best <= cfg.tol * (1.0 + best.abs()) {
77            let (i, &value) =
78                values.iter().enumerate().min_by(|a, b| a.1.total_cmp(b.1)).expect("non-empty");
79            return OptimResult { x: pop[i].clone(), value, iterations: gen, converged: true };
80        }
81    }
82    let (i, &value) =
83        values.iter().enumerate().min_by(|a, b| a.1.total_cmp(b.1)).expect("non-empty");
84    OptimResult { x: pop[i].clone(), value, iterations: cfg.max_iter, converged: false }
85}
86
87/// Small deterministic RNG (xorshift64*), so runs are reproducible for a
88/// given seed without pulling in an external crate.
89struct Xorshift64Star {
90    state: u64,
91}
92
93impl Xorshift64Star {
94    fn new(seed: u64) -> Self {
95        Self { state: seed.max(1) }
96    }
97    fn next(&mut self) -> u64 {
98        let mut x = self.state;
99        x ^= x >> 12;
100        x ^= x << 25;
101        x ^= x >> 27;
102        self.state = x;
103        x.wrapping_mul(0x2545F4914F6CDD1D)
104    }
105    fn uniform(&mut self) -> f64 {
106        (self.next() >> 11) as f64 / (1u64 << 53) as f64
107    }
108}
109
110#[cfg(test)]
111mod tests {
112    use super::*;
113    use std::f64::consts::PI;
114
115    #[test]
116    fn finds_the_global_minimum_of_rastrigin() {
117        // many local minima; the global one is 0 at the origin
118        let f = |x: &[f64]| {
119            20.0 + x.iter().map(|xi| xi * xi - 10.0 * (2.0 * PI * xi).cos()).sum::<f64>()
120        };
121        let bounds = [(-5.12, 5.12), (-5.12, 5.12)];
122        let r = differential_evolution(&OptimConfig::new(1e-10, 600), &f, &bounds, 7);
123        assert!(r.value < 1e-6, "stuck at a local minimum: {r:?}");
124        assert!(r.x.iter().all(|xi| xi.abs() < 1e-3), "{:?}", r.x);
125    }
126
127    #[test]
128    fn is_deterministic_for_a_seed_and_respects_bounds() {
129        let f = |x: &[f64]| (x[0] - 0.5).powi(2) + (x[1] - 0.25).powi(2);
130        let bounds = [(0.0, 1.0), (0.0, 1.0)];
131        let cfg = OptimConfig::new(1e-12, 300);
132        let a = differential_evolution(&cfg, &f, &bounds, 123);
133        let b = differential_evolution(&cfg, &f, &bounds, 123);
134        assert_eq!(a.x, b.x, "same seed must reproduce the same run");
135        assert!(a.x.iter().all(|&v| (0.0..=1.0).contains(&v)));
136        assert!((a.x[0] - 0.5).abs() < 1e-5 && (a.x[1] - 0.25).abs() < 1e-5, "{a:?}");
137        // a different seed still finds the optimum
138        let c = differential_evolution(&cfg, &f, &bounds, 999);
139        assert!((c.x[0] - 0.5).abs() < 1e-5, "{c:?}");
140    }
141}
rustyqlib/core/optimization/differential_evolution.rs

rustyqlib/core/optimization/
differential_evolution.rs
