metalforge 0.3.0

//! PA-DDS — Pareto-Archived Dynamically Dimensioned Search
//! (Asadzadeh & Tolson 2013, Engineering Optimization 45(12), 1489–1509).
//!
//! The multi-objective extension of [`Dds`](super::Dds), built for the same
//! niche: **budget-limited water-resources calibration** where each
//! evaluation is a model run. It keeps DDS's defining property — a single
//! algorithm parameter (`r`) and a perturbation schedule that narrows with
//! the remaining budget — and adds an unbounded **external archive** of
//! non-dominated solutions:
//!
//! 1. Each iteration selects a *current* solution from the archive by a
//!    selection metric; the paper's best-performing variant (and the default
//!    here) weights members by their exclusive **hypervolume contribution**,
//!    so under-covered regions of the front get refined preferentially.
//! 2. The current solution is perturbed exactly like DDS: each variable with
//!    probability `P(i) = 1 − ln(i)/ln(m)` (at least one), reflecting at the
//!    bounds.
//! 3. A candidate that is not dominated by the archive is inserted and the
//!    members it dominates are removed.
//!
//! Objectives are normalized by the archive's per-objective ranges before
//! computing contributions (reference point 1.1 per axis), so differently
//! scaled objectives (e.g. NSE vs log-flow bias) do not distort selection.
//! As with [`SmsEmoa`](super::SmsEmoa), hypervolume selection is intended
//! for **2-3 objectives**; use [`PaddsSelection::Crowding`] (or `Random`)
//! beyond that. Deterministic for a given seed.

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

/// How PA-DDS selects the current solution from the archive (the paper
/// studies exactly these three).
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[non_exhaustive]
pub enum PaddsSelection {
    /// Uniformly random archive member.
    Random,
    /// Probability proportional to NSGA-II crowding distance (boundary
    /// members, whose crowding is infinite, are weighted as twice the largest
    /// finite crowding).
    Crowding,
    /// Probability proportional to the exclusive hypervolume contribution on
    /// range-normalized objectives (the paper's recommended metric).
    Hypervolume,
}

/// PA-DDS configuration.
#[derive(Debug, Clone, Copy)]
pub struct Padds {
    /// DDS neighborhood size as a fraction of each variable's range
    /// (the paper inherits DDS's robust default `0.2`).
    pub r: f64,
    /// Archive selection metric.
    pub selection: PaddsSelection,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for Padds {
    fn default() -> Self {
        Padds {
            r: 0.2,
            selection: PaddsSelection::Hypervolume,
            seed: 42,
        }
    }
}

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

impl Padds {
    /// Approximates the Pareto front of `problem` within the evaluation
    /// budget of `term` (its `target` is ignored). Returns the full external
    /// archive.
    pub fn optimize(&self, problem: &dyn MultiProblem, term: &Termination) -> ParetoFront {
        crate::problem::validate_multi(problem)
            .unwrap_or_else(|e| panic!("Padds: invalid problem: {e}"));
        let bounds = problem.bounds();
        let dim = bounds.len();
        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;
                }
            }
            o
        };

        let mut evaluations = 0usize;
        if term.max_evaluations == 0 {
            return ParetoFront {
                solutions: Vec::new(),
                evaluations,
            };
        }

        // Initial solution seeds the archive.
        let x0: Vec<f64> = bounds
            .iter()
            .map(|&(lo, hi)| rng.uniform_in(lo, hi))
            .collect();
        let o0 = eval(&x0);
        evaluations += 1;
        let mut archive: Vec<MultiSolution> = vec![MultiSolution {
            x: x0,
            objectives: o0,
        }];

        // DDS schedule over the candidate evaluations after the initial one.
        let m = (term.max_evaluations.saturating_sub(1)).max(2) as f64;
        let mut i = 1usize;
        let mut candidate = vec![0.0; dim];
        while evaluations < term.max_evaluations {
            let cur = self.select(&archive, &mut rng);
            let base = &archive[cur].x;

            let p = 1.0 - (i as f64).ln() / m.ln();
            candidate.copy_from_slice(base);
            let mut perturbed = 0;
            for (j, &(lo, hi)) in bounds.iter().enumerate() {
                if rng.uniform() < p {
                    candidate[j] = perturb(base[j], lo, hi, self.r, &mut rng);
                    perturbed += 1;
                }
            }
            if perturbed == 0 {
                let j = rng.index(dim);
                let (lo, hi) = bounds[j];
                candidate[j] = perturb(base[j], lo, hi, self.r, &mut rng);
            }

            let objs = eval(&candidate);
            evaluations += 1;
            i += 1;

            // Archive update: insert if no member dominates the candidate;
            // drop the members the candidate dominates. Exact duplicates of
            // an archive member are discarded (they add nothing).
            let dominated_or_dup = archive
                .iter()
                .any(|s| dominates(&s.objectives, &objs) || s.objectives == objs);
            if !dominated_or_dup {
                archive.retain(|s| !dominates(&objs, &s.objectives));
                archive.push(MultiSolution {
                    x: candidate.clone(),
                    objectives: objs,
                });
            }
        }

        ParetoFront {
            solutions: archive,
            evaluations,
        }
    }

    /// Selects the current solution's archive index by the configured metric.
    fn select(&self, archive: &[MultiSolution], rng: &mut Rng) -> usize {
        let n = archive.len();
        if n == 1 {
            return 0;
        }
        let weights: Vec<f64> = match self.selection {
            PaddsSelection::Random => return rng.index(n),
            PaddsSelection::Crowding => {
                let all: Vec<usize> = (0..n).collect();
                let mut c = crowding_distance(&all, archive);
                // Boundary members have infinite crowding: weight them as 2×
                // the largest finite crowding so they stay favored but finite.
                let max_finite = c
                    .iter()
                    .copied()
                    .filter(|v| v.is_finite())
                    .fold(0.0, f64::max);
                let cap = if max_finite > 0.0 {
                    2.0 * max_finite
                } else {
                    1.0
                };
                for v in &mut c {
                    if !v.is_finite() {
                        *v = cap;
                    }
                }
                c
            }
            PaddsSelection::Hypervolume => {
                // Range-normalize the archive so scale differences between
                // objectives do not distort the contributions.
                let mobj = archive[0].objectives.len();
                let mut lo = vec![f64::INFINITY; mobj];
                let mut hi = vec![f64::NEG_INFINITY; mobj];
                for s in archive {
                    for (j, &v) in s.objectives.iter().enumerate() {
                        lo[j] = lo[j].min(v);
                        hi[j] = hi[j].max(v);
                    }
                }
                let normalized: Vec<Vec<f64>> = archive
                    .iter()
                    .map(|s| {
                        (0..mobj)
                            .map(|j| {
                                let span = hi[j] - lo[j];
                                if span > 0.0 && span.is_finite() {
                                    (s.objectives[j] - lo[j]) / span
                                } else {
                                    0.0
                                }
                            })
                            .collect()
                    })
                    .collect();
                let r = vec![1.1; mobj];
                crate::indicators::hv_contributions(&normalized, &r)
            }
        };

        // Roulette over the weights; degenerate weights fall back to uniform.
        let total: f64 = weights.iter().sum();
        if !(total.is_finite() && total > 0.0) {
            return rng.index(n);
        }
        let u = rng.uniform() * total;
        let mut acc = 0.0;
        for (idx, w) in weights.iter().enumerate() {
            acc += w;
            if u < acc {
                return idx;
            }
        }
        n - 1
    }
}

/// One-dimensional DDS neighborhood move with boundary reflection
/// (identical to [`Dds`](super::Dds)'s rule, Tolson & Shoemaker 2007 eq. 4).
fn perturb(x: f64, lo: f64, hi: f64, r: f64, rng: &mut Rng) -> f64 {
    let range = hi - lo;
    let mut xn = x + range * r * rng.normal();
    if xn < lo {
        xn = lo + (lo - xn);
        if xn > hi {
            xn = lo;
        }
    } else if xn > hi {
        xn = hi - (xn - hi);
        if xn < lo {
            xn = hi;
        }
    }
    xn.clamp(lo, hi)
}