diff --git a/src/operators/mod.rs b/src/operators/mod.rs index f40dbda..faa390d 100644 --- a/src/operators/mod.rs +++ b/src/operators/mod.rs @@ -4,8 +4,10 @@ pub mod binary; pub mod composite; pub mod permutation; pub mod real; +pub mod repair; pub use binary::*; pub use composite::*; pub use permutation::*; pub use real::*; +pub use repair::*; diff --git a/src/operators/repair.rs b/src/operators/repair.rs new file mode 100644 index 0000000..9401aac --- /dev/null +++ b/src/operators/repair.rs @@ -0,0 +1,171 @@ +//! Repair operators: in-place projections that restore decisions to +//! feasibility. + +use crate::traits::Repair; + +/// Clamp every variable of a `Vec` to per-axis inclusive bounds. +/// +/// The simplest possible repair — pair with `GaussianMutation` (which +/// doesn't enforce bounds in v1) to produce a bounds-respecting variant +/// without writing a custom Variation impl. +#[derive(Debug, Clone)] +pub struct ClampToBounds { + /// Per-variable inclusive bounds. + pub bounds: Vec<(f64, f64)>, +} + +impl ClampToBounds { + /// Construct a `ClampToBounds`. + /// + /// # Panics + /// If any `(lo, hi)` has `lo > hi`. + pub fn new(bounds: Vec<(f64, f64)>) -> Self { + for (i, &(lo, hi)) in bounds.iter().enumerate() { + assert!( + lo <= hi, + "ClampToBounds bound at index {i} has lo > hi: ({lo}, {hi})", + ); + } + Self { bounds } + } +} + +impl Repair> for ClampToBounds { + fn repair(&mut self, decision: &mut Vec) { + for (j, x) in decision.iter_mut().enumerate() { + if let Some(&(lo, hi)) = self.bounds.get(j) { + *x = x.clamp(lo, hi); + } + } + } +} + +/// Project a `Vec` onto the simplex `{ x : x ≥ 0, Σ x = total }`. +/// +/// Implements the standard O(n log n) projection algorithm of Wang & Carreira- +/// Perpiñán 2013. Useful for portfolio-style problems where the +/// decision must sum to a budget, and for normalizing reference +/// directions onto the unit simplex. +#[derive(Debug, Clone)] +pub struct ProjectToSimplex { + /// Target sum (the simplex's "size"). Standard probability simplex + /// uses `total = 1.0`. + pub total: f64, +} + +impl ProjectToSimplex { + /// Construct a `ProjectToSimplex`. + /// + /// # Panics + /// If `total <= 0.0`. + pub fn new(total: f64) -> Self { + assert!(total > 0.0, "ProjectToSimplex total must be > 0"); + Self { total } + } +} + +impl Repair> for ProjectToSimplex { + fn repair(&mut self, decision: &mut Vec) { + let n = decision.len(); + if n == 0 { + return; + } + // Sort copy descending. + let mut sorted: Vec = decision.clone(); + sorted.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal)); + + // Find ρ = max{ j : sorted[j-1] - (Σ_{i<=j} sorted[i] - total) / j > 0 }. + let mut cumsum = 0.0; + let mut rho = 0; + let mut tau_at_rho = 0.0; + for (j, &val) in sorted.iter().enumerate() { + cumsum += val; + let tau = (cumsum - self.total) / (j as f64 + 1.0); + if val - tau > 0.0 { + rho = j + 1; + tau_at_rho = tau; + } + } + let _ = rho; + // Apply: x_i ← max(x_i - τ, 0). + for x in decision.iter_mut() { + *x = (*x - tau_at_rho).max(0.0); + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + fn approx_eq(a: f64, b: f64, tol: f64) -> bool { + (a - b).abs() < tol + } + + #[test] + fn clamp_to_bounds_clips() { + let mut r = ClampToBounds::new(vec![(-1.0, 1.0); 3]); + let mut x = vec![-2.5, 0.5, 5.0]; + r.repair(&mut x); + assert_eq!(x, vec![-1.0, 0.5, 1.0]); + } + + #[test] + fn clamp_passthrough_when_already_in_bounds() { + let mut r = ClampToBounds::new(vec![(-1.0, 1.0); 3]); + let mut x = vec![-0.3, 0.0, 0.7]; + let original = x.clone(); + r.repair(&mut x); + assert_eq!(x, original); + } + + #[test] + #[should_panic(expected = "lo > hi")] + fn clamp_invalid_bounds_panics() { + let _ = ClampToBounds::new(vec![(1.0, -1.0)]); + } + + #[test] + fn project_to_unit_simplex_sums_to_total() { + let mut r = ProjectToSimplex::new(1.0); + let mut x = vec![0.5, 0.3, 0.2, -0.5]; + r.repair(&mut x); + let s: f64 = x.iter().sum(); + assert!(approx_eq(s, 1.0, 1e-12)); + for &v in &x { + assert!(v >= 0.0); + } + } + + #[test] + fn project_already_on_simplex_unchanged() { + let mut r = ProjectToSimplex::new(1.0); + let mut x = vec![0.5, 0.3, 0.2]; + r.repair(&mut x); + let s: f64 = x.iter().sum(); + assert!(approx_eq(s, 1.0, 1e-12)); + // Within tolerance, the values should be roughly preserved (no + // clipping needed). + assert!(approx_eq(x[0], 0.5, 1e-12)); + assert!(approx_eq(x[1], 0.3, 1e-12)); + assert!(approx_eq(x[2], 0.2, 1e-12)); + } + + #[test] + fn project_arbitrary_total() { + let mut r = ProjectToSimplex::new(10.0); + let mut x = vec![100.0, 50.0, -20.0, 30.0]; + r.repair(&mut x); + let s: f64 = x.iter().sum(); + assert!(approx_eq(s, 10.0, 1e-9)); + for &v in &x { + assert!(v >= 0.0); + } + } + + #[test] + #[should_panic(expected = "total must be > 0")] + fn project_non_positive_total_panics() { + let _ = ProjectToSimplex::new(0.0); + } +} diff --git a/src/prelude.rs b/src/prelude.rs index c980688..416c653 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -9,7 +9,7 @@ pub use crate::core::{ Population, Problem, Rng, rng_from_seed, }; -pub use crate::traits::{Initializer, Optimizer, Variation}; +pub use crate::traits::{Initializer, Optimizer, Repair, Variation}; pub use crate::pareto::{ Dominance, ParetoArchive, best_candidate, crowding_distance, das_dennis, @@ -17,8 +17,9 @@ pub use crate::pareto::{ }; pub use crate::operators::{ - BitFlipMutation, BoundedGaussianMutation, CompositeVariation, GaussianMutation, - LevyMutation, PolynomialMutation, RealBounds, SimulatedBinaryCrossover, SwapMutation, + BitFlipMutation, BoundedGaussianMutation, ClampToBounds, CompositeVariation, + GaussianMutation, LevyMutation, PolynomialMutation, ProjectToSimplex, RealBounds, + SimulatedBinaryCrossover, SwapMutation, }; pub use crate::algorithms::{ diff --git a/src/traits/mod.rs b/src/traits/mod.rs index faeb203..0829e19 100644 --- a/src/traits/mod.rs +++ b/src/traits/mod.rs @@ -2,8 +2,10 @@ pub mod initializer; pub mod optimizer; +pub mod repair; pub mod variation; pub use initializer::*; pub use optimizer::*; +pub use repair::*; pub use variation::*; diff --git a/src/traits/repair.rs b/src/traits/repair.rs new file mode 100644 index 0000000..07852ee --- /dev/null +++ b/src/traits/repair.rs @@ -0,0 +1,14 @@ +//! Trait for restoring decisions to feasibility. + +/// Transforms an infeasible (or possibly-infeasible) decision into a +/// feasible one, in place. +/// +/// `Repair` is the projection-style alternative to penalty-style +/// constraint handling (which uses `Evaluation::constraint_violation`). +/// Wrap a `Variation` operator's output through a `Repair` to guarantee +/// feasibility; or call `repair()` inside your `Problem::evaluate` if +/// the constraint structure is best handled at evaluation time. +pub trait Repair { + /// Mutate `decision` in place to satisfy the repair's constraints. + fn repair(&mut self, decision: &mut D); +}