metalforge 0.3.0

//! 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: both control knobs
//! are driven by the fraction of trials that improved in the previous
//! generation (`SR`). The published rules are
//!
//! ```text
//! mF = 0.4 + 0.25 · tanh(5 · SR)          F  ~ Normal(mF, 0.02), resampled to [0, 1]
//! pb = max(2, ⌊0.7 · exp(−7 · SR) · N⌋)   (elite pool for the p-best donor)
//! ```
//!
//! so a high success rate raises `F` (bolder steps) while *shrinking* the
//! elite pool (more selective pressure); a low success rate lowers `F` and
//! widens the pool. Structurally L-SRTDE inherits **L-NTADE**'s two-population
//! scheme (Stanovov & Semenkin 2022) rather than L-SHADE's single population:
//!
//! - a **front** population that trials compete against (a random front member
//!   is the target each step) and into which successful trials are inserted
//!   circularly, and
//! - a **newest** population that grows with successful trials during the
//!   generation and is truncated back to the best `N` at the end.
//!
//! The mutant is `front[k] + F·(newest[pbest] − front[k]) + F·(front[r1] −
//! newest[r2])` with `r1` drawn from the front by a rank-based exponential
//! distribution (`weight ∝ exp(−3·rank/N)`) and `r2` uniform from the newest
//! population. There is **no external archive**. `CR` uses a small
//! success-history memory (size 5, init 1.0) sampled `Normal(m_CR, 0.05)`,
//! stores the *repaired* CR (realized crossover fraction), and updates by a
//! smoothed fitness-delta-weighted Lehmer mean. Out-of-bounds components are
//! re-drawn uniformly in the box. Population size follows LPSR down to 4.
//!
//! Fidelity note: constants and structure follow the author's public
//! reference implementation (github.com/VladimirStanovov/L-SRTDE_CEC-2024,
//! reimplemented from the algorithm description — no code copied). Known
//! deviations: the budget is enforced per evaluation (the crate contract)
//! instead of per generation, and the circular front-insertion index is
//! wrapped after LPSR shrinks the front. 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 the paper's `20 · dim`.
    pub init_pop: Option<usize>,
    /// Success-history memory size `H` for the crossover rate.
    pub memory: usize,
    /// 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,
            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 {
        crate::problem::validate(problem)
            .unwrap_or_else(|e| panic!("LSrtde: invalid problem: {e}"));
        let bounds = problem.bounds();
        let dim = bounds.len();
        let mut rng = Rng::new(self.seed);

        let n_init = self.init_pop.unwrap_or(20 * 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
            }
        };

        // CR success-history memory, init 1.0 (F is success-rate driven).
        let mut m_cr = vec![1.0f64; h];
        let mut k_pos = 0usize;

        // Newest population (grows with successes, truncated per generation).
        let mut pop: Vec<Vec<f64>> = Vec::with_capacity(2 * n_init);
        let mut fit: Vec<f64> = Vec::with_capacity(2 * n_init);
        let mut best = Solution {
            x: vec![0.0; dim],
            value: f64::INFINITY,
        };
        let mut nfe = 0usize;
        for _ in 0..n_init {
            if term.reason(nfe, best.value).is_some() {
                break;
            }
            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);
        }

        // Front population: the initial individuals sorted best-first.
        let mut order: Vec<usize> = (0..pop.len()).collect();
        order.sort_by(|&a, &b| {
            fit[a]
                .partial_cmp(&fit[b])
                .unwrap_or(std::cmp::Ordering::Equal)
        });
        let mut front: Vec<Vec<f64>> = order.iter().map(|&i| pop[i].clone()).collect();
        let mut front_fit: Vec<f64> = order.iter().map(|&i| fit[i]).collect();
        // Keep the newest population sorted too (matches the reference's
        // post-truncation state at the top of each generation).
        pop = front.clone();
        fit = front_fit.clone();

        let mut n = front.len(); // current front size (LPSR shrinks it)
        let n_init_actual = n.max(1);
        let mut pf_index = 0usize; // circular front-insertion cursor
        let mut success_rate: f64 = 0.5;
        let mut trial = vec![0.0; dim];

        'outer: while n >= N_MIN && term.reason(nfe, best.value).is_none() {
            let m_f = 0.4 + 0.25 * (5.0 * success_rate).tanh();
            // Elite-pool size shrinks as the success rate rises.
            let p_num = ((0.7 * (-7.0 * success_rate).exp() * n as f64) as usize).clamp(2, n);

            // Rank orders: `ranked` over the newest population, `ranked_f`
            // over the front (best first).
            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 mut ranked_f: Vec<usize> = (0..n).collect();
            ranked_f.sort_by(|&a, &b| {
                front_fit[a]
                    .partial_cmp(&front_fit[b])
                    .unwrap_or(std::cmp::Ordering::Equal)
            });
            // Rank-based exponential weights for the front donor: cumulative
            // distribution over exp(-3·rank/n).
            let cum: Vec<f64> = {
                let mut acc = 0.0;
                let w: Vec<f64> = (0..n).map(|i| (-3.0 * i as f64 / n as f64).exp()).collect();
                let total: f64 = w.iter().sum();
                w.iter()
                    .map(|v| {
                        acc += v / total;
                        acc
                    })
                    .collect()
            };
            let rank_sample = |rng: &mut Rng| -> usize {
                let u = rng.uniform();
                match cum.iter().position(|&c| u < c) {
                    Some(r) => r,
                    None => n - 1,
                }
            };

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

            for _ in 0..n {
                if term.reason(nfe, best.value).is_some() {
                    break 'outer;
                }
                // Random front member is the target this step.
                let chosen = rng.index(n);
                let cell = rng.index(h);
                let prand = loop {
                    let z = ranked[rng.index(p_num)];
                    if z != chosen {
                        break z;
                    }
                };
                let r1 = loop {
                    let z = ranked_f[rank_sample(&mut rng)];
                    if z != prand {
                        break z;
                    }
                };
                let r2 = loop {
                    let z = ranked[rng.index(n)];
                    if z != prand && z != r1 {
                        break z;
                    }
                };

                // F ~ Normal(mF, 0.02) resampled into [0, 1].
                let scale = loop {
                    let v = m_f + 0.02 * rng.normal();
                    if (0.0..=1.0).contains(&v) {
                        break v;
                    }
                };
                let cr = (m_cr[cell] + 0.05 * rng.normal()).clamp(0.0, 1.0);

                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 = front[chosen][j]
                            + scale * (pop[prand][j] - front[chosen][j])
                            + scale * (front[r1][j] - pop[r2][j]);
                        if v < lo || v > hi {
                            // Reference rule: re-draw uniformly in the box.
                            v = rng.uniform_in(lo, hi);
                        }
                        trial[j] = v;
                        from_mutant += 1;
                    } else {
                        trial[j] = front[chosen][j];
                    }
                }

                let tf = eval(&trial, problem);
                nfe += 1;
                if tf < best.value {
                    best = Solution {
                        x: trial.clone(),
                        value: tf,
                    };
                }
                // Ties count as successes (reference: `<=`).
                if tf <= front_fit[chosen] {
                    // Repaired CR: the realized crossover fraction.
                    succ_cr.push(from_mutant as f64 / dim as f64);
                    let d = (front_fit[chosen] - tf).abs();
                    delta.push(if d.is_finite() { d } else { f64::MAX });
                    successes += 1;
                    // Grow the newest population; insert circularly into the front.
                    pop.push(trial.clone());
                    fit.push(tf);
                    front[pf_index].copy_from_slice(&trial);
                    front_fit[pf_index] = tf;
                    pf_index = (pf_index + 1) % n;
                }
            }

            success_rate = successes as f64 / n as f64;

            // CR memory: smoothed weighted Lehmer mean (fallback 1.0), only
            // advancing the cursor when the generation produced successes.
            if !succ_cr.is_empty() {
                let total: f64 = delta.iter().sum();
                let w: Vec<f64> = if total.is_finite() && total > 0.0 {
                    delta.iter().map(|d| d / total).collect()
                } else {
                    vec![1.0 / succ_cr.len() as f64; succ_cr.len()]
                };
                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();
                let lehmer = if den.abs() > 1e-8 { num / den } else { 1.0 };
                m_cr[k_pos] = 0.5 * (lehmer + m_cr[k_pos]);
                k_pos = (k_pos + 1) % h;
            }

            // LPSR on the front (truncation toward N_MIN = 4), dropping the
            // worst front members; the newest population is truncated to the
            // best `n` of (previous newest + this generation's successes).
            let target = ((N_MIN as f64 - n_init_actual as f64) / max_nfe as f64) * nfe as f64
                + n_init_actual as f64;
            let n_new = (target as usize).clamp(N_MIN, n);
            if n_new < n {
                // Remove the worst front members, preserving insertion order.
                let mut keep: Vec<usize> = (0..n).collect();
                keep.sort_by(|&a, &b| {
                    front_fit[a]
                        .partial_cmp(&front_fit[b])
                        .unwrap_or(std::cmp::Ordering::Equal)
                });
                keep.truncate(n_new);
                keep.sort_unstable();
                front = keep
                    .iter()
                    .map(|&i| std::mem::take(&mut front[i]))
                    .collect();
                front_fit = keep.iter().map(|&i| front_fit[i]).collect();
                n = n_new;
                if pf_index >= n {
                    pf_index = 0;
                }
            }
            // Truncate the newest population to the best `n`.
            if pop.len() > n {
                let mut order: Vec<usize> = (0..pop.len()).collect();
                order.sort_by(|&a, &b| {
                    fit[a]
                        .partial_cmp(&fit[b])
                        .unwrap_or(std::cmp::Ordering::Equal)
                });
                order.truncate(n);
                pop = order.iter().map(|&i| std::mem::take(&mut pop[i])).collect();
                fit = order.iter().map(|&i| fit[i]).collect();
            }
        }

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