pare 0.2.1

//! Pareto frontier and skyline query primitives for multi-objective optimization.
//!
//! Maintains a set of non-dominated points across multiple dimensions,
//! supporting maximization/minimization per dimension and crowding distance
//! for diversity maintenance.

#![warn(missing_docs)]

use std::cmp::Ordering;

use thiserror::Error;

pub mod sensitivity;

/// Direction of optimization for a dimension.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum Direction {
    /// Higher values are better.
    Maximize,
    /// Lower values are better.
    Minimize,
}

/// A point in multi-objective space.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Point<V> {
    /// The actual values across all objectives.
    pub values: Vec<f64>,
    /// Associated data (e.g. an ID or candidate metadata).
    pub data: V,
}

/// Statistics for normalizing objective values.
#[derive(Debug, Clone, Copy, Default)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct NormalizationStats {
    /// Population mean of the values.
    pub mean: f64,
    /// Population standard deviation (divides by `n`, not `n-1`).
    pub std: f64,
    /// Minimum observed value.
    pub min: f64,
    /// Maximum observed value.
    pub max: f64,
}

impl NormalizationStats {
    /// Compute normalization statistics from a slice of values.
    pub fn from_values(values: &[f64]) -> Self {
        if values.is_empty() {
            return Self::default();
        }
        let n = values.len() as f64;
        let mean = values.iter().sum::<f64>() / n;
        let var = values.iter().map(|v| (v - mean) * (v - mean)).sum::<f64>() / n;
        let std = var.sqrt();
        let min = values.iter().copied().fold(f64::INFINITY, f64::min);
        let max = values.iter().copied().fold(f64::NEG_INFINITY, f64::max);
        Self {
            mean,
            std,
            min,
            max,
        }
    }
}

/// Pareto frontier maintaining a set of non-dominated points.
#[derive(Debug, Clone)]
pub struct ParetoFrontier<V> {
    points: Vec<Point<V>>,
    directions: Vec<Direction>,
    stats: Vec<NormalizationStats>,
    labels: Vec<String>,
    eps: f64,
}

/// Errors returned when constructing a frontier from raw vectors.
#[derive(Debug, Clone, PartialEq, Eq, Error)]
pub enum FrontierError {
    /// No points were provided.
    #[error("no points provided")]
    Empty,
    /// Points had inconsistent dimensionality (or dimension was 0).
    #[error("points must be non-empty and have consistent dimensionality")]
    InconsistentDimensions,
    /// A value was NaN or infinite.
    #[error("point value at [{point_idx}][{dim_idx}] is not finite")]
    NonFinite {
        /// Index of the point containing the non-finite value.
        point_idx: usize,
        /// Index of the dimension containing the non-finite value.
        dim_idx: usize,
    },
}

impl<V> ParetoFrontier<V> {
    /// Create a new frontier with specified directions for each objective.
    pub fn new(directions: Vec<Direction>) -> Self {
        let dim = directions.len();
        Self {
            points: Vec::new(),
            directions,
            stats: vec![NormalizationStats::default(); dim],
            labels: (0..dim).map(|i| i.to_string()).collect(),
            eps: 1e-9,
        }
    }

    /// Set labels for objectives (e.g. `["accuracy", "latency"]`).
    ///
    /// Labels are stored for caller use and can be retrieved with [`labels`](Self::labels).
    /// They are not used internally by any scoring or dominance method.
    ///
    /// # Panics
    ///
    /// Panics if `labels.len() != directions.len()`.
    pub fn with_labels(mut self, labels: Vec<String>) -> Self {
        assert_eq!(
            labels.len(),
            self.directions.len(),
            "labels len ({}) must match directions len ({})",
            labels.len(),
            self.directions.len(),
        );
        self.labels = labels;
        self
    }

    /// Get the objective labels.
    ///
    /// Returns the labels set by [`with_labels`](Self::with_labels), or default
    /// numeric labels `["0", "1", ...]` if none were set.
    pub fn labels(&self) -> &[String] {
        &self.labels
    }

    /// Get the optimization directions.
    pub fn directions(&self) -> &[Direction] {
        &self.directions
    }

    /// Set the epsilon tolerance for dominance comparisons.
    ///
    /// The epsilon controls how strict dominance checks are:
    ///
    /// - Point `a` is considered "strictly better" than `b` in a dimension only
    ///   if `a > b + eps` (for Maximize) or `a + eps < b` (for Minimize).
    /// - Point `a` is considered "not worse" than `b` only if `a + eps >= b`
    ///   (for Maximize) or `a <= b + eps` (for Minimize).
    ///
    /// **Larger epsilon** = more tolerance = fewer domination events = larger frontiers.
    /// **Smaller epsilon** = less tolerance = more domination = smaller frontiers.
    ///
    /// The default is `1e-9`, which provides floating-point stability without
    /// materially changing results for objectives with reasonable magnitude.
    /// Set to `0.0` for exact comparisons.
    ///
    /// This is a *numerical tolerance*, not epsilon-dominance archiving (which
    /// deliberately coarsens objective space to bound archive size).
    pub fn with_eps(mut self, eps: f64) -> Self {
        self.eps = eps;
        self
    }

    /// Get the current epsilon tolerance.
    pub fn eps(&self) -> f64 {
        self.eps
    }

    /// Get the normalization statistics for each objective.
    pub fn stats(&self) -> &[NormalizationStats] {
        &self.stats
    }

    /// Number of points currently on the frontier.
    pub fn len(&self) -> usize {
        self.points.len()
    }

    /// Is the frontier empty?
    pub fn is_empty(&self) -> bool {
        self.points.is_empty()
    }

    /// Get all points on the frontier.
    pub fn points(&self) -> &[Point<V>] {
        &self.points
    }

    /// Mutable access to points.
    ///
    /// **Invariant warning:** callers may modify point data/values but must not introduce
    /// dominated points. Structural mutations (push/remove) that could violate the
    /// non-dominance invariant should go through [`push`](Self::push) instead.
    pub fn points_mut(&mut self) -> &mut [Point<V>] {
        &mut self.points
    }

    /// Add a point to the frontier if it is non-dominated.
    ///
    /// Returns `true` if the point was added (non-dominated), `false` if it was
    /// dominated by an existing point.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Maximize]);
    /// assert!(f.push(vec![0.9, 0.1], "a"));   // added
    /// assert!(f.push(vec![0.1, 0.9], "b"));   // added (trade-off)
    /// assert!(!f.push(vec![0.5, 0.05], "c"));  // dominated by "a"
    /// assert_eq!(f.len(), 2);
    /// ```
    ///
    /// # Panics
    ///
    /// Panics if `values.len() != directions.len()` or any value is NaN/infinite.
    pub fn push(&mut self, values: Vec<f64>, data: V) -> bool {
        assert_eq!(
            values.len(),
            self.directions.len(),
            "values len ({}) must match directions len ({})",
            values.len(),
            self.directions.len(),
        );
        assert!(
            values.iter().all(|v| v.is_finite()),
            "all values must be finite (no NaN or Inf)",
        );

        // Check if any existing point dominates the new point.
        for existing in &self.points {
            if dominates(&self.directions, self.eps, &existing.values, &values) {
                return false;
            }
        }

        // Remove points dominated by the new point.
        let directions = &self.directions;
        let eps = self.eps;
        self.points
            .retain(|existing| !dominates(directions, eps, &values, &existing.values));

        self.points.push(Point { values, data });
        self.update_stats();
        true
    }

    fn update_stats(&mut self) {
        for i in 0..self.directions.len() {
            let dim_values: Vec<f64> = self.points.iter().map(|p| p.values[i]).collect();
            self.stats[i] = NormalizationStats::from_values(&dim_values);
        }
    }

    /// Returns the normalized score for a point across all dimensions (higher is better).
    ///
    /// **Warning: dynamic normalization.** This method normalizes each dimension to `[0,1]`
    /// using the min/max of the *current frontier*. Adding or removing points changes the
    /// normalization, which means scores are **not stable across different frontier snapshots**.
    /// If you need stable scores, normalize objectives against fixed (static) bounds before
    /// insertion and use uniform weights on the pre-normalized values.
    ///
    /// # Panics
    ///
    /// Panics if `point_idx >= len()` or `weights.len() > directions.len()`.
    pub fn scalar_score(&self, point_idx: usize, weights: &[f64]) -> f64 {
        assert!(
            point_idx < self.points.len(),
            "point_idx ({point_idx}) out of bounds (len={})",
            self.points.len(),
        );
        assert!(
            weights.len() <= self.directions.len(),
            "weights len ({}) exceeds directions len ({})",
            weights.len(),
            self.directions.len(),
        );
        let p = &self.points[point_idx];
        let mut score = 0.0;
        let mut w_sum = 0.0;

        for (i, &w) in weights.iter().enumerate() {
            let val = p.values[i];
            let stat = self.stats[i];
            let range = stat.max - stat.min;
            let norm = if range > self.eps {
                (val - stat.min) / range
            } else {
                0.5
            };

            let oriented = match self.directions[i] {
                Direction::Maximize => norm,
                Direction::Minimize => 1.0 - norm,
            };

            score += oriented * w;
            w_sum += w;
        }

        if w_sum > 0.0 {
            score / w_sum
        } else {
            0.0
        }
    }

    /// Like [`scalar_score`][Self::scalar_score], but normalizes against caller-supplied bounds
    /// instead of the frontier's current min/max.
    ///
    /// This avoids the instability where adding a new arm to the frontier changes the
    /// normalized score of every existing arm.  Callers supply `bounds[i] = (min_i, max_i)`;
    /// values outside the range are clamped.  If `min == max` for a dimension, that
    /// dimension contributes `0` to the score.
    ///
    /// # Example
    ///
    /// ```rust
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Minimize]);
    /// f.push(vec![0.9, 10.0], "a");
    /// f.push(vec![0.7, 5.0],  "b");
    ///
    /// // ok_rate ∈ [0,1], cost ∈ [0,100]
    /// let bounds = vec![(0.0_f64, 1.0_f64), (0.0_f64, 100.0_f64)];
    /// let s0 = f.scalar_score_static(0, &[1.0, 1.0], &bounds);
    /// let s1 = f.scalar_score_static(1, &[1.0, 1.0], &bounds);
    /// assert!(s0 > s1, "arm a should score higher: better ok_rate, lower cost");
    /// ```
    pub fn scalar_score_static(
        &self,
        point_idx: usize,
        weights: &[f64],
        bounds: &[(f64, f64)],
    ) -> f64 {
        assert!(
            point_idx < self.points.len(),
            "point_idx ({point_idx}) out of bounds (len={})",
            self.points.len(),
        );
        let p = &self.points[point_idx];
        let dims = p.values.len().min(weights.len()).min(bounds.len());
        let mut score = 0.0;
        let mut w_sum = 0.0;

        for i in 0..dims {
            let (lo, hi) = bounds[i];
            let w = weights[i];
            let range = hi - lo;
            let norm = if range.abs() > self.eps {
                ((p.values[i] - lo) / range).clamp(0.0, 1.0)
            } else {
                0.0
            };
            let oriented = match self.directions[i] {
                Direction::Maximize => norm,
                Direction::Minimize => 1.0 - norm,
            };
            score += oriented * w;
            w_sum += w;
        }

        if w_sum > 0.0 {
            score / w_sum
        } else {
            0.0
        }
    }

    /// Find the index of the point with the highest weighted scalar score.
    pub fn best_index(&self, weights: &[f64]) -> Option<usize> {
        if self.is_empty() {
            return None;
        }
        (0..self.points.len()).max_by(|&a, &b| {
            self.scalar_score(a, weights)
                .partial_cmp(&self.scalar_score(b, weights))
                .unwrap_or(Ordering::Equal)
        })
    }

    /// Check if point A dominates point B.
    pub fn point_dominates(&self, a: &[f64], b: &[f64]) -> bool {
        dominates(&self.directions, self.eps, a, b)
    }

    /// Compute crowding distances for all points on the frontier (NSGA-II style).
    ///
    /// Boundary points get `f64::INFINITY`. Interior points get a finite score
    /// proportional to the spacing around them.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Maximize]);
    /// f.push(vec![0.0, 1.0], "a");
    /// f.push(vec![0.5, 0.5], "b");
    /// f.push(vec![1.0, 0.0], "c");
    ///
    /// let cd = f.crowding_distances();
    /// assert!(cd[0].is_infinite()); // boundary
    /// assert!(cd[2].is_infinite()); // boundary
    /// assert!(cd[1].is_finite());   // interior
    /// ```
    pub fn crowding_distances(&self) -> Vec<f64> {
        let n = self.points.len();
        if n == 0 {
            return Vec::new();
        }
        if n == 1 {
            return vec![f64::INFINITY];
        }

        let mut distances = vec![0.0; n];
        let mut indices: Vec<usize> = (0..n).collect();

        for i in 0..self.directions.len() {
            indices.sort_by(|&a, &b| {
                let av = self.points[a].values[i];
                let bv = self.points[b].values[i];
                av.partial_cmp(&bv).unwrap_or(Ordering::Equal)
            });

            let min_val = self.points[indices[0]].values[i];
            let max_val = self.points[indices[n - 1]].values[i];
            let range = max_val - min_val;

            // Fortin-Parizeau (GECCO 2013): all points tied at the boundary
            // value get infinite distance, not just the first/last index.
            let mut lo = 0;
            while lo < n && (self.points[indices[lo]].values[i] - min_val).abs() <= self.eps {
                distances[indices[lo]] = f64::INFINITY;
                lo += 1;
            }
            let mut hi = n - 1;
            while hi > 0 && (self.points[indices[hi]].values[i] - max_val).abs() <= self.eps {
                distances[indices[hi]] = f64::INFINITY;
                hi -= 1;
            }

            if range > self.eps {
                for window in indices.windows(3) {
                    let prev = window[0];
                    let curr = window[1];
                    let next = window[2];
                    if distances[curr].is_infinite() {
                        continue;
                    }
                    let val_prev = self.points[prev].values[i];
                    let val_next = self.points[next].values[i];
                    distances[curr] += (val_next - val_prev) / range;
                }
            }
        }
        distances
    }

    /// The ideal point: best observed value per objective on the frontier.
    ///
    /// For Maximize objectives, this is the maximum; for Minimize, the minimum.
    /// Returns `None` if the frontier is empty.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Minimize]);
    /// f.push(vec![0.9, 10.0], "a");
    /// f.push(vec![0.7, 5.0],  "b");
    ///
    /// let ideal = f.ideal_point().unwrap();
    /// assert!((ideal[0] - 0.9).abs() < 1e-9);  // best accuracy
    /// assert!((ideal[1] - 5.0).abs() < 1e-9);   // best (lowest) latency
    /// ```
    pub fn ideal_point(&self) -> Option<Vec<f64>> {
        if self.is_empty() {
            return None;
        }
        Some(
            (0..self.directions.len())
                .map(|i| {
                    let vals = self.points.iter().map(|p| p.values[i]);
                    match self.directions[i] {
                        Direction::Maximize => vals.fold(f64::NEG_INFINITY, f64::max),
                        Direction::Minimize => vals.fold(f64::INFINITY, f64::min),
                    }
                })
                .collect(),
        )
    }

    /// The nadir point: worst observed value per objective on the frontier.
    ///
    /// For Maximize objectives, this is the minimum; for Minimize, the maximum.
    /// Returns `None` if the frontier is empty.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Minimize]);
    /// f.push(vec![0.9, 10.0], "a");
    /// f.push(vec![0.7, 5.0],  "b");
    ///
    /// let nadir = f.nadir_point().unwrap();
    /// assert!((nadir[0] - 0.7).abs() < 1e-9);  // worst accuracy
    /// assert!((nadir[1] - 10.0).abs() < 1e-9);  // worst (highest) latency
    /// ```
    pub fn nadir_point(&self) -> Option<Vec<f64>> {
        if self.is_empty() {
            return None;
        }
        Some(
            (0..self.directions.len())
                .map(|i| {
                    let vals = self.points.iter().map(|p| p.values[i]);
                    match self.directions[i] {
                        Direction::Maximize => vals.fold(f64::INFINITY, f64::min),
                        Direction::Minimize => vals.fold(f64::NEG_INFINITY, f64::max),
                    }
                })
                .collect(),
        )
    }

    /// Normalized objective values for a point, mapped to `[0, 1]`.
    ///
    /// Each dimension is normalized using the frontier's ideal and nadir points:
    /// `0.0` = nadir (worst on the frontier), `1.0` = ideal (best on the frontier).
    /// Direction is accounted for: both Maximize and Minimize objectives produce
    /// values where higher = better after normalization.
    ///
    /// Returns `None` if the frontier is empty. Returns `0.5` for dimensions
    /// where ideal == nadir (no spread).
    ///
    /// This is the normalization that ASF and other MCDM methods assume.
    /// For feeding normalized data to downstream systems (visualization, export),
    /// these values are unit-free and directly comparable across objectives.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Minimize]);
    /// f.push(vec![0.9, 10.0], "a");  // best accuracy, worst latency
    /// f.push(vec![0.7, 5.0],  "b");  // worst accuracy, best latency
    ///
    /// let norm_a = f.normalized_values(0).unwrap();
    /// assert!((norm_a[0] - 1.0).abs() < 1e-9);  // best in dim 0
    /// assert!((norm_a[1] - 0.0).abs() < 1e-9);  // worst in dim 1
    ///
    /// let norm_b = f.normalized_values(1).unwrap();
    /// assert!((norm_b[0] - 0.0).abs() < 1e-9);  // worst in dim 0
    /// assert!((norm_b[1] - 1.0).abs() < 1e-9);  // best in dim 1
    /// ```
    pub fn normalized_values(&self, point_idx: usize) -> Option<Vec<f64>> {
        if self.is_empty() {
            return None;
        }
        assert!(
            point_idx < self.points.len(),
            "point_idx ({point_idx}) out of bounds (len={})",
            self.points.len(),
        );
        let ideal = self.ideal_point()?;
        let nadir = self.nadir_point()?;
        let p = &self.points[point_idx];

        Some(
            (0..self.directions.len())
                .map(|i| {
                    let range = (ideal[i] - nadir[i]).abs();
                    if range < self.eps {
                        return 0.5;
                    }
                    match self.directions[i] {
                        Direction::Maximize => (p.values[i] - nadir[i]) / range,
                        Direction::Minimize => (nadir[i] - p.values[i]) / range,
                    }
                })
                .collect(),
        )
    }

    /// Suggest a reference point for hypervolume computation.
    ///
    /// Returns the nadir point shifted outward by a margin (default: 10% of
    /// the objective range per dimension). This follows the standard practice
    /// from Guerreiro et al. (2021): choose a point "clearly worse than
    /// acceptable" as the reference.
    ///
    /// The margin prevents zero-volume contributions from boundary points.
    /// Use the same reference point across runs for comparable hypervolume values.
    ///
    /// Returns `None` if the frontier is empty.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Minimize]);
    /// f.push(vec![0.9, 10.0], ());
    /// f.push(vec![0.7, 5.0],  ());
    ///
    /// let ref_pt = f.suggest_ref_point(0.1).unwrap();
    /// // For Maximize: nadir - margin * range = 0.7 - 0.1 * 0.2 = 0.68
    /// // For Minimize: nadir + margin * range = 10.0 + 0.1 * 5.0 = 10.5
    /// let hv = f.hypervolume(&ref_pt);
    /// assert!(hv > 0.0);
    /// ```
    pub fn suggest_ref_point(&self, margin: f64) -> Option<Vec<f64>> {
        if self.is_empty() {
            return None;
        }
        let ideal = self.ideal_point()?;
        let nadir = self.nadir_point()?;

        Some(
            (0..self.directions.len())
                .map(|i| {
                    let range = (ideal[i] - nadir[i]).abs();
                    let shift = margin * range;
                    match self.directions[i] {
                        // Nadir is the worst; shift further in the "worse" direction
                        Direction::Maximize => nadir[i] - shift,
                        Direction::Minimize => nadir[i] + shift,
                    }
                })
                .collect(),
        )
    }

    /// Return all frontier point indices sorted by descending weighted score.
    ///
    /// This is the ranked version of [`best_index`](Self::best_index): instead of
    /// returning only the top-scoring point, it returns all points ordered from
    /// best to worst according to the weighted linear score.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Maximize]);
    /// f.push(vec![0.9, 0.1], "A");
    /// f.push(vec![0.5, 0.5], "B");
    /// f.push(vec![0.1, 0.9], "C");
    ///
    /// let ranked = f.ranked_indices(&[1.0, 0.0]); // only care about dim 0
    /// assert_eq!(ranked[0], 0); // A wins
    /// assert_eq!(ranked[2], 2); // C last
    /// ```
    pub fn ranked_indices(&self, weights: &[f64]) -> Vec<usize> {
        let mut indices: Vec<usize> = (0..self.points.len()).collect();
        indices.sort_by(|&a, &b| {
            let sa = self.scalar_score(a, weights);
            let sb = self.scalar_score(b, weights);
            sb.partial_cmp(&sa).unwrap_or(Ordering::Equal)
        });
        indices
    }

    /// Achievement Scalarizing Function (ASF) score for a point.
    ///
    /// ASF is the standard MCDM method for selecting a point on the Pareto front
    /// closest to an ideal direction. Unlike weighted-sum scoring, ASF works on
    /// non-convex fronts.
    ///
    /// $$
    /// \text{ASF}(f, w, z^*) = \max_i \frac{f_i' - z_i^{*\prime}}{w_i}
    /// $$
    ///
    /// where primed values are oriented so that lower is better (Minimize directions
    /// are kept as-is; Maximize directions are negated). The point with the
    /// **lowest** ASF value is the best compromise.
    ///
    /// `ideal` should be the ideal point (use [`ideal_point`](Self::ideal_point)).
    /// Weights control the preferred tradeoff direction; equal weights give equal
    /// importance to all objectives.
    ///
    /// Returns `f64::INFINITY` if any weight is zero (division by zero).
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Minimize]);
    /// f.push(vec![0.9, 10.0], "a");
    /// f.push(vec![0.7, 5.0],  "b");
    ///
    /// let ideal = f.ideal_point().unwrap();
    /// let sa = f.asf(0, &[1.0, 1.0], &ideal);
    /// let sb = f.asf(1, &[1.0, 1.0], &ideal);
    /// // Both are on the front; scores depend on distance from ideal
    /// assert!(sa >= 0.0);
    /// assert!(sb >= 0.0);
    /// ```
    pub fn asf(&self, point_idx: usize, weights: &[f64], ideal: &[f64]) -> f64 {
        assert!(
            point_idx < self.points.len(),
            "point_idx ({point_idx}) out of bounds (len={})",
            self.points.len(),
        );
        let p = &self.points[point_idx];
        let dims = p.values.len().min(weights.len()).min(ideal.len());
        let mut max_val = f64::NEG_INFINITY;
        for i in 0..dims {
            let w = weights[i];
            if w.abs() < 1e-30 {
                return f64::INFINITY;
            }
            // Orient so lower is better
            let (fi, zi) = match self.directions[i] {
                Direction::Maximize => (-p.values[i], -ideal[i]),
                Direction::Minimize => (p.values[i], ideal[i]),
            };
            let val = (fi - zi) / w;
            if val > max_val {
                max_val = val;
            }
        }
        max_val
    }

    /// Find the index of the point with the lowest ASF score (best compromise).
    ///
    /// This is the MCDM selection method: given weights expressing your preference
    /// direction and the ideal point, returns the frontier point closest to that
    /// ideal along the weighted direction.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Maximize]);
    /// f.push(vec![0.9, 0.1], "A");
    /// f.push(vec![0.5, 0.5], "B");
    /// f.push(vec![0.1, 0.9], "C");
    ///
    /// let ideal = f.ideal_point().unwrap();
    /// let best = f.best_asf(&[1.0, 1.0], &ideal).unwrap();
    /// // B is the most balanced point
    /// assert_eq!(f.points()[best].data, "B");
    /// ```
    pub fn best_asf(&self, weights: &[f64], ideal: &[f64]) -> Option<usize> {
        if self.is_empty() {
            return None;
        }
        (0..self.points.len()).min_by(|&a, &b| {
            let sa = self.asf(a, weights, ideal);
            let sb = self.asf(b, weights, ideal);
            sa.partial_cmp(&sb).unwrap_or(Ordering::Equal)
        })
    }

    /// Remove points that do not satisfy a predicate, then recompute stats.
    ///
    /// This is useful for applying post-hoc constraints (e.g., budget limits)
    /// to a frontier without rebuilding it from scratch.
    ///
    /// **Warning:** after `retain`, the remaining points are still mutually
    /// non-dominated (removing points cannot create new dominance relationships),
    /// but they may no longer form the Pareto front of the original candidate set.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Minimize]);
    /// f.push(vec![0.9, 100.0], "expensive");
    /// f.push(vec![0.7, 20.0],  "cheap");
    /// f.push(vec![0.5, 5.0],   "cheapest");
    ///
    /// f.retain(|p| p.values[1] < 50.0); // budget constraint on cost
    /// assert_eq!(f.len(), 2);
    /// ```
    pub fn retain<F>(&mut self, f: F)
    where
        F: FnMut(&Point<V>) -> bool,
    {
        self.points.retain(f);
        self.update_stats();
    }

    /// Calculate the hypervolume of the frontier relative to a reference point.
    ///
    /// **Reference point selection matters.** Hypervolume values are only comparable for
    /// the same `ref_point` and the same objective scaling.  Choose `ref_point` as a
    /// "clearly worse than acceptable" vector (a nadir point), not merely the origin
    /// unless your objectives are literally improvements over zero.  Keep `ref_point`
    /// fixed when comparing progress across runs or algorithms.
    ///
    /// For mixed maximize/minimize objectives, `ref_point` is interpreted in the original
    /// objective units (this method handles the orientation internally).
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Maximize]);
    /// f.push(vec![1.0, 1.0], ());
    ///
    /// // Hypervolume relative to origin: 1.0 * 1.0 = 1.0
    /// let hv = f.hypervolume(&[0.0, 0.0]);
    /// assert!((hv - 1.0).abs() < 1e-6);
    /// ```
    pub fn hypervolume(&self, ref_point: &[f64]) -> f64 {
        let dim = self.directions.len();
        if dim == 0 || self.is_empty() {
            return 0.0;
        }
        assert_eq!(
            ref_point.len(),
            dim,
            "ref_point len ({}) must match directions len ({dim})",
            ref_point.len(),
        );

        // Convert to an oriented, all-maximize coordinate system rooted at `ref_point`.
        // For each dimension, we measure "improvement over reference", clamped at 0.
        //
        // This lets us compute hypervolume against the origin in a consistent way:
        // - Maximize:  oriented = value - ref
        // - Minimize:  oriented = ref - value
        //
        // Hypervolume is then the Lebesgue measure of the union of axis-aligned boxes
        // [0, p] for each point p in oriented space.
        let mut oriented: Vec<Vec<f64>> = self
            .points
            .iter()
            .map(|p| {
                p.values
                    .iter()
                    .enumerate()
                    .map(|(i, &v)| match self.directions[i] {
                        Direction::Maximize => (v - ref_point[i]).max(0.0),
                        Direction::Minimize => (ref_point[i] - v).max(0.0),
                    })
                    .collect::<Vec<_>>()
            })
            .collect();

        // Remove points that contribute no volume (any dim <= 0 means their box has 0 volume).
        oriented.retain(|p| p.iter().all(|&x| x > self.eps));
        if oriented.is_empty() {
            return 0.0;
        }

        // In oriented space, dominance is pure-maximize.
        let oriented = nondominated_max(&oriented, self.eps);
        hypervolume_max_exact(&oriented, dim, self.eps)
    }

    /// Hypervolume contribution of each point on the frontier.
    ///
    /// The contribution of point `i` is `HV(all) - HV(all \ {i})`: how much
    /// hypervolume would be lost if that point were removed. Points with high
    /// contributions are most valuable to the frontier's quality.
    ///
    /// This is the building block for indicator-based selection (SMS-EMOA, IBEA).
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Maximize]);
    /// f.push(vec![1.0, 0.5], ());
    /// f.push(vec![0.5, 1.0], ());
    ///
    /// let contribs = f.hypervolume_contributions(&[0.0, 0.0]);
    /// assert_eq!(contribs.len(), 2);
    /// // Both points contribute positive volume
    /// assert!(contribs.iter().all(|&c| c > 0.0));
    /// // Sum of contributions >= total HV (overlap possible in general,
    /// // but for 2 non-dominated points they partition the space)
    /// ```
    pub fn hypervolume_contributions(&self, ref_point: &[f64]) -> Vec<f64> {
        let n = self.points.len();
        if n == 0 {
            return Vec::new();
        }
        let total = self.hypervolume(ref_point);
        let dim = self.directions.len();

        (0..n)
            .map(|skip| {
                // Build oriented points excluding `skip`
                let mut oriented: Vec<Vec<f64>> = self
                    .points
                    .iter()
                    .enumerate()
                    .filter(|&(j, _)| j != skip)
                    .map(|(_, p)| {
                        p.values
                            .iter()
                            .enumerate()
                            .map(|(i, &v)| match self.directions[i] {
                                Direction::Maximize => (v - ref_point[i]).max(0.0),
                                Direction::Minimize => (ref_point[i] - v).max(0.0),
                            })
                            .collect::<Vec<_>>()
                    })
                    .collect();

                oriented.retain(|p| p.iter().all(|&x| x > self.eps));
                if oriented.is_empty() {
                    return total;
                }
                let oriented = nondominated_max(&oriented, self.eps);
                let hv_without = hypervolume_max_exact(&oriented, dim, self.eps);
                (total - hv_without).max(0.0)
            })
            .collect()
    }

    /// Identify the high-tradeoff (knee) point on the frontier.
    ///
    /// The knee point is where the marginal tradeoff between objectives is
    /// steepest -- small movement along the front causes large changes in
    /// objective values. This is typically the most "balanced" compromise.
    ///
    /// Uses the maximum normalized distance from the hyperplane connecting
    /// the extreme points. For 2D, this is the point farthest from the line
    /// between the endpoints. For higher dimensions, it generalizes to the
    /// point farthest from the convex hull of the extremes.
    ///
    /// Returns `None` if the frontier has fewer than 3 points (no interior
    /// points to evaluate).
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let mut f = ParetoFrontier::new(vec![Direction::Maximize, Direction::Maximize]);
    /// f.push(vec![1.0, 0.0], "extreme_a");
    /// f.push(vec![0.6, 0.6], "knee");
    /// f.push(vec![0.0, 1.0], "extreme_b");
    ///
    /// let knee = f.knee_index().unwrap();
    /// assert_eq!(f.points()[knee].data, "knee");
    /// ```
    pub fn knee_index(&self) -> Option<usize> {
        let n = self.points.len();
        if n < 3 {
            return None;
        }
        let dim = self.directions.len();

        // Normalize all points to [0,1] using ideal/nadir
        let ideal = self.ideal_point()?;
        let nadir = self.nadir_point()?;

        let normalize = |p: &Point<V>| -> Vec<f64> {
            (0..dim)
                .map(|i| {
                    let range = (ideal[i] - nadir[i]).abs();
                    if range < self.eps {
                        return 0.5;
                    }
                    match self.directions[i] {
                        Direction::Maximize => (p.values[i] - nadir[i]) / range,
                        Direction::Minimize => (nadir[i] - p.values[i]) / range,
                    }
                })
                .collect()
        };

        let normed: Vec<Vec<f64>> = self.points.iter().map(normalize).collect();

        // For each point, compute the distance to the line/hyperplane
        // connecting the extreme points. Use the standard formula:
        // distance = |sum_i(p_i) - 1| / sqrt(d)  for the unit simplex normal.
        // This works because normalized extreme points are near the axes.
        let inv_sqrt_d = 1.0 / (dim as f64).sqrt();

        let mut best_idx = 0;
        let mut best_dist = f64::NEG_INFINITY;

        for (idx, p) in normed.iter().enumerate() {
            // Distance from the hyperplane sum(x_i) = 1 (the simplex face)
            let sum: f64 = p.iter().sum();
            let dist = (sum - 1.0).abs() * inv_sqrt_d;
            // Prefer points above the hyperplane (convex knee)
            let signed_dist = if sum >= 1.0 { dist } else { -dist };
            if signed_dist > best_dist {
                best_dist = signed_dist;
                best_idx = idx;
            }
        }

        Some(best_idx)
    }

    /// Construct a frontier from directions and an iterator of `(values, data)` pairs.
    ///
    /// Each item is passed to [`push`](Self::push); dominated points are discarded
    /// as usual. This is the batch version of calling `new` + repeated `push`.
    ///
    /// ```
    /// use pare::{Direction, ParetoFrontier};
    ///
    /// let items = vec![
    ///     (vec![0.9, 0.1], "A"),
    ///     (vec![0.5, 0.5], "B"),
    ///     (vec![0.4, 0.4], "C"), // dominated by B
    /// ];
    ///
    /// let f = ParetoFrontier::from_points(
    ///     vec![Direction::Maximize, Direction::Maximize],
    ///     items,
    /// );
    /// assert_eq!(f.len(), 2);
    /// ```
    pub fn from_points(
        directions: Vec<Direction>,
        items: impl IntoIterator<Item = (Vec<f64>, V)>,
    ) -> Self {
        let mut frontier = Self::new(directions);
        for (values, data) in items {
            frontier.push(values, data);
        }
        frontier
    }
}

impl<V> Extend<(Vec<f64>, V)> for ParetoFrontier<V> {
    /// Extend the frontier with an iterator of `(values, data)` pairs.
    ///
    /// Each item is passed to [`push`](Self::push); dominated points are
    /// discarded as usual.
    fn extend<I: IntoIterator<Item = (Vec<f64>, V)>>(&mut self, iter: I) {
        for (values, data) in iter {
            self.push(values, data);
        }
    }
}

// ============================================================================
// Epsilon-dominance archiving (Laumanns et al. 2002)
// ============================================================================

/// Grid-based epsilon-dominance archive with bounded size.
///
/// Unlike [`ParetoFrontier::with_eps`] (numerical tolerance for float comparison),
/// this type deliberately coarsens objective space into grid cells of width `eps`
/// per dimension. At most one point is kept per grid cell, guaranteeing an
/// archive size bounded by `prod_i(range_i / eps_i)`.
///
/// Use this for streaming / online optimization where memory must be bounded
/// while maintaining an approximate Pareto front.
///
/// ```
/// use pare::{Direction, EpsilonArchive};
///
/// let mut archive = EpsilonArchive::new_uniform(
///     vec![Direction::Maximize, Direction::Maximize],
///     0.1,
/// );
///
/// archive.push(vec![0.91, 0.12], "a");
/// archive.push(vec![0.92, 0.11], "b"); // same grid cell as "a", replaces if better
/// archive.push(vec![0.50, 0.50], "c"); // different cell, kept
/// assert!(archive.len() <= 2);
/// ```
#[derive(Debug, Clone)]
pub struct EpsilonArchive<V> {
    grid_eps: Vec<f64>,
    directions: Vec<Direction>,
    archive: Vec<Point<V>>,
}

impl<V> EpsilonArchive<V> {
    /// Create an archive with per-dimension grid epsilon values.
    ///
    /// # Panics
    ///
    /// Panics if `grid_eps.len() != directions.len()` or any epsilon is not positive.
    pub fn new(directions: Vec<Direction>, grid_eps: Vec<f64>) -> Self {
        assert_eq!(
            directions.len(),
            grid_eps.len(),
            "directions len ({}) must match grid_eps len ({})",
            directions.len(),
            grid_eps.len(),
        );
        assert!(
            grid_eps.iter().all(|&e| e > 0.0 && e.is_finite()),
            "all grid_eps values must be positive and finite"
        );
        Self {
            grid_eps,
            directions,
            archive: Vec::new(),
        }
    }

    /// Create an archive with a uniform grid epsilon for all dimensions.
    ///
    /// # Panics
    ///
    /// Panics if `eps` is not positive.
    pub fn new_uniform(directions: Vec<Direction>, eps: f64) -> Self {
        let d = directions.len();
        Self::new(directions, vec![eps; d])
    }

    /// Map a value to its grid cell index for a given dimension.
    fn cell(&self, value: f64, dim: usize) -> i64 {
        let oriented = match self.directions[dim] {
            Direction::Maximize => value,
            Direction::Minimize => -value,
        };
        (oriented / self.grid_eps[dim]).floor() as i64
    }

    /// Map a point to its grid cell vector.
    fn cell_vec(&self, values: &[f64]) -> Vec<i64> {
        (0..self.directions.len())
            .map(|i| self.cell(values[i], i))
            .collect()
    }

    /// Returns true if cell `a` epsilon-dominates cell `b`:
    /// `a[i] >= b[i]` for all i, and `a[i] > b[i]` for at least one i.
    fn cell_dominates(a: &[i64], b: &[i64]) -> bool {
        let mut strictly_better = false;
        for (&ai, &bi) in a.iter().zip(b.iter()) {
            if ai < bi {
                return false;
            }
            if ai > bi {
                strictly_better = true;
            }
        }
        strictly_better
    }

    /// Insert a point into the archive if it is not epsilon-dominated.
    ///
    /// Returns `true` if the point was added (or replaced an existing point
    /// in the same cell), `false` if it was rejected.
    ///
    /// # Panics
    ///
    /// Panics if `values.len() != directions.len()` or any value is non-finite.
    pub fn push(&mut self, values: Vec<f64>, data: V) -> bool {
        assert_eq!(
            values.len(),
            self.directions.len(),
            "values len ({}) must match directions len ({})",
            values.len(),
            self.directions.len(),
        );
        assert!(
            values.iter().all(|v| v.is_finite()),
            "all values must be finite"
        );

        let new_cell = self.cell_vec(&values);

        // Check if any existing point's cell epsilon-dominates the new cell
        for existing in &self.archive {
            let ex_cell = self.cell_vec(&existing.values);
            if Self::cell_dominates(&ex_cell, &new_cell) {
                return false;
            }
            // Same cell: keep the one that's better in the raw values
            if ex_cell == new_cell {
                // Compare raw values: count dimensions where new is strictly better
                let new_better = self.raw_dominates(&values, &existing.values);
                if !new_better {
                    return false; // existing is at least as good
                }
                // New point is better in raw values; will replace below
            }
        }

        // Pre-compute cell vectors to avoid borrowing self in the closure.
        let existing_cells: Vec<Vec<i64>> = self
            .archive
            .iter()
            .map(|p| self.cell_vec(&p.values))
            .collect();

        // Remove points whose cells are epsilon-dominated by the new cell,
        // or that share the same cell (being replaced).
        let mut idx = 0;
        self.archive.retain(|_| {
            let keep = !Self::cell_dominates(&new_cell, &existing_cells[idx])
                && existing_cells[idx] != new_cell;
            idx += 1;
            keep
        });

        self.archive.push(Point { values, data });
        true
    }

    /// Raw dominance check (oriented by directions, no epsilon).
    fn raw_dominates(&self, a: &[f64], b: &[f64]) -> bool {
        let mut strictly_better = false;
        for (i, (&av, &bv)) in a.iter().zip(b.iter()).enumerate() {
            match self.directions[i] {
                Direction::Maximize => {
                    if av < bv {
                        return false;
                    }
                    if av > bv {
                        strictly_better = true;
                    }
                }
                Direction::Minimize => {
                    if av > bv {
                        return false;
                    }
                    if av < bv {
                        strictly_better = true;
                    }
                }
            }
        }
        strictly_better
    }

    /// Number of points in the archive.
    pub fn len(&self) -> usize {
        self.archive.len()
    }

    /// Is the archive empty?
    pub fn is_empty(&self) -> bool {
        self.archive.is_empty()
    }

    /// Get all points in the archive.
    pub fn points(&self) -> &[Point<V>] {
        &self.archive
    }

    /// Get the grid epsilon values.
    pub fn grid_eps(&self) -> &[f64] {
        &self.grid_eps
    }

    /// Get the optimization directions.
    pub fn directions(&self) -> &[Direction] {
        &self.directions
    }

    /// Convert into a standard Pareto frontier for downstream analysis
    /// (crowding distance, hypervolume, scoring, etc.).
    ///
    /// The resulting frontier uses the default numerical epsilon (`1e-9`),
    /// not the grid epsilon. All archived points are inserted; since they
    /// are mutually non-epsilon-dominated, most will survive standard
    /// dominance filtering (though some edge cases near cell boundaries
    /// may be filtered).
    pub fn into_frontier(self) -> ParetoFrontier<V> {
        let mut frontier = ParetoFrontier::new(self.directions);
        for point in self.archive {
            frontier.push(point.values, point.data);
        }
        frontier
    }
}

impl ParetoFrontier<usize> {
    /// Construct a Pareto frontier from raw objective vectors.
    ///
    /// This is a convenience constructor used by the README quickstart:
    /// - the attached `data` is the original index of each point
    /// - all objectives are assumed to be [`Direction::Maximize`]
    ///
    /// For mixed maximize/minimize objectives, build a frontier with
    /// [`ParetoFrontier::new`] and insert points via [`ParetoFrontier::push`].
    ///
    /// ```rust
    /// use pare::ParetoFrontier;
    ///
    /// // [Relevance, Recency] - higher is better
    /// let candidates = vec![
    ///     vec![0.9, 0.1], // A
    ///     vec![0.5, 0.5], // B
    ///     vec![0.1, 0.9], // C
    ///     vec![0.4, 0.4], // D (dominated by B)
    /// ];
    ///
    /// let frontier = ParetoFrontier::try_new(&candidates).unwrap().indices();
    /// assert_eq!(frontier, vec![0, 1, 2]);
    /// ```
    pub fn try_new(points: &[Vec<f64>]) -> Result<Self, FrontierError> {
        if points.is_empty() {
            return Err(FrontierError::Empty);
        }
        let d = points[0].len();
        if d == 0 || points.iter().any(|p| p.len() != d) {
            return Err(FrontierError::InconsistentDimensions);
        }
        for (pi, p) in points.iter().enumerate() {
            for (di, &v) in p.iter().enumerate() {
                if !v.is_finite() {
                    return Err(FrontierError::NonFinite {
                        point_idx: pi,
                        dim_idx: di,
                    });
                }
            }
        }

        let mut frontier = ParetoFrontier::new(vec![Direction::Maximize; d]);
        for (i, p) in points.iter().enumerate() {
            frontier.push(p.clone(), i);
        }
        Ok(frontier)
    }

    /// Return the original indices of points on the frontier.
    ///
    /// This is only meaningful for frontiers constructed with `data = usize`,
    /// such as those produced by [`ParetoFrontier::try_new`].
    ///
    /// **Order note.** The returned indices are in the frontier's internal insertion
    /// order, not sorted by any objective.  If you need a stable ordering, sort the
    /// result explicitly (e.g. by one objective value or by `data`).
    pub fn indices(&self) -> Vec<usize> {
        self.points.iter().map(|p| p.data).collect()
    }
}

/// Return a non-dominated subset in maximize space (all objectives are "higher is better").
fn nondominated_max(points: &[Vec<f64>], eps: f64) -> Vec<Vec<f64>> {
    let mut out: Vec<Vec<f64>> = Vec::new();
    'outer: for p in points {
        // If an existing point dominates p, drop p.
        for q in &out {
            if dominates_max(q, p, eps) {
                continue 'outer;
            }
        }
        // Otherwise, remove points dominated by p and insert p.
        out.retain(|q| !dominates_max(p, q, eps));
        out.push(p.clone());
    }
    out
}

fn dominates_max(a: &[f64], b: &[f64], eps: f64) -> bool {
    let mut strictly_better = false;
    for (&av, &bv) in a.iter().zip(b.iter()) {
        if av + eps < bv {
            return false;
        }
        if av > bv + eps {
            strictly_better = true;
        }
    }
    strictly_better
}

/// Exact hypervolume in maximize space, reference at the origin.
///
/// This uses recursive slicing on the last dimension:
/// $$
/// HV_d(P) = \sum_{k} (y_k - y_{k+1}) \cdot HV_{d-1}(\pi(P_{y \ge y_k}))
/// $$
/// where $y_k$ are the unique last-coordinate levels (descending), and $\pi$ drops the last coordinate.
fn hypervolume_max_exact(points: &[Vec<f64>], dim: usize, eps: f64) -> f64 {
    debug_assert!(dim >= 1);
    if points.is_empty() {
        return 0.0;
    }
    if dim == 1 {
        return points.iter().map(|p| p[0]).fold(0.0, f64::max);
    }
    if dim == 2 {
        return hypervolume_max_2d(points, eps);
    }
    if dim == 3 {
        return hypervolume_max_3d(points, eps);
    }

    // Collect unique slice levels (positive only).
    let mut levels: Vec<f64> = points
        .iter()
        .map(|p| p[dim - 1])
        .filter(|&v| v > eps)
        .collect();
    levels.sort_by(|a, b| b.partial_cmp(a).unwrap_or(Ordering::Equal)); // descending
    levels.dedup_by(|a, b| (*a - *b).abs() <= eps);

    let mut hv = 0.0;
    for (idx, &level) in levels.iter().enumerate() {
        let next = if idx + 1 < levels.len() {
            levels[idx + 1]
        } else {
            0.0
        };
        let thickness = (level - next).max(0.0);
        if thickness <= eps {
            continue;
        }

        let mut projected: Vec<Vec<f64>> = points
            .iter()
            .filter(|p| p[dim - 1] + eps >= level)
            .map(|p| p[0..dim - 1].to_vec())
            .collect();
        if projected.is_empty() {
            continue;
        }
        projected = nondominated_max(&projected, eps);
        let cross = hypervolume_max_exact(&projected, dim - 1, eps);
        hv += thickness * cross;
    }
    hv
}

fn hypervolume_max_2d(points: &[Vec<f64>], eps: f64) -> f64 {
    // Sort by x ascending. For a non-dominated set, y should be non-increasing with x.
    let mut idxs: Vec<usize> = (0..points.len()).collect();
    idxs.sort_by(|&i, &j| {
        let xi = points[i][0];
        let xj = points[j][0];
        match xi.partial_cmp(&xj).unwrap_or(Ordering::Equal) {
            Ordering::Equal => points[j][1]
                .partial_cmp(&points[i][1])
                .unwrap_or(Ordering::Equal),
            ord => ord,
        }
    });

    let mut area = 0.0;
    let mut prev_x = 0.0;
    for i in idxs {
        let x = points[i][0].max(0.0);
        let y = points[i][1].max(0.0);
        let dx = (x - prev_x).max(0.0);
        if dx > eps && y > eps {
            area += dx * y;
        }
        prev_x = x;
    }
    area
}

/// Specialized 3D hypervolume using z-sweep with incremental 2D front.
///
/// O(n^2) worst case (dominated by 2D front maintenance), but avoids the
/// overhead of recursive slicing (projection, re-filtering, re-sorting at
/// each recursion level). For typical non-dominated sets this is significantly
/// faster than the generic recursive path.
fn hypervolume_max_3d(points: &[Vec<f64>], eps: f64) -> f64 {
    if points.is_empty() {
        return 0.0;
    }

    // Sort by z descending.
    let mut sorted: Vec<usize> = (0..points.len()).collect();
    sorted.sort_by(|&i, &j| {
        points[j][2]
            .partial_cmp(&points[i][2])
            .unwrap_or(Ordering::Equal)
    });

    // 2D front maintained as (x, y) pairs sorted by x ascending.
    // Invariant: y is non-increasing with x (since non-dominated).
    let mut front_2d: Vec<(f64, f64)> = Vec::new();
    let mut hv = 0.0;
    let mut prev_z = points[sorted[0]][2];

    // Process first point.
    front_2d.push((points[sorted[0]][0], points[sorted[0]][1]));

    for &idx in sorted.iter().skip(1) {
        let z = points[idx][2];
        let thickness = (prev_z - z).max(0.0);

        if thickness > eps {
            // Compute 2D HV from the current front.
            let hv_2d = sweep_2d_front(&front_2d, eps);
            hv += thickness * hv_2d;
        }

        prev_z = z;

        // Insert (x, y) into the 2D front.
        let x = points[idx][0];
        let y = points[idx][1];
        insert_2d_front(&mut front_2d, x, y, eps);
    }

    // Final slab from last z to 0.
    if prev_z > eps {
        let hv_2d = sweep_2d_front(&front_2d, eps);
        hv += prev_z * hv_2d;
    }

    hv
}

/// Insert (x, y) into a sorted 2D non-dominated front (x ascending, y non-increasing).
/// Removes dominated points after insertion.
fn insert_2d_front(front: &mut Vec<(f64, f64)>, x: f64, y: f64, eps: f64) {
    // Check if dominated: any existing point with x >= new_x and y >= new_y
    for &(fx, fy) in front.iter() {
        if fx + eps >= x && fy + eps >= y {
            return; // dominated
        }
    }

    // Remove points dominated by (x, y)
    front.retain(|&(fx, fy)| !(x + eps >= fx && y + eps >= fy));

    // Insert in sorted position by x
    let pos = front.partition_point(|&(fx, _)| fx < x);
    front.insert(pos, (x, y));
}

/// Compute 2D HV from a sorted non-dominated front (x ascending, y non-increasing).
/// Reference at origin.
fn sweep_2d_front(front: &[(f64, f64)], eps: f64) -> f64 {
    let mut area = 0.0;
    let mut prev_x = 0.0;
    for &(x, y) in front {
        let x = x.max(0.0);
        let y = y.max(0.0);
        let dx = (x - prev_x).max(0.0);
        if dx > eps && y > eps {
            area += dx * y;
        }
        prev_x = x;
    }
    area
}

/// Standalone dominance check.
///
/// Returns `true` if point `a` dominates point `b` (at least as good in all
/// dimensions and strictly better in at least one).
///
/// `a`, `b`, and `directions` must all have the same length; mismatched lengths
/// cause a panic in debug builds and silent truncation in release builds.
///
/// ```
/// use pare::{dominates, Direction};
///
/// let dirs = vec![Direction::Maximize, Direction::Maximize];
/// assert!(dominates(&dirs, 1e-9, &[0.9, 0.8], &[0.5, 0.5]));
/// assert!(!dominates(&dirs, 1e-9, &[0.9, 0.3], &[0.5, 0.5])); // trade-off
/// ```
pub fn dominates(directions: &[Direction], eps: f64, a: &[f64], b: &[f64]) -> bool {
    debug_assert_eq!(
        a.len(),
        b.len(),
        "dominates: a.len() ({}) != b.len() ({})",
        a.len(),
        b.len()
    );
    debug_assert_eq!(
        a.len(),
        directions.len(),
        "dominates: a.len() ({}) != directions.len() ({})",
        a.len(),
        directions.len()
    );
    let mut strictly_better = false;
    for (i, (&av, &bv)) in a.iter().zip(b.iter()).enumerate() {
        let dir = directions[i];
        match dir {
            Direction::Maximize => {
                if av + eps < bv {
                    return false;
                }
                if av > bv + eps {
                    strictly_better = true;
                }
            }
            Direction::Minimize => {
                if av > bv + eps {
                    return false;
                }
                if av + eps < bv {
                    strictly_better = true;
                }
            }
        }
    }
    strictly_better
}

// ============================================================================
// Back-compat helpers (bench + older call sites)
// ============================================================================

/// Return indices of the Pareto frontier (all objectives maximized).
///
/// This is a convenience wrapper intended for benchmarks and legacy code.
/// The canonical API is [`ParetoFrontier`].
///
/// ```
/// use pare::pareto_indices;
///
/// let points = vec![
///     vec![0.9f32, 0.1],
///     vec![0.5, 0.5],
///     vec![0.1, 0.9],
///     vec![0.4, 0.4], // dominated by [0.5, 0.5]
/// ];
///
/// let front = pareto_indices(&points).unwrap();
/// assert!(front.contains(&0));
/// assert!(front.contains(&1));
/// assert!(front.contains(&2));
/// assert!(!front.contains(&3)); // dominated
/// ```
pub fn pareto_indices(points: &[Vec<f32>]) -> Option<Vec<usize>> {
    if points.is_empty() {
        return Some(Vec::new());
    }
    let d = points[0].len();
    if d == 0 || points.iter().any(|p| p.len() != d) {
        return None;
    }
    // Check for non-finite values.
    if points.iter().any(|p| p.iter().any(|v| !v.is_finite())) {
        return None;
    }

    let directions = vec![Direction::Maximize; d];
    let eps = 1e-9;

    let as_f64: Vec<Vec<f64>> = points
        .iter()
        .map(|p| p.iter().map(|&x| x as f64).collect())
        .collect();

    let mut keep = vec![true; points.len()];
    for i in 0..as_f64.len() {
        if !keep[i] {
            continue;
        }
        for j in 0..as_f64.len() {
            if i == j || !keep[j] {
                continue;
            }
            if dominates(&directions, eps, &as_f64[j], &as_f64[i]) {
                keep[i] = false;
                break; // i is dominated; no need to check further
            }
        }
    }

    Some(
        keep.into_iter()
            .enumerate()
            .filter_map(|(i, ok)| ok.then_some(i))
            .collect(),
    )
}

/// Return indices of the Pareto frontier using a simple 2D sweep (maximize/maximize).
///
/// Falls back to [`pareto_indices`] if `points` are not 2D.
pub fn pareto_indices_2d(points: &[Vec<f32>]) -> Option<Vec<usize>> {
    if points.is_empty() {
        return Some(Vec::new());
    }
    if points.iter().any(|p| p.len() != 2) {
        return pareto_indices(points);
    }
    // Check for non-finite values.
    if points.iter().any(|p| p.iter().any(|v| !v.is_finite())) {
        return None;
    }

    // Sort by x descending, then y descending.
    let mut idxs: Vec<usize> = (0..points.len()).collect();
    idxs.sort_by(|&i, &j| {
        let xi = points[i][0];
        let xj = points[j][0];
        xj.partial_cmp(&xi)
            .unwrap_or(Ordering::Equal)
            .then_with(|| {
                points[j][1]
                    .partial_cmp(&points[i][1])
                    .unwrap_or(Ordering::Equal)
            })
    });

    let mut best_y = f32::NEG_INFINITY;
    let mut out = Vec::new();
    for i in idxs {
        let y = points[i][1];
        if y > best_y {
            out.push(i);
            best_y = y;
        }
    }
    Some(out)
}

/// Partition points into successive non-dominated layers (Pareto layers).
///
/// Layer 0 is the Pareto front. Layer 1 is the front of the remaining points
/// after removing layer 0, and so on. All objectives are maximized.
///
/// Returns `None` if any point has inconsistent dimensions or non-finite values.
///
/// ```
/// use pare::pareto_layers;
///
/// let points = vec![
///     vec![0.9f32, 0.1], // layer 0
///     vec![0.5, 0.5],    // layer 0
///     vec![0.4, 0.4],    // layer 1 (dominated by [0.5, 0.5])
///     vec![0.3, 0.3],    // layer 2
/// ];
///
/// let layers = pareto_layers(&points).unwrap();
/// assert!(layers[0].contains(&0));
/// assert!(layers[0].contains(&1));
/// assert_eq!(layers[1], vec![2]);
/// assert_eq!(layers[2], vec![3]);
/// ```
pub fn pareto_layers(points: &[Vec<f32>]) -> Option<Vec<Vec<usize>>> {
    if points.is_empty() {
        return Some(Vec::new());
    }
    let d = points[0].len();
    if d == 0 || points.iter().any(|p| p.len() != d) {
        return None;
    }
    if points.iter().any(|p| p.iter().any(|v| !v.is_finite())) {
        return None;
    }

    let directions = vec![Direction::Maximize; d];
    let eps = 1e-9_f64;
    let as_f64: Vec<Vec<f64>> = points
        .iter()
        .map(|p| p.iter().map(|&x| x as f64).collect())
        .collect();

    let mut remaining: Vec<usize> = (0..points.len()).collect();
    let mut layers = Vec::new();

    while !remaining.is_empty() {
        // Find the non-dominated set among remaining
        let mut keep = vec![true; remaining.len()];
        for i in 0..remaining.len() {
            if !keep[i] {
                continue;
            }
            for j in 0..remaining.len() {
                if i == j || !keep[j] {
                    continue;
                }
                if dominates(
                    &directions,
                    eps,
                    &as_f64[remaining[j]],
                    &as_f64[remaining[i]],
                ) {
                    keep[i] = false;
                    break;
                }
            }
        }

        let layer: Vec<usize> = remaining
            .iter()
            .zip(keep.iter())
            .filter_map(|(&idx, &k)| k.then_some(idx))
            .collect();

        remaining = remaining
            .into_iter()
            .zip(keep.iter())
            .filter_map(|(idx, &k)| (!k).then_some(idx))
            .collect();

        layers.push(layer);
    }

    Some(layers)
}

/// Return indices of the k-dominant frontier (maximize/maximize/...).
///
/// Definition used here: point `a` k-dominates `b` if `a` is strictly better
/// in at least `k` dimensions and is not worse in any dimension.
pub fn pareto_indices_k_dominance(points: &[Vec<f32>], k: usize) -> Option<Vec<usize>> {
    if points.is_empty() {
        return Some(Vec::new());
    }
    let d = points[0].len();
    if d == 0 || points.iter().any(|p| p.len() != d) {
        return None;
    }
    // Check for non-finite values.
    if points.iter().any(|p| p.iter().any(|v| !v.is_finite())) {
        return None;
    }
    let k = k.min(d);
    let eps = 1e-9_f64;

    let as_f64: Vec<Vec<f64>> = points
        .iter()
        .map(|p| p.iter().map(|&x| x as f64).collect())
        .collect();

    let mut keep = vec![true; points.len()];
    for i in 0..as_f64.len() {
        if !keep[i] {
            continue;
        }
        'cmp: for j in 0..as_f64.len() {
            if i == j || !keep[i] {
                continue;
            }
            let mut better = 0usize;
            for (&aj, &ai) in as_f64[j].iter().zip(as_f64[i].iter()).take(d) {
                if aj + eps < ai {
                    continue 'cmp; // j is worse in this dim -> cannot dominate i
                }
                if aj > ai + eps {
                    better += 1;
                }
            }
            if better >= k {
                keep[i] = false;
            }
        }
    }

    Some(
        keep.into_iter()
            .enumerate()
            .filter_map(|(i, ok)| ok.then_some(i))
            .collect(),
    )
}

// ============================================================================
// Quality indicators
// ============================================================================

/// Euclidean distance between two points.
fn euclidean_dist(a: &[f64], b: &[f64]) -> f64 {
    a.iter()
        .zip(b.iter())
        .map(|(&ai, &bi)| (ai - bi) * (ai - bi))
        .sum::<f64>()
        .sqrt()
}

/// Generational Distance (GD) from a front to a reference set.
///
/// Measures how far each point in `front` is from the nearest point in
/// `reference`. Lower is better; zero means every front point lies on the
/// reference set.
///
/// $$
/// \text{GD}(F, R) = \frac{1}{|F|} \sum_{f \in F} \min_{r \in R} \| f - r \|
/// $$
///
/// Returns `None` if either set is empty or dimensions are inconsistent.
///
/// ```
/// use pare::generational_distance;
///
/// let front = vec![vec![0.5, 0.5]];
/// let reference = vec![vec![0.5, 0.5], vec![1.0, 0.0]];
/// let gd = generational_distance(&front, &reference).unwrap();
/// assert!(gd < 1e-9); // front point is on the reference set
/// ```
pub fn generational_distance(front: &[Vec<f64>], reference: &[Vec<f64>]) -> Option<f64> {
    if front.is_empty() || reference.is_empty() {
        return None;
    }
    let d = front[0].len();
    if d == 0 || front.iter().any(|p| p.len() != d) || reference.iter().any(|p| p.len() != d) {
        return None;
    }

    let sum: f64 = front
        .iter()
        .map(|f| {
            reference
                .iter()
                .map(|r| euclidean_dist(f, r))
                .fold(f64::INFINITY, f64::min)
        })
        .sum();

    Some(sum / front.len() as f64)
}

/// Inverted Generational Distance (IGD) from a front to a reference set.
///
/// Measures how well the front covers the reference set. For each reference
/// point, finds the nearest front point. Lower is better; zero means every
/// reference point is covered by the front.
///
/// $$
/// \text{IGD}(F, R) = \frac{1}{|R|} \sum_{r \in R} \min_{f \in F} \| r - f \|
/// $$
///
/// IGD is generally preferred over GD because it captures both convergence
/// and spread. A front that covers only one region of the reference set
/// will have low GD but high IGD.
///
/// Returns `None` if either set is empty or dimensions are inconsistent.
///
/// ```
/// use pare::inverted_generational_distance;
///
/// let front = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
/// let reference = vec![vec![1.0, 0.0], vec![0.5, 0.5], vec![0.0, 1.0]];
/// let igd = inverted_generational_distance(&front, &reference).unwrap();
/// // front covers the extremes but not the middle
/// assert!(igd > 0.0);
/// ```
pub fn inverted_generational_distance(front: &[Vec<f64>], reference: &[Vec<f64>]) -> Option<f64> {
    // IGD(F, R) = GD(R, F) -- swap the roles
    generational_distance(reference, front)
}