feat(operators): add PolynomialMutation

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

@@ -157,6 +157,84 @@ impl Variation<Vec<f64>> 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<Vec<f64>> for PolynomialMutation {
+    fn vary(&mut self, parents: &[Vec<f64>], rng: &mut Rng) -> Vec<Vec<f64>> {
+        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);
+    }
 }

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::{