metalforge 0.3.0

//! Quality indicators for multi-objective fronts.
//!
//! The standard toolkit for measuring how good an approximation front is
//! (minimization convention throughout, matching [`MultiProblem`]):
//!
//! - [`hypervolume`] — the dominated volume w.r.t. a reference point
//!   (Zitzler & Thiele 1999): the only strictly Pareto-compliant unary
//!   indicator. Exact `O(n log n)` sweep in 2D; WFG exclusive decomposition
//!   (While, Bradstreet & Barone 2012) for 3+ objectives.
//! - [`hv_contributions`] — each point's exclusive hypervolume (what
//!   [`SmsEmoa`] selects on).
//! - [`gd`] / [`igd`] — (inverted) generational distance to a reference
//!   front (Van Veldhuizen 1999; Coello & Reyes-Sierra 2004).
//! - [`gd_plus`] / [`igd_plus`] — the weakly-Pareto-compliant `+` variants
//!   (Ishibuchi, Masuda, Tanigaki & Nojima, EMO 2015), which replace the
//!   Euclidean distance with the *dominance-aware* distance
//!   `d⁺(z, a) = ‖max(a − z, 0)‖` so that dominating a reference point never
//!   reads as an error. Prefer IGD+ over plain IGD when comparing algorithms.
//!
//! All functions are pure, deterministic, and dependency-free. Reference
//! fronts are ordinary `&[Vec<f64>]` slices, so analytical fronts, sampled
//! fronts, or another run's [`ParetoFront::objective_vectors`] all work.
//!
//! [`MultiProblem`]: crate::problem::MultiProblem
//! [`SmsEmoa`]: crate::algo::SmsEmoa
//! [`ParetoFront::objective_vectors`]: crate::solution::ParetoFront::objective_vectors

/// Exact dominated hypervolume of `points` w.r.t. the reference point `r`
/// (minimization: every point should weakly dominate `r`; components beyond
/// `r` contribute zero). Duplicates and dominated points are handled
/// correctly (they simply add no exclusive volume).
///
/// 2D uses an `O(n log n)` sweep; 3+ objectives use the WFG exclusive
/// decomposition — exact but exponential in the worst case, intended for the
/// front sizes typical of population-based MOEAs (≤ a few hundred points,
/// 2–4 objectives).
pub fn hypervolume(points: &[Vec<f64>], r: &[f64]) -> f64 {
    if points.is_empty() {
        return 0.0;
    }
    if r.len() == 2 {
        return hypervolume_2d(points, r);
    }
    wfg(points, r)
}

/// The exclusive hypervolume contribution of every point: `contribution[i] =
/// HV(points) − HV(points \ {i})`. Dominated or duplicate points get 0.
///
/// 2D is computed in one `O(n log n)` sweep; 3+ objectives fall back to
/// leave-one-out WFG (`n + 1` hypervolume evaluations).
pub fn hv_contributions(points: &[Vec<f64>], r: &[f64]) -> Vec<f64> {
    if points.is_empty() {
        return Vec::new();
    }
    if r.len() == 2 {
        return contributions_2d(points, r);
    }
    let total = wfg(points, r);
    (0..points.len())
        .map(|i| {
            let without: Vec<Vec<f64>> = points
                .iter()
                .enumerate()
                .filter(|&(j, _)| j != i)
                .map(|(_, p)| p.clone())
                .collect();
            (total - wfg(&without, r)).max(0.0)
        })
        .collect()
}

/// Generational Distance: the mean Euclidean distance from each point of
/// `front` to its nearest neighbor in `reference` (Van Veldhuizen 1999).
/// Measures convergence only; `NaN` for an empty front.
pub fn gd(front: &[Vec<f64>], reference: &[Vec<f64>]) -> f64 {
    mean_min_distance(front, reference, euclidean)
}

/// Inverted Generational Distance: the mean Euclidean distance from each
/// *reference* point to its nearest neighbor in `front`. Measures both
/// convergence and coverage; `NaN` for an empty reference set. Sensitive to
/// the reference-front resolution — prefer [`igd_plus`] for algorithm
/// comparisons.
pub fn igd(front: &[Vec<f64>], reference: &[Vec<f64>]) -> f64 {
    mean_min_distance(reference, front, euclidean)
}

/// GD⁺ (Ishibuchi et al. 2015): like [`gd`] but with the dominance-aware
/// distance `‖max(a − z, 0)‖` from front point `a` to reference point `z`,
/// so components where the front point is *better* than the reference do not
/// count as error.
pub fn gd_plus(front: &[Vec<f64>], reference: &[Vec<f64>]) -> f64 {
    mean_min_distance(front, reference, d_plus)
}

/// IGD⁺ (Ishibuchi et al. 2015): like [`igd`] but with the dominance-aware
/// distance `‖max(a − z, 0)‖` from each reference point `z` to front point
/// `a`. Weakly Pareto-compliant, unlike plain IGD — the recommended default
/// for comparing fronts against a reference.
pub fn igd_plus(front: &[Vec<f64>], reference: &[Vec<f64>]) -> f64 {
    mean_min_distance(reference, front, |z, a| d_plus(a, z))
}

/// Mean over `from` of the minimum `dist(p, q)` over `to`.
fn mean_min_distance(
    from: &[Vec<f64>],
    to: &[Vec<f64>],
    dist: impl Fn(&[f64], &[f64]) -> f64,
) -> f64 {
    if from.is_empty() || to.is_empty() {
        return f64::NAN;
    }
    let total: f64 = from
        .iter()
        .map(|p| to.iter().map(|q| dist(p, q)).fold(f64::INFINITY, f64::min))
        .sum();
    total / from.len() as f64
}

fn euclidean(a: &[f64], b: &[f64]) -> f64 {
    a.iter()
        .zip(b)
        .map(|(x, y)| (x - y) * (x - y))
        .sum::<f64>()
        .sqrt()
}

/// Dominance-aware distance from front point `a` to reference point `z`
/// (minimization): `‖max(a − z, 0)‖`.
fn d_plus(a: &[f64], z: &[f64]) -> f64 {
    a.iter()
        .zip(z)
        .map(|(ai, zi)| (ai - zi).max(0.0).powi(2))
        .sum::<f64>()
        .sqrt()
}

/// 2D hypervolume by a sweep over the non-dominated staircase.
fn hypervolume_2d(points: &[Vec<f64>], r: &[f64]) -> f64 {
    // Sort by f1 ascending, f2 ascending; walk keeping the running best f2.
    let mut sorted: Vec<&Vec<f64>> = points.iter().collect();
    sorted.sort_by(|a, b| {
        a[0].partial_cmp(&b[0])
            .unwrap_or(std::cmp::Ordering::Equal)
            .then(a[1].partial_cmp(&b[1]).unwrap_or(std::cmp::Ordering::Equal))
    });
    let mut vol = 0.0;
    let mut best_f2 = r[1];
    for p in sorted {
        if p[1] < best_f2 && p[0] < r[0] {
            vol += (r[0] - p[0]) * (best_f2 - p[1]);
            best_f2 = p[1];
        }
    }
    vol
}

/// 2D exclusive contributions in one sweep over the sorted staircase.
fn contributions_2d(points: &[Vec<f64>], r: &[f64]) -> Vec<f64> {
    let k = points.len();
    let mut order: Vec<usize> = (0..k).collect();
    order.sort_by(|&a, &b| {
        points[a][0]
            .partial_cmp(&points[b][0])
            .unwrap_or(std::cmp::Ordering::Equal)
            .then(
                points[a][1]
                    .partial_cmp(&points[b][1])
                    .unwrap_or(std::cmp::Ordering::Equal),
            )
    });
    // Keep only the staircase (strictly improving f2 as f1 grows); everything
    // else is dominated or duplicate and contributes 0.
    let mut contr = vec![0.0; k];
    let mut stair: Vec<usize> = Vec::with_capacity(k);
    let mut best_f2 = f64::INFINITY;
    for &i in &order {
        if points[i][1] < best_f2 {
            stair.push(i);
            best_f2 = points[i][1];
        }
    }
    for (pos, &i) in stair.iter().enumerate() {
        let x_next = if pos + 1 < stair.len() {
            points[stair[pos + 1]][0]
        } else {
            r[0]
        };
        let y_prev = if pos > 0 {
            points[stair[pos - 1]][1]
        } else {
            r[1]
        };
        contr[i] = ((x_next - points[i][0]) * (y_prev - points[i][1])).max(0.0);
    }
    contr
}

/// WFG exclusive-decomposition hypervolume (While et al. 2012):
/// `HV(S) = Σᵢ [incl(pᵢ) − HV(nd(limit(pᵢ, p_{i+1..})))]`.
fn wfg(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 - wfg(&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]];
        assert!((hypervolume(&pts, &[3.0, 4.0]) - 3.0).abs() < 1e-9);
        // A dominated third point must not change the volume.
        let with_dominated = vec![vec![1.0, 3.0], vec![2.0, 2.0], vec![2.5, 3.5]];
        assert!((hypervolume(&with_dominated, &[3.0, 4.0]) - 3.0).abs() < 1e-9);
        // Single 3D point: box volume.
        let one = vec![vec![1.0, 1.0, 1.0]];
        assert!((hypervolume(&one, &[2.0, 3.0, 4.0]) - 6.0).abs() < 1e-9);
    }

    #[test]
    fn wfg_and_2d_sweep_agree() {
        // Same 2D front evaluated by both paths must match exactly.
        let pts = vec![
            vec![0.1, 0.9],
            vec![0.4, 0.5],
            vec![0.4, 0.5], // duplicate
            vec![0.8, 0.2],
            vec![0.9, 0.8], // dominated
        ];
        let r = [1.1, 1.1];
        let sweep = hypervolume_2d(&pts, &r);
        let recursive = wfg(&pts, &r);
        assert!(
            (sweep - recursive).abs() < 1e-12,
            "sweep {sweep} vs wfg {recursive}"
        );
    }

    #[test]
    fn contributions_match_leave_one_out() {
        let pts = vec![vec![0.1, 0.9], vec![0.4, 0.5], vec![0.8, 0.2]];
        let r = [1.0, 1.0];
        let fast = hv_contributions(&pts, &r);
        let total = hypervolume(&pts, &r);
        for (i, f) in fast.iter().enumerate() {
            let without: Vec<Vec<f64>> = pts
                .iter()
                .enumerate()
                .filter(|&(j, _)| j != i)
                .map(|(_, p)| p.clone())
                .collect();
            let slow = total - hypervolume(&without, &r);
            assert!(
                (f - slow).abs() < 1e-12,
                "point {i}: fast {f} vs slow {slow}"
            );
        }
        // Duplicates and dominated points contribute zero.
        let with_extra = vec![
            vec![0.1, 0.9],
            vec![0.1, 0.9],
            vec![0.5, 0.95], // dominated
        ];
        let c = hv_contributions(&with_extra, &r);
        assert_eq!(c[1], 0.0);
        assert_eq!(c[2], 0.0);
    }

    #[test]
    fn igd_plus_ignores_dominating_deviation() {
        // A front point that DOMINATES the reference point: plain IGD reads a
        // positive error, IGD+ correctly reads zero.
        let reference = vec![vec![0.5, 0.5]];
        let dominating_front = vec![vec![0.4, 0.4]];
        assert!(igd(&dominating_front, &reference) > 0.0);
        assert_eq!(igd_plus(&dominating_front, &reference), 0.0);
        // A front point dominated BY the reference reads the same in both.
        let dominated_front = vec![vec![0.6, 0.6]];
        let d_igd = igd(&dominated_front, &reference);
        let d_plusv = igd_plus(&dominated_front, &reference);
        assert!((d_igd - d_plusv).abs() < 1e-12);
    }

    #[test]
    fn gd_igd_basics() {
        let front = vec![vec![0.0, 1.0], vec![1.0, 0.0]];
        // Front == reference: all indicators are zero.
        assert_eq!(gd(&front, &front), 0.0);
        assert_eq!(igd(&front, &front), 0.0);
        assert_eq!(gd_plus(&front, &front), 0.0);
        assert_eq!(igd_plus(&front, &front), 0.0);
        // Empty inputs give NaN, not a panic.
        assert!(gd(&[], &front).is_nan());
        assert!(igd(&front, &[]).is_nan());
    }
}