forge_core/

constraint.rs

//! Nonlinear constraint handling.
//!
//! forge's optimizers are box-constrained by design; real calibration and
//! design problems add **nonlinear** constraints (parameter relationships in
//! hydrological models, adjacency/budget rules in spatial planning). This
//! module provides the two standard parameter-free techniques from the CEC
//! constrained-optimization lineage as **adapters**: a [`ConstrainedProblem`]
//! is converted into an ordinary [`Problem`], so *every* forge optimizer —
//! DE, L-SHADE, L-SRTDE, CMA-ES, DDS, SCE-UA, PSO, restarts, ensembles —
//! handles constraints without modification.
//!
//! - [`DebRules`] — Deb's feasibility rules (Deb 2000, CMAME 186:311–338):
//!   feasible beats infeasible; feasible candidates compare by objective;
//!   infeasible candidates compare by total violation. Encoded exactly as
//!   Deb's original fitness transformation: an infeasible point scores
//!   `FEASIBLE_CEILING + total_violation`, above every feasible objective.
//! - [`EpsilonLShade`](crate::algo::EpsilonLShade) — the ε-constrained
//!   method (Takahama & Sakai lineage, the backbone of CEC winners such as
//!   LSHADE44 and εMAg-ES) lives *inside* the algorithm, not here: violations
//!   up to a shrinking tolerance ε(t) count as feasible, and the comparison
//!   is re-applied with the **current** ε each generation. An adapter cannot
//!   express that faithfully — optimizers cache fitness values, so a point
//!   accepted under an early, generous ε would keep its stale score after ε
//!   contracts and block genuinely feasible solutions.
//!
//! **Caveat:** the [`DebRules`] transformation is exact as long as every
//! feasible objective value is below [`FEASIBLE_CEILING`] (`1e100`).

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

/// Objective values at or above this are reserved for infeasible candidates
/// (Deb's transformation): `score = FEASIBLE_CEILING · (1 + violation)`.
/// The encoding is multiplicative because an additive `1e100 + violation`
/// would be absorbed by the ulp of `1e100` (~1.9e84) and erase the ranking
/// between different violations.
pub const FEASIBLE_CEILING: f64 = 1e100;

/// Deb-style score for an infeasible candidate: always ≥ [`FEASIBLE_CEILING`]
/// and strictly increasing in the violation (down to ~1e-16 relative).
#[inline]
fn infeasible_score(violation: f64) -> f64 {
    FEASIBLE_CEILING * (1.0 + violation)
}

/// A box-constrained problem with additional nonlinear constraints.
///
/// Report constraints through [`violations`](ConstrainedProblem::violations):
/// one entry per constraint, `0.0` when satisfied and the (positive)
/// magnitude of the breach otherwise. For an inequality `g(x) ≤ 0` return
/// `g(x).max(0.0)`; for an equality `h(x) = 0` with tolerance `δ` return
/// `(h(x).abs() − δ).max(0.0)`.
pub trait ConstrainedProblem {
    /// Number of decision variables.
    fn dim(&self) -> usize;

    /// Per-variable inclusive bounds, length [`ConstrainedProblem::dim`].
    fn bounds(&self) -> &[Bound];

    /// Objective value to **minimize** (non-finite ⇒ infeasible, as usual).
    fn objective(&self, x: &[f64]) -> f64;

    /// Constraint violations: `0.0` per satisfied constraint, positive
    /// magnitudes otherwise. Negative entries are clamped to `0.0`.
    fn violations(&self, x: &[f64]) -> Vec<f64>;

    /// Total violation `Σ max(vᵢ, 0)` — the quantity Deb's rules compare.
    fn total_violation(&self, x: &[f64]) -> f64 {
        self.violations(x).iter().map(|v| v.max(0.0)).sum()
    }
}

/// Deb's feasibility rules (Deb 2000) as a [`Problem`] adapter.
///
/// Ranking induced: feasible by objective, then infeasible by total
/// violation — every feasible point beats every infeasible one. Works with
/// any forge optimizer; exact whenever feasible objectives stay below
/// [`FEASIBLE_CEILING`].
///
/// ```
/// use forge_core::constraint::{constrained_func, DebRules};
/// use forge_core::{De, Optimizer, Termination};
///
/// // Minimize x + y subject to x·y ≥ 1 on [0, 10]².
/// let p = constrained_func(
///     vec![(0.0, 10.0); 2],
///     |x| x[0] + x[1],
///     |x| vec![(1.0 - x[0] * x[1]).max(0.0)],
/// );
/// let report = De::default().optimize(&DebRules(p), &Termination::budget(6000));
/// assert!((report.best_value() - 2.0).abs() < 1e-2); // optimum at (1, 1)
/// ```
pub struct DebRules<P>(pub P);

impl<P: ConstrainedProblem> Problem for DebRules<P> {
    fn dim(&self) -> usize {
        self.0.dim()
    }
    fn bounds(&self) -> &[Bound] {
        self.0.bounds()
    }
    fn objective(&self, x: &[f64]) -> f64 {
        let violation = self.0.total_violation(x);
        if violation > 0.0 {
            // Infeasible: rank by violation only, above every feasible value.
            infeasible_score(violation)
        } else {
            self.0.objective(x)
        }
    }
}

/// A [`ConstrainedProblem`] defined inline by closures.
pub fn constrained_func<F, G>(bounds: Vec<Bound>, f: F, g: G) -> ConstrainedFunc<F, G>
where
    F: Fn(&[f64]) -> f64,
    G: Fn(&[f64]) -> Vec<f64>,
{
    ConstrainedFunc { bounds, f, g }
}

/// Closure-backed constrained problem produced by [`constrained_func`].
pub struct ConstrainedFunc<F, G> {
    bounds: Vec<Bound>,
    f: F,
    g: G,
}

impl<F, G> ConstrainedProblem for ConstrainedFunc<F, G>
where
    F: Fn(&[f64]) -> f64,
    G: Fn(&[f64]) -> Vec<f64>,
{
    fn dim(&self) -> usize {
        self.bounds.len()
    }
    fn bounds(&self) -> &[Bound] {
        &self.bounds
    }
    fn objective(&self, x: &[f64]) -> f64 {
        (self.f)(x)
    }
    fn violations(&self, x: &[f64]) -> Vec<f64> {
        (self.g)(x)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::algo::{LShade, Optimizer};
    use crate::termination::Termination;

    /// Minimize x² + y² subject to x + y ≥ 1: optimum at (0.5, 0.5), f = 0.5.
    #[allow(clippy::type_complexity)]
    fn ring() -> ConstrainedFunc<impl Fn(&[f64]) -> f64, impl Fn(&[f64]) -> Vec<f64>> {
        constrained_func(
            vec![(-2.0, 2.0); 2],
            |x| x[0] * x[0] + x[1] * x[1],
            |x| vec![(1.0 - x[0] - x[1]).max(0.0)],
        )
    }

    #[test]
    fn deb_rules_find_the_constrained_optimum() {
        let report = LShade::default().optimize(&DebRules(ring()), &Termination::budget(8000));
        assert!(
            (report.best_value() - 0.5).abs() < 1e-3,
            "got {}",
            report.best_value()
        );
        assert!((report.best()[0] - 0.5).abs() < 0.03);
        assert!((report.best()[1] - 0.5).abs() < 0.03);
    }

    #[test]
    fn infeasible_ranks_by_violation_and_below_feasible() {
        let p = DebRules(ring());
        let feasible = p.objective(&[1.0, 1.0]); // violation 0, f = 2
        let mildly_infeasible = p.objective(&[0.4, 0.4]); // violation 0.2
        let badly_infeasible = p.objective(&[-1.0, -1.0]); // violation 3
        assert!(feasible < mildly_infeasible);
        assert!(mildly_infeasible < badly_infeasible);
        assert!(mildly_infeasible >= FEASIBLE_CEILING);
    }
}