metalforge 0.1.1

//! Differential Evolution — DE/rand/1/bin (Storn & Price 1997, J. Global
//! Optimization 11, 341–359).
//!
//! A population-based optimizer: each target vector is challenged by a trial
//! built from a scaled difference of two random members added to a third
//! (`rand/1`), recombined with the target by binomial crossover (`bin`), and
//! kept only if it improves. Simple, robust on continuous landscapes, and a
//! natural complement to the geoscientific DDS/SCE-UA pair.

use super::{clamp, sample, Evaluator, Optimizer};
use crate::problem::Problem;
use crate::rng::Rng;
use crate::solution::{Report, Solution};
use crate::termination::Termination;

/// Differential Evolution configuration.
#[derive(Debug, Clone, Copy)]
pub struct De {
    /// Population size (≥ 4; `rand/1` needs three distinct donors per target).
    pub pop_size: usize,
    /// Differential weight `F`, typically in `[0.4, 1.0]`.
    pub f: f64,
    /// Crossover probability `CR` in `[0, 1]`.
    pub cr: f64,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for De {
    fn default() -> Self {
        De {
            pop_size: 30,
            f: 0.8,
            cr: 0.9,
            seed: 42,
        }
    }
}

impl Optimizer for De {
    fn with_seed(&self, seed: u64) -> Self {
        De { seed, ..*self }
    }

    /// Minimizes `problem` within its bounds using DE/rand/1/bin.
    fn optimize(&self, problem: &dyn Problem, term: &Termination) -> Report {
        let bounds = problem.bounds();
        let dim = bounds.len();
        let np = self.pop_size.max(4);
        let mut rng = Rng::new(self.seed);

        // Initialize and evaluate the population, seeding the evaluator.
        let first_x = sample(bounds, &mut rng);
        let first_v = finite_or_worst(problem.objective(&first_x));
        let mut ev = Evaluator::new(
            problem,
            term,
            Solution {
                x: first_x.clone(),
                value: first_v,
            },
        );

        let mut pop: Vec<Vec<f64>> = Vec::with_capacity(np);
        let mut fit: Vec<f64> = Vec::with_capacity(np);
        pop.push(first_x);
        fit.push(first_v);
        for _ in 1..np {
            if ev.done() {
                break;
            }
            let x = sample(bounds, &mut rng);
            let v = finite_or_worst(ev.eval(&x));
            pop.push(x);
            fit.push(v);
        }
        let np = pop.len(); // may be smaller if the budget was tiny

        let mut trial = vec![0.0; dim];
        while !ev.done() {
            for i in 0..np {
                if ev.done() {
                    break;
                }
                // Three distinct donors, all different from the target i.
                let (a, b, c) = three_distinct(np, i, &mut rng);
                // Guarantee at least one mutant component (jrand).
                let jrand = rng.index(dim);
                for j in 0..dim {
                    let (lo, hi) = bounds[j];
                    if rng.uniform() < self.cr || j == jrand {
                        let mutant = pop[a][j] + self.f * (pop[b][j] - pop[c][j]);
                        trial[j] = reflect(mutant, lo, hi);
                    } else {
                        trial[j] = pop[i][j];
                    }
                }

                let tv = finite_or_worst(ev.eval(&trial));
                // Greedy selection: replace the target if the trial is no worse.
                if tv <= fit[i] {
                    pop[i].copy_from_slice(&trial);
                    fit[i] = tv;
                }
            }
        }

        ev.finish()
    }
}

#[inline]
fn finite_or_worst(v: f64) -> f64 {
    if v.is_finite() {
        v
    } else {
        f64::INFINITY
    }
}

/// Reflects a value back into `[lo, hi]` once at the violated bound, then clamps
/// — keeps mutants in-domain without bunching them on the boundary.
fn reflect(x: f64, lo: f64, hi: f64) -> f64 {
    let mut v = x;
    if v < lo {
        v = lo + (lo - v);
    } else if v > hi {
        v = hi - (v - hi);
    }
    clamp(v, lo, hi)
}

/// Three indices in `0..np`, mutually distinct and all different from `target`.
/// Requires `np >= 4`.
fn three_distinct(np: usize, target: usize, rng: &mut Rng) -> (usize, usize, usize) {
    let pick = |rng: &mut Rng, exclude: &[usize]| loop {
        let r = rng.index(np);
        if r != target && !exclude.contains(&r) {
            return r;
        }
    };
    let a = pick(rng, &[]);
    let b = pick(rng, &[a]);
    let c = pick(rng, &[a, b]);
    (a, b, c)
}