metalforge 0.3.0

//! Shuffled Complex Evolution — SCE-UA (Duan, Sorooshian & Gupta 1992,
//! Water Resources Research 28, 1015–1031).
//!
//! A robust global optimizer that partitions a random population into
//! *complexes*, evolves each by Competitive Complex Evolution (a
//! reflection/contraction simplex with triangular parent selection), then
//! periodically shuffles the complexes back together. More expensive than DDS
//! but stronger on multimodal surfaces. Migrated from the implementation
//! validated in `rainflow-core`; minimizing here (the paper's sense), whereas
//! rainflow wrapped it to maximize NSE/KGE.

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

/// SCE-UA configuration. The complex geometry (points per complex, sub-complex
/// size, evolution steps) follows the standard `m = 2n+1`, `q = n+1`, `β = m`
/// of Duan et al. 1992, so only the number of complexes is exposed.
#[derive(Debug, Clone, Copy)]
pub struct Sce {
    /// Number of complexes `p` (≥ 1). More complexes ⇒ more global, slower.
    pub complexes: usize,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for Sce {
    fn default() -> Self {
        Sce {
            complexes: 4,
            seed: 42,
        }
    }
}

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

    /// Minimizes `problem` within its bounds using SCE-UA.
    ///
    /// Stops on the [`Termination`] budget/target only; the paper's own
    /// convergence criteria (`kstop`/`pcento`/`peps`) are not implemented.
    fn optimize(&self, problem: &dyn Problem, term: &Termination) -> Report {
        crate::problem::validate(problem).unwrap_or_else(|e| panic!("Sce: invalid problem: {e}"));
        let bounds = problem.bounds();
        let n = bounds.len();
        let complexes = self.complexes.max(1);

        // Standard SCE-UA complex geometry.
        let m = 2 * n + 1; // points per complex
        let q = n + 1; // points per sub-complex (simplex)
        let beta = m; // evolution steps per complex per shuffle
        let pop_size = complexes * m;

        let mut rng = Rng::new(self.seed);

        // Seed the evaluator with the first sampled point, then fill the rest.
        let first_x = sample(bounds, &mut rng);
        let first_v = problem.objective(&first_x);
        let mut ev = Evaluator::new(
            problem,
            term,
            Solution {
                x: first_x.clone(),
                value: finite_or_worst(first_v),
            },
        );

        let mut pop: Vec<(Vec<f64>, f64)> = Vec::with_capacity(pop_size);
        pop.push((first_x, finite_or_worst(first_v)));
        for _ in 1..pop_size {
            if ev.done() {
                break;
            }
            let x = sample(bounds, &mut rng);
            let s = finite_or_worst(ev.eval(&x));
            pop.push((x, s));
        }

        let cmp = |a: &(Vec<f64>, f64), b: &(Vec<f64>, f64)| {
            a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal)
        };

        while !ev.done() && pop.len() == pop_size {
            pop.sort_by(cmp);

            // Partition the rank-sorted population into complexes round-robin.
            for k in 0..complexes {
                let mut complex: Vec<(Vec<f64>, f64)> =
                    (0..m).map(|j| pop[k + j * complexes].clone()).collect();

                // Competitive Complex Evolution: evolve `beta` times.
                for _ in 0..beta {
                    if ev.done() {
                        break;
                    }
                    complex.sort_by(cmp);
                    // Triangular selection biases parents toward better points.
                    let parents = triangular_subset(m, q, &mut rng);
                    let worst_idx = *parents.last().unwrap();

                    // Centroid of the q-1 best parents (excluding the worst).
                    let mut centroid = vec![0.0; n];
                    for &pi in &parents[..q - 1] {
                        for (c, &v) in centroid.iter_mut().zip(&complex[pi].0) {
                            *c += v;
                        }
                    }
                    let qm1 = (q - 1) as f64;
                    for c in &mut centroid {
                        *c /= qm1;
                    }

                    // Reflect the worst point through the centroid.
                    let worst = complex[worst_idx].0.clone();
                    let mut trial: Vec<f64> = (0..n)
                        .map(|d| centroid[d] + (centroid[d] - worst[d]))
                        .collect();
                    let in_bounds = trial
                        .iter()
                        .zip(bounds)
                        .all(|(&v, &(lo, hi))| v >= lo && v <= hi);

                    // One CCE step can cost up to 3 evaluations; the budget is
                    // re-checked before each so the contract is never exceeded.
                    let mut tf;
                    if in_bounds {
                        tf = finite_or_worst(ev.eval(&trial));
                        // Ties are accepted (paper: replace when not worse).
                        if tf > complex[worst_idx].1 {
                            if ev.done() {
                                break;
                            }
                            // Reflection failed: contract toward the centroid.
                            trial = (0..n).map(|d| 0.5 * (centroid[d] + worst[d])).collect();
                            tf = finite_or_worst(ev.eval(&trial));
                            if tf > complex[worst_idx].1 {
                                if ev.done() {
                                    break;
                                }
                                // Contraction failed too: random mutation in hull.
                                trial = sample_in_hull(&complex, bounds, &mut rng);
                                tf = finite_or_worst(ev.eval(&trial));
                            }
                        }
                    } else {
                        // Out of bounds: random mutation in the complex's hull.
                        trial = sample_in_hull(&complex, bounds, &mut rng);
                        tf = finite_or_worst(ev.eval(&trial));
                    }

                    complex[worst_idx] = (trial, tf);
                }

                // Shuffle the evolved complex back into the population.
                for (j, point) in complex.into_iter().enumerate() {
                    pop[k + j * complexes] = point;
                }
            }
        }

        ev.finish()
    }
}

/// Non-finite scores are treated as the worst possible (`+∞`), so degenerate
/// parameter combinations are rejected rather than poisoning the population.
#[inline]
fn finite_or_worst(v: f64) -> f64 {
    if v.is_finite() {
        v
    } else {
        f64::INFINITY
    }
}

/// Picks `q` distinct indices from `0..m` with a triangular probability that
/// favors lower (better, since the complex is rank-sorted) indices, returned
/// sorted ascending (Duan et al. 1992, trapezoidal distribution).
fn triangular_subset(m: usize, q: usize, rng: &mut Rng) -> Vec<usize> {
    let mut chosen = Vec::with_capacity(q);
    let mf = m as f64;
    while chosen.len() < q {
        // P(i) = 2(m - i)/(m(m+1)); inverse-CDF sampling on rank.
        let u = rng.uniform();
        let idx = ((mf + 0.5) - ((mf + 0.5).powi(2) - mf * (mf + 1.0) * u).sqrt()) as usize;
        let idx = idx.min(m - 1);
        if !chosen.contains(&idx) {
            chosen.push(idx);
        }
    }
    chosen.sort_unstable();
    chosen
}

/// Uniform sample inside the smallest axis-aligned box containing the complex,
/// clamped to the global bounds (SCE-UA's mutation step).
fn sample_in_hull(complex: &[(Vec<f64>, f64)], bounds: &[Bound], rng: &mut Rng) -> Vec<f64> {
    let n = bounds.len();
    (0..n)
        .map(|d| {
            let mut lo = complex[0].0[d];
            let mut hi = lo;
            for point in complex {
                lo = lo.min(point.0[d]);
                hi = hi.max(point.0[d]);
            }
            let v = rng.uniform_in(lo, hi);
            clamp(v, bounds[d].0, bounds[d].1)
        })
        .collect()
}