metalforge 0.1.6

//! L-SRTDE — Success Rate-based adaptive Differential Evolution
//! (Stanovov & Semenkin 2024, CEC 2024 bound-constrained competition **winner**).
//!
//! L-SRTDE marks the shift in adaptive DE from *fitness-magnitude* / success-
//! history adaptation toward **success-rate** adaptation: the scaling factor's
//! mean is driven by the fraction of trials that improved in the previous
//! generation (`SR`), not by a per-parameter memory. The published, defining
//! rule is
//!
//! ```text
//! mF = 0.4 + 0.25 · tanh(5 · SR)
//! ```
//!
//! so a high success rate raises `F` (bolder, exploratory steps) and a low one
//! lowers it (finer, exploitative steps). `F` is then drawn `Cauchy(mF, 0.1)`.
//! The pbest fraction is likewise coupled to `SR` (selective pressure rises as
//! `SR` drops). The rest of the scaffold is shared with L-SHADE: `current-to-
//! pbest/1` with an external archive, **success-history** adaptation of `CR`
//! (sampled with the tighter dispersion `0.05` and stored *repaired*), and
//! Linear Population Size Reduction.
//!
//! Fidelity note: the `mF(SR)` rule is the paper's; auxiliary constants (memory
//! size, archive rate) and the exact `pb(SR)` schedule — not all public in the
//! paywalled paper — follow the SHADE family and the paper's described
//! behavior. Implements [`Optimizer`]; 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-SRTDE configuration.
#[derive(Debug, Clone, Copy)]
pub struct LSrtde {
    /// Initial population size `N_init`; `None` uses `18 · dim`.
    pub init_pop: Option<usize>,
    /// Success-history memory size `H` for the crossover rate.
    pub memory: usize,
    /// Archive size as a multiple of the current population.
    pub archive_rate: f64,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for LSrtde {
    fn default() -> Self {
        LSrtde {
            init_pop: None,
            memory: 5,
            archive_rate: 2.1,
            seed: 42,
        }
    }
}

const N_MIN: usize = 4;

impl Optimizer for LSrtde {
    fn with_seed(&self, seed: u64) -> Self {
        LSrtde { 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 memory for CR only (F is success-rate driven). A NaN
        // cell is the terminal CR (⊥), as in L-SHADE.
        let mut m_cr = vec![0.5f64; h];
        let mut k_pos = 0usize;

        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();

        // Success rate of the previous generation (seed at 0.5).
        let mut success_rate: f64 = 0.5;
        let mut trial = vec![0.0; dim];

        'outer: while term.reason(nfe, best.value).is_none() {
            // Success-rate-driven control (Stanovov & Semenkin 2024).
            let m_f = 0.4 + 0.25 * (5.0 * success_rate).tanh();
            // pbest fraction rises with SR (pressure increases as SR drops).
            let pb = 0.1 + 0.2 * success_rate;

            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 = ((pb * n as f64).round() as usize).clamp(2, n);

            let mut succ_cr: Vec<f64> = Vec::new();
            let mut delta: Vec<f64> = Vec::new();
            let mut successes = 0usize;

            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.05 * rng.normal()).clamp(0.0, 1.0)
                };
                let scale = loop {
                    let v = m_f + 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();
                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] };

                let jrand = rng.index(dim);
                let mut from_mutant = 0usize;
                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]);
                        if v < lo {
                            v = (lo + pop[i][j]) / 2.0;
                        } else if v > hi {
                            v = (hi + pop[i][j]) / 2.0;
                        }
                        trial[j] = v;
                        from_mutant += 1;
                    } else {
                        trial[j] = pop[i][j];
                    }
                }

                let tf = eval(&trial, problem);
                nfe += 1;
                if tf < best.value {
                    best = Solution {
                        x: trial.clone(),
                        value: tf,
                    };
                }
                if tf <= fit[i] {
                    if tf < fit[i] {
                        // Repaired CR: the actual fraction taken from the mutant
                        // (Gong et al. 2014), stored instead of the sampled CR.
                        succ_cr.push(from_mutant as f64 / dim as f64);
                        delta.push(fit[i] - tf);
                        archive.push(pop[i].clone());
                        successes += 1;
                    }
                    pop[i].copy_from_slice(&trial);
                    fit[i] = tf;
                }
            }

            // Success rate for the next generation drives mF and pb.
            success_rate = successes as f64 / n as f64;

            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);
            }

            // CR memory update (weighted Lehmer mean of repaired CRs).
            let total: f64 = delta.iter().sum();
            if !succ_cr.is_empty() && total > 0.0 {
                let w: Vec<f64> = delta.iter().map(|d| d / total).collect();
                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 {
                    let num: f64 = w.iter().zip(&succ_cr).map(|(wk, c)| wk * c * c).sum();
                    let den: f64 = w.iter().zip(&succ_cr).map(|(wk, c)| wk * c).sum();
                    if den != 0.0 {
                        num / den
                    } else {
                        0.5
                    }
                };
                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);
                pop = ranked2.iter().map(|&i| pop[i].clone()).collect();
                fit = ranked2.iter().map(|&i| fit[i]).collect();
                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,
        }
    }
}

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