metalforge 0.1.6

//! L-SHADE — Success-History based Adaptive DE with Linear Population Size
//! Reduction (Tanabe & Fukunaga 2014, CEC).
//!
//! The reference adaptive Differential Evolution: the SHADE/L-SHADE lineage won
//! the CEC single-objective competitions repeatedly from 2014 on, and JADE /
//! SaDE / SHADE / L-SHADE are the canonical baselines a continuous optimizer is
//! expected to beat. It upgrades plain `rand/1/bin` ([`De`](super::De)) on three
//! fronts:
//!
//! - **`current-to-pbest/1`** mutation with an **external archive** of replaced
//!   parents (greedy yet diverse).
//! - **Success-history adaptation** of the scaling factor `F` and crossover rate
//!   `CR` from a memory of the values that produced improvements (weighted
//!   Lehmer means), so the user never tunes them.
//! - **Linear Population Size Reduction (LPSR)**: the population shrinks from
//!   `N_init` to 4 over the evaluation budget, shifting explore → exploit.
//!
//! Implements [`Optimizer`] and is deterministic for a given seed.

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

/// L-SHADE configuration.
#[derive(Debug, Clone, Copy)]
pub struct LShade {
    /// Initial population size `N_init`; `None` uses the paper's `18 · dim`.
    pub init_pop: Option<usize>,
    /// Success-history memory size `H` (paper default 6).
    pub memory: usize,
    /// `p` for `current-to-pbest`: the top `p·N` individuals are pbest
    /// candidates (paper default 0.11).
    pub p_best: f64,
    /// Archive size as a multiple of the current population (paper default 2.6).
    pub archive_rate: f64,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for LShade {
    fn default() -> Self {
        LShade {
            init_pop: None,
            memory: 6,
            p_best: 0.11,
            archive_rate: 2.6,
            seed: 42,
        }
    }
}

/// Minimum population size L-SHADE shrinks to (the `current-to-pbest/1` donor
/// count).
const N_MIN: usize = 4;

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

    fn optimize(&self, problem: &dyn Problem, term: &Termination) -> Report {
        let bounds = problem.bounds();
        let dim = bounds.len();
        let mut rng = Rng::new(self.seed);

        let n_init = self.init_pop.unwrap_or(18 * dim).max(N_MIN);
        let h = self.memory.max(1);
        let max_nfe = term.max_evaluations.max(1);

        let eval = |x: &[f64], problem: &dyn Problem| -> f64 {
            let v = problem.objective(x);
            if v.is_finite() {
                v
            } else {
                f64::INFINITY
            }
        };

        // Success-history memories. A NaN cell is the "terminal" CR (⊥), set
        // when the only successful trials used CR = 0 (separable landscapes).
        let mut m_cr = vec![0.5f64; h];
        let mut m_f = vec![0.5f64; h];
        let mut k_pos = 0usize;

        // Initial population.
        let mut pop: Vec<Vec<f64>> = Vec::with_capacity(n_init);
        let mut fit: Vec<f64> = Vec::with_capacity(n_init);
        let mut archive: Vec<Vec<f64>> = Vec::new();
        let mut best = Solution {
            x: vec![0.0; dim],
            value: f64::INFINITY,
        };
        let mut nfe = 0usize;
        for _ in 0..n_init {
            let x: Vec<f64> = bounds
                .iter()
                .map(|&(lo, hi)| rng.uniform_in(lo, hi))
                .collect();
            let f = eval(&x, problem);
            nfe += 1;
            if f < best.value {
                best = Solution {
                    x: x.clone(),
                    value: f,
                };
            }
            pop.push(x);
            fit.push(f);
        }
        let mut n = pop.len();

        let mut trial = vec![0.0; dim];
        'outer: while term.reason(nfe, best.value).is_none() {
            // Rank for pbest selection.
            let mut ranked: Vec<usize> = (0..n).collect();
            ranked.sort_by(|&a, &b| {
                fit[a]
                    .partial_cmp(&fit[b])
                    .unwrap_or(std::cmp::Ordering::Equal)
            });
            let p_num = ((self.p_best * n as f64).round() as usize).clamp(2, n);

            let mut succ_cr: Vec<f64> = Vec::new();
            let mut succ_f: Vec<f64> = Vec::new();
            let mut delta: Vec<f64> = Vec::new();

            for i in 0..n {
                if term.reason(nfe, best.value).is_some() {
                    break 'outer;
                }
                let r = rng.index(h);
                let cr = if m_cr[r].is_nan() {
                    0.0
                } else {
                    (m_cr[r] + 0.1 * rng.normal()).clamp(0.0, 1.0)
                };
                // F ~ Cauchy(M_F, 0.1), resampled while ≤ 0, truncated to 1.
                let scale = loop {
                    let v = m_f[r] + 0.1 * cauchy(&mut rng);
                    if v > 0.0 {
                        break v.min(1.0);
                    }
                };

                let pbest = ranked[rng.index(p_num)];
                let r1 = loop {
                    let z = rng.index(n);
                    if z != i {
                        break z;
                    }
                };
                let na = archive.len();
                // r2 from population ∪ archive, distinct from i and r1.
                let r2 = loop {
                    let z = rng.index(n + na);
                    if z >= n || (z != i && z != r1) {
                        break z;
                    }
                };
                let x_r2: &[f64] = if r2 < n { &pop[r2] } else { &archive[r2 - n] };

                // current-to-pbest/1 mutation + binomial crossover, repaired.
                let jrand = rng.index(dim);
                for j in 0..dim {
                    let (lo, hi) = bounds[j];
                    if rng.uniform() <= cr || j == jrand {
                        let mut v = pop[i][j]
                            + scale * (pop[pbest][j] - pop[i][j])
                            + scale * (pop[r1][j] - x_r2[j]);
                        // L-SHADE bound repair: midpoint toward the parent.
                        if v < lo {
                            v = (lo + pop[i][j]) / 2.0;
                        } else if v > hi {
                            v = (hi + pop[i][j]) / 2.0;
                        }
                        trial[j] = v;
                    } else {
                        trial[j] = pop[i][j];
                    }
                }

                let tf = eval(&trial, problem);
                nfe += 1;
                if tf < best.value {
                    best = Solution {
                        x: trial.clone(),
                        value: tf,
                    };
                }
                // Greedy selection; improvements feed the archive and history.
                if tf <= fit[i] {
                    if tf < fit[i] {
                        succ_cr.push(cr);
                        succ_f.push(scale);
                        delta.push(fit[i] - tf);
                        archive.push(pop[i].clone());
                    }
                    pop[i].copy_from_slice(&trial);
                    fit[i] = tf;
                }
            }

            // Trim the archive to its size cap by random removal.
            let arc_size = (self.archive_rate * n as f64).round() as usize;
            while archive.len() > arc_size {
                let idx = rng.index(archive.len());
                archive.swap_remove(idx);
            }

            // Update one memory cell from this generation's successes.
            let total: f64 = delta.iter().sum();
            if !succ_f.is_empty() && total > 0.0 {
                let w: Vec<f64> = delta.iter().map(|d| d / total).collect();
                m_f[k_pos] = lehmer(&w, &succ_f);
                let max_cr = succ_cr.iter().copied().fold(f64::NEG_INFINITY, f64::max);
                m_cr[k_pos] = if m_cr[k_pos].is_nan() || max_cr == 0.0 {
                    f64::NAN
                } else {
                    lehmer(&w, &succ_cr)
                };
                k_pos = (k_pos + 1) % h;
            }

            // Linear Population Size Reduction.
            let target =
                ((N_MIN as f64 - n_init as f64) / max_nfe as f64) * nfe as f64 + n_init as f64;
            let n_new = (target.round() as usize).clamp(N_MIN, n);
            if n_new < n {
                let mut ranked2: Vec<usize> = (0..n).collect();
                ranked2.sort_by(|&a, &b| {
                    fit[a]
                        .partial_cmp(&fit[b])
                        .unwrap_or(std::cmp::Ordering::Equal)
                });
                ranked2.truncate(n_new);
                let new_pop: Vec<Vec<f64>> = ranked2.iter().map(|&i| pop[i].clone()).collect();
                let new_fit: Vec<f64> = ranked2.iter().map(|&i| fit[i]).collect();
                pop = new_pop;
                fit = new_fit;
                n = n_new;
                let arc2 = (self.archive_rate * n as f64).round() as usize;
                while archive.len() > arc2 {
                    let idx = rng.index(archive.len());
                    archive.swap_remove(idx);
                }
            }
        }

        let stop = term
            .reason(nfe, best.value)
            .unwrap_or(StopReason::BudgetExhausted);
        Report {
            solution: best,
            stop,
            evaluations: nfe,
        }
    }
}

/// Weighted Lehmer mean `Σ wₖ sₖ² / Σ wₖ sₖ` (the SHADE memory update).
fn lehmer(w: &[f64], s: &[f64]) -> f64 {
    let num: f64 = w.iter().zip(s).map(|(wk, sk)| wk * sk * sk).sum();
    let den: f64 = w.iter().zip(s).map(|(wk, sk)| wk * sk).sum();
    if den != 0.0 {
        num / den
    } else {
        0.5
    }
}

/// A standard Cauchy(0, 1) sample.
fn cauchy(rng: &mut Rng) -> f64 {
    (std::f64::consts::PI * (rng.uniform() - 0.5)).tan()
}