metalforge 0.3.0

//! ε-constrained L-SHADE — L-SHADE with the ε-level comparison for nonlinear
//! constraints (Takahama & Sakai's ε-constrained method, the backbone of the
//! CEC constrained-competition lineage: LSHADE44, εMAg-ES, IUDE).
//!
//! The DE machinery is exactly [`LShade`](super::LShade) (current-to-pbest/1
//! with archive, success-history CR/F, LPSR). What changes is *comparison*:
//! each candidate carries `(objective, total_violation)`, and selection and
//! pbest ranking use the **ε-comparison** with the current tolerance:
//!
//! - if both violations are ≤ ε (or exactly equal): compare by objective;
//! - otherwise: compare by violation.
//!
//! The tolerance starts at the violation of the θ-quantile of the initial
//! population (`θ = 0.05`, self-calibrating to the problem's constraint
//! scale) and shrinks as `ε(t) = ε(0)·(1 − t/T_c)^cp` with `cp = 5`,
//! `T_c = 0.8 · budget`, reaching exactly 0 afterwards — early search may
//! traverse slightly-infeasible regions, the endgame is strictly feasible.
//! Crucially the comparison is re-applied with the *current* ε every
//! generation over stored `(f, v)` pairs, which is why this lives inside the
//! algorithm instead of a problem adapter (cached transformed fitnesses
//! would go stale as ε contracts).
//!
//! The reported best is tracked with the **ε = 0** comparison throughout, so
//! the final answer is the best truly-feasible solution seen (or, if nothing
//! feasible was found, the least-violating one — check
//! [`ConstrainedProblem::total_violation`] on the result). Deterministic for
//! a given seed.
//!
//! [`ConstrainedProblem::total_violation`]: crate::constraint::ConstrainedProblem::total_violation

use crate::constraint::ConstrainedProblem;
use crate::rng::Rng;
use crate::solution::{Report, Solution, StopReason};
use crate::termination::Termination;
use std::cmp::Ordering;

/// ε-constrained L-SHADE configuration.
#[derive(Debug, Clone, Copy)]
pub struct EpsilonLShade {
    /// Initial population size `N_init`; `None` uses `18 · dim`.
    pub init_pop: Option<usize>,
    /// Success-history memory size `H`.
    pub memory: usize,
    /// `p` for `current-to-pbest` (top `p·N` are pbest candidates).
    pub p_best: f64,
    /// Archive size as a multiple of the current population.
    pub archive_rate: f64,
    /// Initial ε; `None` self-calibrates to the θ = 0.05 violation quantile of
    /// the initial population (Takahama & Sakai's rule).
    pub epsilon0: Option<f64>,
    /// Fraction of the budget after which ε is exactly 0.
    pub cutoff: f64,
    /// Shrinkage exponent `cp`.
    pub cp: f64,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for EpsilonLShade {
    fn default() -> Self {
        EpsilonLShade {
            init_pop: None,
            memory: 6,
            p_best: 0.11,
            archive_rate: 2.6,
            epsilon0: None,
            cutoff: 0.8,
            cp: 5.0,
            seed: 42,
        }
    }
}

const N_MIN: usize = 4;

/// `(objective, total_violation)` of a candidate.
type Fv = (f64, f64);

/// The ε-comparison: `a` better than `b` under tolerance `eps`.
fn eps_less(a: Fv, b: Fv, eps: f64) -> bool {
    let (fa, va) = a;
    let (fb, vb) = b;
    if (va <= eps && vb <= eps) || va == vb {
        fa < fb
    } else {
        va < vb
    }
}

/// Total order for ranking under tolerance `eps` (ties by objective).
fn eps_cmp(a: Fv, b: Fv, eps: f64) -> Ordering {
    if eps_less(a, b, eps) {
        Ordering::Less
    } else if eps_less(b, a, eps) {
        Ordering::Greater
    } else {
        Ordering::Equal
    }
}

impl EpsilonLShade {
    /// Minimizes `problem` under its nonlinear constraints within the budget.
    pub fn optimize(&self, problem: &dyn ConstrainedProblem, term: &Termination) -> Report {
        struct AsProblem<'a>(&'a dyn ConstrainedProblem);
        impl crate::problem::Problem for AsProblem<'_> {
            fn dim(&self) -> usize {
                self.0.dim()
            }
            fn bounds(&self) -> &[crate::problem::Bound] {
                self.0.bounds()
            }
            fn objective(&self, _x: &[f64]) -> f64 {
                0.0
            }
        }
        crate::problem::validate(&AsProblem(problem))
            .unwrap_or_else(|e| panic!("EpsilonLShade: 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(18 * dim).max(N_MIN);
        let h = self.memory.max(1);
        let max_nfe = term.max_evaluations.max(1);
        let t_c = ((self.cutoff.clamp(0.0, 1.0)) * max_nfe as f64) as usize;

        // Evaluate objective + violation; non-finite objectives become +inf
        // (still comparable by violation), non-finite violations become +inf.
        let eval = |x: &[f64]| -> Fv {
            let f = problem.objective(x);
            let v = problem.total_violation(x);
            (
                if f.is_finite() { f } else { f64::INFINITY },
                if v.is_finite() { v } else { f64::INFINITY },
            )
        };

        let mut m_cr = vec![0.5f64; h];
        let mut m_f = vec![0.5f64; h];
        let mut k_pos = 0usize;

        let mut pop: Vec<Vec<f64>> = Vec::with_capacity(n_init);
        let mut fit: Vec<Fv> = Vec::with_capacity(n_init);
        let mut archive: Vec<Vec<f64>> = Vec::new();
        // The reported best always uses the ε = 0 comparison.
        let mut best_x = vec![0.0; dim];
        let mut best_fv: Fv = (f64::INFINITY, f64::INFINITY);
        let mut nfe = 0usize;

        // Target checks apply only to truly feasible bests.
        let best_for_term = |b: Fv| if b.1 == 0.0 { b.0 } else { f64::INFINITY };

        for _ in 0..n_init {
            if term.reason(nfe, best_for_term(best_fv)).is_some() {
                break;
            }
            let x: Vec<f64> = bounds
                .iter()
                .map(|&(lo, hi)| rng.uniform_in(lo, hi))
                .collect();
            let fv = eval(&x);
            nfe += 1;
            if eps_less(fv, best_fv, 0.0) {
                best_fv = fv;
                best_x = x.clone();
            }
            pop.push(x);
            fit.push(fv);
        }
        let mut n = pop.len();

        // ε(0): configured, or the θ = 0.05 violation quantile of the initial
        // population.
        let eps0 = self.epsilon0.unwrap_or_else(|| {
            let mut viols: Vec<f64> = fit.iter().map(|&(_, v)| v).collect();
            viols.sort_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal));
            let idx = ((0.05 * viols.len() as f64) as usize).min(viols.len().saturating_sub(1));
            let q = viols.get(idx).copied().unwrap_or(0.0);
            if q.is_finite() {
                q
            } else {
                0.0
            }
        });
        let epsilon = |t: usize| -> f64 {
            if t >= t_c || t_c == 0 {
                0.0
            } else {
                eps0 * (1.0 - t as f64 / t_c as f64).powf(self.cp)
            }
        };

        let mut trial = vec![0.0; dim];
        'outer: while n >= N_MIN && term.reason(nfe, best_for_term(best_fv)).is_none() {
            let eps = epsilon(nfe);

            // Rank for pbest selection under the current ε.
            let mut ranked: Vec<usize> = (0..n).collect();
            ranked.sort_by(|&a, &b| eps_cmp(fit[a], fit[b], eps));
            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_for_term(best_fv)).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)
                };
                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();
                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);
                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;
                    } else {
                        trial[j] = pop[i][j];
                    }
                }

                let tfv = eval(&trial);
                nfe += 1;
                if eps_less(tfv, best_fv, 0.0) {
                    best_fv = tfv;
                    best_x.copy_from_slice(&trial);
                }
                // ε-selection: replace when not ε-worse; strict ε-improvement
                // feeds the archive and the success history.
                if !eps_less(fit[i], tfv, eps) {
                    if eps_less(tfv, fit[i], eps) {
                        succ_cr.push(cr);
                        succ_f.push(scale);
                        // Improvement magnitude: objective and violation gains
                        // both count (LSHADE44-style), capped to stay finite.
                        let df = (fit[i].0 - tfv.0).max(0.0);
                        let dv = (fit[i].1 - tfv.1).max(0.0);
                        let d = df + dv;
                        delta.push(if d.is_finite() {
                            d.max(1e-12)
                        } else {
                            f64::MAX
                        });
                        let arc_size = (self.archive_rate * n as f64).round() as usize;
                        if archive.len() < arc_size.max(1) {
                            archive.push(pop[i].clone());
                        } else if arc_size > 0 {
                            let idx = rng.index(archive.len());
                            archive[idx] = pop[i].clone();
                        }
                    }
                    pop[i].copy_from_slice(&trial);
                    fit[i] = tfv;
                }
            }

            let total: f64 = delta.iter().sum();
            if !succ_f.is_empty() && total.is_finite() && 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;
            }

            // LPSR, dropping the ε-worst individuals.
            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| eps_cmp(fit[a], fit[b], eps));
                ranked2.truncate(n_new);
                pop = ranked2
                    .iter()
                    .map(|&i| std::mem::take(&mut pop[i]))
                    .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_for_term(best_fv))
            .unwrap_or(StopReason::BudgetExhausted);
        Report {
            solution: Solution {
                x: best_x,
                // Report the objective when the best is feasible; otherwise
                // +inf (the decision vector still holds the least-violating
                // point found — re-check `total_violation` on it).
                value: best_for_term(best_fv),
            },
            stop,
            evaluations: nfe,
        }
    }

    /// Returns the same configuration with its RNG seed replaced.
    pub fn with_seed(&self, seed: u64) -> Self {
        EpsilonLShade { seed, ..*self }
    }
}

/// Weighted Lehmer mean (as in [`LShade`](super::LShade)).
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()
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::constraint::constrained_func;

    #[test]
    fn solves_a_constrained_quadratic() {
        // Minimize x² + y² subject to x + y ≥ 1: optimum (0.5, 0.5), f = 0.5.
        let p = constrained_func(
            vec![(-2.0, 2.0); 2],
            |x| x[0] * x[0] + x[1] * x[1],
            |x| vec![(1.0 - x[0] - x[1]).max(0.0)],
        );
        let r = EpsilonLShade::default().optimize(&p, &Termination::budget(8000));
        assert!(
            (r.best_value() - 0.5).abs() < 1e-4,
            "got {}",
            r.best_value()
        );
    }

    #[test]
    fn is_deterministic_and_budget_exact() {
        let p = constrained_func(
            vec![(-5.0, 5.0); 3],
            |x| x.iter().map(|v| v * v).sum(),
            |x| vec![(1.0 - x[0]).max(0.0)],
        );
        let t = Termination::budget(3000);
        let a = EpsilonLShade::default().optimize(&p, &t);
        let b = EpsilonLShade::default().optimize(&p, &t);
        assert_eq!(a.solution, b.solution);
        assert_eq!(a.evaluations, 3000);
    }

    #[test]
    fn equality_constraint_via_tolerance() {
        // Minimize x + y subject to x·y = 1 (tolerance 1e-4): optimum (1, 1).
        let p = constrained_func(
            vec![(0.01, 10.0); 2],
            |x| x[0] + x[1],
            |x| vec![((x[0] * x[1] - 1.0).abs() - 1e-4).max(0.0)],
        );
        let r = EpsilonLShade::default().optimize(&p, &Termination::budget(12_000));
        assert!(
            (r.best_value() - 2.0).abs() < 1e-2,
            "got {}",
            r.best_value()
        );
    }
}