Skip to main content

forge_core/algo/
de.rs

1//! Differential Evolution — DE/rand/1/bin (Storn & Price 1997, J. Global
2//! Optimization 11, 341–359).
3//!
4//! A population-based optimizer: each target vector is challenged by a trial
5//! built from a scaled difference of two random members added to a third
6//! (`rand/1`), recombined with the target by binomial crossover (`bin`), and
7//! kept only if it improves. Simple, robust on continuous landscapes, and a
8//! natural complement to the geoscientific DDS/SCE-UA pair.
9
10use super::{clamp, sample, Evaluator, Optimizer};
11use crate::problem::Problem;
12use crate::rng::Rng;
13use crate::solution::{Report, Solution};
14use crate::termination::Termination;
15
16/// Differential Evolution configuration.
17#[derive(Debug, Clone, Copy)]
18pub struct De {
19    /// Population size (≥ 4; `rand/1` needs three distinct donors per target).
20    pub pop_size: usize,
21    /// Differential weight `F`, typically in `[0.4, 1.0]`.
22    pub f: f64,
23    /// Crossover probability `CR` in `[0, 1]`.
24    pub cr: f64,
25    /// RNG seed; same seed + same problem + same budget ⇒ same result.
26    pub seed: u64,
27}
28
29impl Default for De {
30    fn default() -> Self {
31        De {
32            pop_size: 30,
33            f: 0.8,
34            cr: 0.9,
35            seed: 42,
36        }
37    }
38}
39
40impl Optimizer for De {
41    fn with_seed(&self, seed: u64) -> Self {
42        De { seed, ..*self }
43    }
44
45    /// Minimizes `problem` within its bounds using DE/rand/1/bin.
46    fn optimize(&self, problem: &dyn Problem, term: &Termination) -> Report {
47        crate::problem::validate(problem).unwrap_or_else(|e| panic!("De: invalid problem: {e}"));
48        let bounds = problem.bounds();
49        let dim = bounds.len();
50        let np = self.pop_size.max(4);
51        let mut rng = Rng::new(self.seed);
52
53        // Initialize and evaluate the population, seeding the evaluator.
54        let first_x = sample(bounds, &mut rng);
55        let first_v = finite_or_worst(problem.objective(&first_x));
56        let mut ev = Evaluator::new(
57            problem,
58            term,
59            Solution {
60                x: first_x.clone(),
61                value: first_v,
62            },
63        );
64
65        let mut pop: Vec<Vec<f64>> = Vec::with_capacity(np);
66        let mut fit: Vec<f64> = Vec::with_capacity(np);
67        pop.push(first_x);
68        fit.push(first_v);
69        for _ in 1..np {
70            if ev.done() {
71                break;
72            }
73            let x = sample(bounds, &mut rng);
74            let v = finite_or_worst(ev.eval(&x));
75            pop.push(x);
76            fit.push(v);
77        }
78        let np = pop.len(); // may be smaller if the budget was tiny
79
80        let mut trial = vec![0.0; dim];
81        while !ev.done() {
82            for i in 0..np {
83                if ev.done() {
84                    break;
85                }
86                // Three distinct donors, all different from the target i.
87                let (a, b, c) = three_distinct(np, i, &mut rng);
88                // Guarantee at least one mutant component (jrand).
89                let jrand = rng.index(dim);
90                for j in 0..dim {
91                    let (lo, hi) = bounds[j];
92                    if rng.uniform() < self.cr || j == jrand {
93                        let mutant = pop[a][j] + self.f * (pop[b][j] - pop[c][j]);
94                        trial[j] = reflect(mutant, lo, hi);
95                    } else {
96                        trial[j] = pop[i][j];
97                    }
98                }
99
100                let tv = finite_or_worst(ev.eval(&trial));
101                // Greedy selection: replace the target if the trial is no worse.
102                if tv <= fit[i] {
103                    pop[i].copy_from_slice(&trial);
104                    fit[i] = tv;
105                }
106            }
107        }
108
109        ev.finish()
110    }
111}
112
113#[inline]
114fn finite_or_worst(v: f64) -> f64 {
115    if v.is_finite() {
116        v
117    } else {
118        f64::INFINITY
119    }
120}
121
122/// Reflects a value back into `[lo, hi]` once at the violated bound, then clamps
123/// — keeps mutants in-domain without bunching them on the boundary.
124fn reflect(x: f64, lo: f64, hi: f64) -> f64 {
125    let mut v = x;
126    if v < lo {
127        v = lo + (lo - v);
128    } else if v > hi {
129        v = hi - (v - hi);
130    }
131    clamp(v, lo, hi)
132}
133
134/// Three indices in `0..np`, mutually distinct and all different from `target`.
135/// Requires `np >= 4`.
136fn three_distinct(np: usize, target: usize, rng: &mut Rng) -> (usize, usize, usize) {
137    let pick = |rng: &mut Rng, exclude: &[usize]| loop {
138        let r = rng.index(np);
139        if r != target && !exclude.contains(&r) {
140            return r;
141        }
142    };
143    let a = pick(rng, &[]);
144    let b = pick(rng, &[a]);
145    let c = pick(rng, &[a, b]);
146    (a, b, c)
147}