From 464617b1f8e642ff3751b12ed6837d6dd4ecad5d Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Mon, 4 May 2026 19:44:05 -0600 Subject: [PATCH] feat(operators): add PolynomialMutation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Deb's standard real-valued mutation pair to SBX, used together by canonical NSGA-II. For each variable, with probability `per_variable_probability` (typical: 1/n where n is dim), perturb the parent value by a polynomial-distributed delta scaled by the bound range, then clamp. Per-dim formula: - `u ~ U[0, 1)` - `δ = (2u)^(1/(η+1)) − 1` if `u < 0.5` else `1 − (2(1−u))^(1/(η+1))` - `child[j] = parent[j] + δ · (hi − lo)`, clamped to bounds `eta` is the distribution index (typical 20; smaller → more spread). This is the simple bound-rescale form; the bound-aware δ_q variant from the full paper is left as a future refinement. Always returns one child. Tests cover: child stays in bounds with high sigma-equivalent eta, per_variable_probability=0 returns the parent unchanged, and standard panics. --- src/operators/real.rs | 110 ++++++++++++++++++++++++++++++++++++++++++ src/prelude.rs | 4 +- 2 files changed, 112 insertions(+), 2 deletions(-) diff --git a/src/operators/real.rs b/src/operators/real.rs index c703827..3366fb4 100644 --- a/src/operators/real.rs +++ b/src/operators/real.rs @@ -157,6 +157,84 @@ impl Variation> for SimulatedBinaryCrossover { } } +/// Deb's polynomial mutation — the canonical NSGA-II mutation operator and +/// the standard pair to [`SimulatedBinaryCrossover`]. +/// +/// Per dimension, with probability `per_variable_probability`, perturb the +/// parent's value by a polynomial-distributed delta scaled by the bound +/// range: +/// +/// - `u ~ U[0, 1)` +/// - `δ = (2u)^(1/(η+1)) − 1` if `u < 0.5` else `1 − (2(1−u))^(1/(η+1))` +/// - `child[j] = parent[j] + δ · (hi − lo)`, clamped to `[lo, hi]` +/// +/// `eta` is the distribution index (typical 20; smaller → more spread). +/// The convention for `per_variable_probability` is `1.0 / dim`. +/// +/// This is the simple bound-rescale form; the bound-aware `δ_q` variant from +/// the full paper is left as a future refinement. +#[derive(Debug, Clone)] +pub struct PolynomialMutation { + /// Per-variable inclusive bounds. Length must match the parent decision. + pub bounds: Vec<(f64, f64)>, + /// Distribution index `η_m`. Must be `>= 0.0`. Default 20.0. + pub eta: f64, + /// Per-variable mutation probability. Conventional: `1.0 / dim`. + pub per_variable_probability: f64, +} + +impl PolynomialMutation { + /// Construct a `PolynomialMutation`. + /// + /// # Panics + /// If any bound has `lo > hi`, `eta < 0.0`, or + /// `per_variable_probability` is outside `[0.0, 1.0]`. + pub fn new(bounds: Vec<(f64, f64)>, eta: f64, per_variable_probability: f64) -> Self { + for (i, &(lo, hi)) in bounds.iter().enumerate() { + assert!( + lo <= hi, + "PolynomialMutation bound at index {i} has lo > hi: ({lo}, {hi})", + ); + } + assert!(eta >= 0.0, "PolynomialMutation eta must be >= 0.0"); + assert!( + (0.0..=1.0).contains(&per_variable_probability), + "PolynomialMutation per_variable_probability must be in [0.0, 1.0]", + ); + Self { bounds, eta, per_variable_probability } + } +} + +impl Variation> for PolynomialMutation { + fn vary(&mut self, parents: &[Vec], rng: &mut Rng) -> Vec> { + assert!( + !parents.is_empty(), + "PolynomialMutation requires at least one parent", + ); + assert_eq!( + parents[0].len(), + self.bounds.len(), + "PolynomialMutation parent length must match bounds length", + ); + let exponent = 1.0 / (self.eta + 1.0); + let mut child = parents[0].clone(); + for j in 0..self.bounds.len() { + if !rng.random_bool(self.per_variable_probability) { + continue; + } + let u: f64 = rng.random(); + let delta = if u < 0.5 { + (2.0 * u).powf(exponent) - 1.0 + } else { + 1.0 - (2.0 * (1.0 - u)).powf(exponent) + }; + let (lo, hi) = self.bounds[j]; + child[j] = (child[j] + delta * (hi - lo)).clamp(lo, hi); + } + vec![child] + } +} + /// Bounded variant of [`GaussianMutation`]: add `Normal(0, sigma)` noise to /// every variable of the first parent, then clamp each variable to its /// per-dimension inclusive bound. @@ -346,4 +424,36 @@ mod tests { fn sbx_negative_eta_panics() { let _ = SimulatedBinaryCrossover::new(vec![(0.0, 1.0)], -1.0, 0.5); } + + #[test] + fn polynomial_mutation_keeps_child_in_bounds() { + let mut m = PolynomialMutation::new(vec![(-1.0, 1.0); 5], 5.0, 1.0); + let mut rng = rng_from_seed(99); + let parent = vec![0.0_f64; 5]; + for _ in 0..100 { + let children = m.vary(std::slice::from_ref(&parent), &mut rng); + assert_eq!(children.len(), 1); + assert_eq!(children[0].len(), 5); + for &x in &children[0] { + assert!(x >= -1.0 && x <= 1.0, "out of bounds: {x}"); + } + } + } + + #[test] + fn polynomial_mutation_zero_probability_returns_parent() { + let mut m = PolynomialMutation::new(vec![(-10.0, 10.0); 3], 20.0, 0.0); + let mut rng = rng_from_seed(0); + let parent = vec![1.0, -2.0, 3.0]; + let children = m.vary(std::slice::from_ref(&parent), &mut rng); + assert_eq!(children[0], parent); + } + + #[test] + #[should_panic(expected = "must match bounds length")] + fn polynomial_mutation_mismatched_length_panics() { + let mut m = PolynomialMutation::new(vec![(0.0, 1.0); 3], 20.0, 0.1); + let mut rng = rng_from_seed(0); + m.vary(&[vec![0.5; 2]], &mut rng); + } } diff --git a/src/prelude.rs b/src/prelude.rs index 1ce45b9..c3dd7a7 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -17,8 +17,8 @@ pub use crate::pareto::{ }; pub use crate::operators::{ - BitFlipMutation, BoundedGaussianMutation, GaussianMutation, RealBounds, - SimulatedBinaryCrossover, SwapMutation, + BitFlipMutation, BoundedGaussianMutation, GaussianMutation, PolynomialMutation, + RealBounds, SimulatedBinaryCrossover, SwapMutation, }; pub use crate::algorithms::{