metalforge 0.1.6

//! 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, from the worst
//! non-dominated front, the solution 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::{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 {
        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 = 1e30; // 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 {
            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 });
        }

        while 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: drop the least-contributing member of the worst
            // front (or the lone member if that front is a singleton).
            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 {
                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.
fn least_contributor(objs: &[&[f64]]) -> usize {
    let m = objs[0].len();
    let r = reference_point(objs);
    if m == 2 {
        least_contributor_2d(objs, &r)
    } else {
        least_contributor_general(objs, &r)
    }
}

/// Reference point dominated by every solution: per-objective max plus a small
/// range-scaled offset.
fn reference_point(objs: &[&[f64]]) -> Vec<f64> {
    let m = objs[0].len();
    (0..m)
        .map(|j| {
            let mut lo = f64::INFINITY;
            let mut hi = f64::NEG_INFINITY;
            for o in objs {
                lo = lo.min(o[j]);
                hi = hi.max(o[j]);
            }
            hi + (0.1 * (hi - lo)).max(1e-6)
        })
        .collect()
}

/// Bi-objective exact contributions via a sweep (`O(n log n)`): the exclusive
/// rectangle of each front point.
fn least_contributor_2d(objs: &[&[f64]], r: &[f64]) -> usize {
    let k = objs.len();
    let mut order: Vec<usize> = (0..k).collect();
    order.sort_by(|&a, &b| {
        objs[a][0]
            .partial_cmp(&objs[b][0])
            .unwrap_or(std::cmp::Ordering::Equal)
            .then(
                objs[a][1]
                    .partial_cmp(&objs[b][1])
                    .unwrap_or(std::cmp::Ordering::Equal),
            )
    });
    let mut best = 0usize;
    let mut best_contr = f64::INFINITY;
    for pos in 0..k {
        let i = order[pos];
        let xi = objs[i][0];
        let yi = objs[i][1];
        let x_next = if pos + 1 < k {
            objs[order[pos + 1]][0]
        } else {
            r[0]
        };
        let y_prev = if pos > 0 {
            objs[order[pos - 1]][1]
        } else {
            r[1]
        };
        let contr = (x_next - xi) * (y_prev - yi);
        if contr < best_contr {
            best_contr = contr;
            best = i;
        }
    }
    best
}

/// General-M contributions: `HV(F) − HV(F \ {p})` for each `p`, via WFG.
fn least_contributor_general(objs: &[&[f64]], r: &[f64]) -> usize {
    let full: Vec<Vec<f64>> = objs.iter().map(|o| o.to_vec()).collect();
    let total = hypervolume(&full, r);
    let mut best = 0usize;
    let mut best_contr = f64::INFINITY;
    for i in 0..objs.len() {
        let without: Vec<Vec<f64>> = (0..objs.len())
            .filter(|&j| j != i)
            .map(|j| objs[j].to_vec())
            .collect();
        let contr = total - hypervolume(&without, r);
        if contr < best_contr {
            best_contr = contr;
            best = i;
        }
    }
    best
}

/// Exact dominated hypervolume of `points` w.r.t. reference `r` (minimization;
/// every point assumed to dominate `r`), by the WFG exclusive decomposition.
fn hypervolume(points: &[Vec<f64>], r: &[f64]) -> f64 {
    if points.is_empty() {
        return 0.0;
    }
    let mut vol = 0.0;
    for i in 0..points.len() {
        let incl: f64 = r
            .iter()
            .zip(&points[i])
            .map(|(rj, pj)| (rj - pj).max(0.0))
            .product();
        // Limit set: the overlap of each later point with points[i].
        let limited: Vec<Vec<f64>> = points[i + 1..]
            .iter()
            .map(|q| (0..r.len()).map(|j| points[i][j].max(q[j])).collect())
            .collect();
        let nd = non_dominated(&limited);
        vol += incl - hypervolume(&nd, r);
    }
    vol
}

/// Keeps the non-dominated subset (minimization).
fn non_dominated(set: &[Vec<f64>]) -> Vec<Vec<f64>> {
    let mut out = Vec::new();
    for (i, p) in set.iter().enumerate() {
        let dominated = set.iter().enumerate().any(|(j, q)| {
            i != j && q.iter().zip(p).all(|(a, b)| a <= b) && q.iter().zip(p).any(|(a, b)| a < b)
        });
        if !dominated {
            out.push(p.clone());
        }
    }
    out
}

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

    #[test]
    fn hypervolume_known_2d_and_3d() {
        // Two non-dominated 2D points vs r=(3,4): union area = 3.
        let pts = vec![vec![1.0, 3.0], vec![2.0, 2.0]];
        let hv = hypervolume(&pts, &[3.0, 4.0]);
        assert!((hv - 3.0).abs() < 1e-9, "2D HV = {hv}");

        // Single 3D point: box volume.
        let one = vec![vec![1.0, 1.0, 1.0]];
        let hv3 = hypervolume(&one, &[2.0, 3.0, 4.0]);
        assert!((hv3 - (1.0 * 2.0 * 3.0)).abs() < 1e-9, "3D HV = {hv3}");
    }

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