metalforge 0.1.0

//! Standard optimization test functions, as [`Problem`]s.
//!
//! These are the classic benchmarks used to validate global optimizers: the
//! convex unimodal **sphere**, the narrow-valley **Rosenbrock**, and the highly
//! multimodal **Rastrigin**. All have a known global minimum of `0`, which the
//! integration tests check each algorithm can reach.

use crate::problem::{Bound, MultiProblem, Problem};

/// Sphere: `f(x) = Σ xᵢ²`. Convex, unimodal, minimum `0` at the origin.
pub struct Sphere {
    bounds: Vec<Bound>,
}

impl Sphere {
    /// `dim`-dimensional sphere on `[-5.12, 5.12]^dim`.
    pub fn new(dim: usize) -> Self {
        Sphere {
            bounds: vec![(-5.12, 5.12); dim],
        }
    }
}

impl Problem for Sphere {
    fn dim(&self) -> usize {
        self.bounds.len()
    }
    fn bounds(&self) -> &[Bound] {
        &self.bounds
    }
    fn objective(&self, x: &[f64]) -> f64 {
        x.iter().map(|v| v * v).sum()
    }
}

/// Rosenbrock: `f(x) = Σ [100(x_{i+1} − xᵢ²)² + (1 − xᵢ)²]`. A long, curved,
/// nearly-flat valley; minimum `0` at the all-ones vector.
pub struct Rosenbrock {
    bounds: Vec<Bound>,
}

impl Rosenbrock {
    /// `dim`-dimensional Rosenbrock (`dim >= 2`) on `[-5, 10]^dim`.
    pub fn new(dim: usize) -> Self {
        assert!(dim >= 2, "Rosenbrock needs at least 2 dimensions");
        Rosenbrock {
            bounds: vec![(-5.0, 10.0); dim],
        }
    }
}

impl Problem for Rosenbrock {
    fn dim(&self) -> usize {
        self.bounds.len()
    }
    fn bounds(&self) -> &[Bound] {
        &self.bounds
    }
    fn objective(&self, x: &[f64]) -> f64 {
        x.windows(2)
            .map(|w| 100.0 * (w[1] - w[0] * w[0]).powi(2) + (1.0 - w[0]).powi(2))
            .sum()
    }
}

/// Rastrigin: `f(x) = 10n + Σ [xᵢ² − 10 cos(2π xᵢ)]`. Highly multimodal with a
/// regular lattice of local minima; global minimum `0` at the origin.
pub struct Rastrigin {
    bounds: Vec<Bound>,
}

impl Rastrigin {
    /// `dim`-dimensional Rastrigin on `[-5.12, 5.12]^dim`.
    pub fn new(dim: usize) -> Self {
        Rastrigin {
            bounds: vec![(-5.12, 5.12); dim],
        }
    }
}

impl Problem for Rastrigin {
    fn dim(&self) -> usize {
        self.bounds.len()
    }
    fn bounds(&self) -> &[Bound] {
        &self.bounds
    }
    fn objective(&self, x: &[f64]) -> f64 {
        let n = x.len() as f64;
        10.0 * n
            + x.iter()
                .map(|v| v * v - 10.0 * (2.0 * std::f64::consts::PI * v).cos())
                .sum::<f64>()
    }
}

/// ZDT1 (Zitzler, Deb & Thiele 2000): a two-objective benchmark with a known,
/// convex Pareto front. On `[0, 1]^n`,
///
/// ```text
/// f1 = x1
/// g  = 1 + 9 · (Σ_{i≥2} xi) / (n − 1)
/// f2 = g · (1 − √(f1 / g))
/// ```
///
/// The Pareto-optimal set has `x2 … xn = 0` (so `g = 1`), giving the analytical
/// front `f2 = 1 − √f1` for `f1 ∈ [0, 1]` — the ground truth the NSGA-II test
/// measures distance to.
pub struct Zdt1 {
    bounds: Vec<Bound>,
}

impl Zdt1 {
    /// `dim`-dimensional ZDT1 (`dim >= 2`) on `[0, 1]^dim`.
    pub fn new(dim: usize) -> Self {
        assert!(dim >= 2, "ZDT1 needs at least 2 dimensions");
        Zdt1 {
            bounds: vec![(0.0, 1.0); dim],
        }
    }

    /// The analytical Pareto front value `f2 = 1 − √f1` for a given `f1`.
    pub fn front_f2(f1: f64) -> f64 {
        1.0 - f1.sqrt()
    }
}

impl MultiProblem for Zdt1 {
    fn dim(&self) -> usize {
        self.bounds.len()
    }
    fn bounds(&self) -> &[Bound] {
        &self.bounds
    }
    fn n_objectives(&self) -> usize {
        2
    }
    fn objectives(&self, x: &[f64]) -> Vec<f64> {
        let n = x.len();
        let f1 = x[0];
        let g = 1.0 + 9.0 * x[1..].iter().sum::<f64>() / (n - 1) as f64;
        let f2 = g * (1.0 - (f1 / g).sqrt());
        vec![f1, f2]
    }
}

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

    #[test]
    fn known_minima_are_zero() {
        assert_eq!(Sphere::new(3).objective(&[0.0, 0.0, 0.0]), 0.0);
        assert!(Rosenbrock::new(3).objective(&[1.0, 1.0, 1.0]).abs() < 1e-12);
        assert!(Rastrigin::new(4).objective(&[0.0; 4]).abs() < 1e-12);
    }

    #[test]
    fn zdt1_pareto_set_lies_on_the_analytical_front() {
        // x2..xn = 0 ⇒ g = 1 ⇒ f2 = 1 − √f1 for any f1 = x1.
        let p = Zdt1::new(5);
        for &f1 in &[0.0, 0.25, 0.5, 1.0] {
            let mut x = vec![0.0; 5];
            x[0] = f1;
            let o = p.objectives(&x);
            assert!((o[0] - f1).abs() < 1e-12);
            assert!((o[1] - Zdt1::front_f2(f1)).abs() < 1e-12);
        }
    }
}