metalforge 0.3.0

//! SMS-EMOA — S-Metric Selection Evolutionary Multi-objective Algorithm
//! (Beume, Naujoks & Emmerich 2007, EJOR 181, 1653–1669).
//!
//! A steady-state MOEA whose selection maximizes the **dominated hypervolume**:
//! each generation it breeds a single offspring, then discards one member of
//! the worst non-dominated front — by the paper's two-criterion Reduce: when
//! the population holds **several** fronts, the last-front member dominated by
//! the most points (cheap, and better-founded for dominated points); when it is
//! a **single** front, the member with the smallest **hypervolume
//! contribution**. Driving the population by the hypervolume indicator gives
//! tighter convergence and spread than dominance + crowding at low objective
//! counts — SMS-EMOA empirically beats NSGA-II / SPEA2 / ε-MOEA on 2-3
//! objectives.
//!
//! Hypervolume is exact: an `O(n log n)` sweep for the common **bi-objective**
//! case, and the WFG recursion (While et al. 2012) for 3+ objectives. The latter
//! is correct but cost grows with objectives and front size, so SMS-EMOA is
//! intended for **2-3 objectives** — use [`NsgaIII`](super::NsgaIII) for
//! many-objective problems.
//!
//! Returns a [`ParetoFront`]; deterministic for a given seed.

use super::nsga2::{dominates, fast_non_dominated_sort, mutate, sbx};
use crate::problem::MultiProblem;
use crate::rng::Rng;
use crate::solution::{MultiSolution, ParetoFront};
use crate::termination::Termination;

/// SMS-EMOA configuration.
#[derive(Debug, Clone, Copy)]
pub struct SmsEmoa {
    /// Population size.
    pub pop_size: usize,
    /// SBX distribution index `η_c`.
    pub crossover_eta: f64,
    /// Probability of applying SBX to the parent pair.
    pub crossover_prob: f64,
    /// Polynomial-mutation distribution index `η_m`.
    pub mutation_eta: f64,
    /// Per-variable mutation probability; `None` uses `1/dim`.
    pub mutation_prob: Option<f64>,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for SmsEmoa {
    fn default() -> Self {
        SmsEmoa {
            pop_size: 100,
            crossover_eta: 15.0,
            crossover_prob: 0.9,
            mutation_eta: 20.0,
            mutation_prob: None,
            seed: 42,
        }
    }
}

impl SmsEmoa {
    /// Approximates the Pareto front of `problem` within the evaluation budget
    /// of `term` (its `target` is ignored).
    pub fn optimize(&self, problem: &dyn MultiProblem, term: &Termination) -> ParetoFront {
        crate::problem::validate_multi(problem)
            .unwrap_or_else(|e| panic!("SmsEmoa: invalid problem: {e}"));
        let bounds = problem.bounds();
        let dim = bounds.len();
        let n = self.pop_size.max(2);
        let pm = self.mutation_prob.unwrap_or(1.0 / dim as f64);
        let mut rng = Rng::new(self.seed);

        let eval = |x: &[f64]| -> Vec<f64> {
            let mut o = problem.objectives(x);
            for v in &mut o {
                if !v.is_finite() {
                    *v = INFEASIBLE; // sanitized worst, keeps the HV geometry finite
                }
            }
            o
        };

        let mut pop: Vec<MultiSolution> = Vec::with_capacity(n + 1);
        let mut evaluations = 0;
        for _ in 0..n {
            if evaluations >= term.max_evaluations {
                break;
            }
            let x: Vec<f64> = bounds
                .iter()
                .map(|&(lo, hi)| rng.uniform_in(lo, hi))
                .collect();
            let objectives = eval(&x);
            evaluations += 1;
            pop.push(MultiSolution { x, objectives });
        }
        if pop.is_empty() {
            return ParetoFront {
                solutions: Vec::new(),
                evaluations,
            };
        }

        while pop.len() == n && evaluations < term.max_evaluations {
            // Breed one offspring (random parents; SBX keeps the first child).
            let a = rng.index(n);
            let b = rng.index(n);
            let (mut child, _) = sbx(
                &pop[a].x,
                &pop[b].x,
                bounds,
                self.crossover_eta,
                self.crossover_prob,
                &mut rng,
            );
            mutate(&mut child, bounds, self.mutation_eta, pm, &mut rng);
            let objectives = eval(&child);
            evaluations += 1;
            pop.push(MultiSolution {
                x: child,
                objectives,
            });

            // Reduce by one, per the paper's two-criterion rule.
            let fronts = fast_non_dominated_sort(&pop);
            let last = fronts.last().expect("at least one front");
            let drop = if last.len() == 1 {
                last[0]
            } else if let Some(&bad) = last
                .iter()
                .find(|&&i| pop[i].objectives.iter().any(|&v| v >= INFEASIBLE))
            {
                // A sanitized-infeasible point must never survive by inflating
                // the reference point; it goes first.
                bad
            } else if fronts.len() > 1 {
                // Several fronts: drop the last-front member dominated by the
                // most points (Beume et al. 2007, §3.1) — cheaper than HV and
                // better-founded for dominated points.
                *last
                    .iter()
                    .max_by_key(|&&i| {
                        pop.iter()
                            .filter(|q| dominates(&q.objectives, &pop[i].objectives))
                            .count()
                    })
                    .expect("last front is non-empty")
            } else {
                // Single front: drop the smallest hypervolume contributor.
                let objs: Vec<&[f64]> =
                    last.iter().map(|&i| pop[i].objectives.as_slice()).collect();
                last[least_contributor(&objs)]
            };
            pop.swap_remove(drop);
        }

        let fronts = fast_non_dominated_sort(&pop);
        let solutions = fronts[0].iter().map(|&i| pop[i].clone()).collect();
        ParetoFront {
            solutions,
            evaluations,
        }
    }
}

/// Index (into `objs`) of the solution with the smallest hypervolume
/// contribution, relative to a nadir-based reference point. Contributions
/// come from the shared [`crate::indicators::hv_contributions`] (2D sweep /
/// WFG); ties break to the lowest index deterministically.
fn least_contributor(objs: &[&[f64]]) -> usize {
    let r = reference_point(objs);
    let points: Vec<Vec<f64>> = objs.iter().map(|o| o.to_vec()).collect();
    let contr = crate::indicators::hv_contributions(&points, &r);
    contr
        .iter()
        .enumerate()
        .min_by(|a, b| a.1.partial_cmp(b.1).unwrap_or(std::cmp::Ordering::Equal))
        .map(|(i, _)| i)
        .unwrap_or(0)
}

/// Sanitized objective value standing in for a non-finite (infeasible) result.
const INFEASIBLE: f64 = 1e30;

/// Reference point dominated by every solution: per-objective max plus the
/// paper's offset of 1 (Beume et al. 2007), which protects the exclusive
/// volume of boundary solutions better than a range-scaled sliver.
fn reference_point(objs: &[&[f64]]) -> Vec<f64> {
    let m = objs[0].len();
    (0..m)
        .map(|j| {
            let hi = objs.iter().map(|o| o[j]).fold(f64::NEG_INFINITY, f64::max);
            hi + 1.0
        })
        .collect()
}

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

    #[test]
    fn least_contributor_drops_the_redundant_point() {
        // Three points; the middle one is nearly redundant (tiny exclusive area).
        let a = [0.0, 2.0];
        let b = [1.0, 1.99];
        let c = [2.0, 0.0];
        let objs: Vec<&[f64]> = vec![&a, &b, &c];
        assert_eq!(least_contributor(&objs), 1);
    }
}