From 1b8070476b6e3f60cae742046f5922b08a2dff43 Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Tue, 5 May 2026 09:06:28 -0600 Subject: [PATCH] feat(operators): add LevyMutation real-valued heavy-tailed mutation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Lévy-flight perturbation: each variable receives a step drawn from a heavy-tailed Lévy(α) distribution rather than a Normal. The result is "mostly small steps with rare big jumps," which gives a more exploratory mutation than Gaussian without abandoning local search. Decision type: Vec, with optional bounds (clamped per-axis if `bounds` is non-empty). The step is sampled via Mantegna's algorithm which generates Lévy(α) by combining two Normal samples and taking the right power, controlled by the tail exponent `alpha` (typical 1.5; 1 is heavy, 2 collapses to Normal). This is the only genuinely-different mutation kernel from Cuckoo Search and other Lévy-flight metaheuristics; ship it as a Variation operator usable from any algorithm rather than as a separate algorithm. --- src/operators/real.rs | 141 ++++++++++++++++++++++++++++++++++++++++++ src/prelude.rs | 2 +- 2 files changed, 142 insertions(+), 1 deletion(-) diff --git a/src/operators/real.rs b/src/operators/real.rs index 37af42e..514a5ca 100644 --- a/src/operators/real.rs +++ b/src/operators/real.rs @@ -289,6 +289,111 @@ impl Variation> for BoundedGaussianMutation { } } +/// Heavy-tailed Lévy-flight mutation for `Vec` decisions. +/// +/// Adds a Lévy(α)-distributed step to every variable, optionally clamped to +/// per-variable bounds. Compared with `GaussianMutation`, the Lévy +/// distribution has a heavy tail — most steps are small and local but +/// occasional steps are very large, giving a single mutation operator +/// that does both refinement and exploration. This is the kernel that +/// powers Cuckoo Search and other Lévy-flight metaheuristics. +/// +/// Implementation: Mantegna's algorithm combines two Gaussians to +/// produce a Lévy(α) sample. `alpha` is the tail exponent in `(0, 2]`; +/// typical value is `1.5`. `1.0` gives the Cauchy distribution (very +/// heavy); `2.0` collapses to the Normal. +#[derive(Debug, Clone)] +pub struct LevyMutation { + /// Tail exponent `α ∈ (0, 2]`. Smaller = heavier tail. + pub alpha: f64, + /// Step scale. + pub scale: f64, + /// Optional per-variable bounds. Empty `Vec` → no clamping. + pub bounds: Vec<(f64, f64)>, +} + +impl LevyMutation { + /// Construct a `LevyMutation`. + /// + /// # Panics + /// If `alpha` is not in `(0, 2]`, `scale <= 0.0`, or any bound has + /// `lo > hi`. + pub fn new(alpha: f64, scale: f64, bounds: Vec<(f64, f64)>) -> Self { + assert!( + alpha > 0.0 && alpha <= 2.0, + "LevyMutation alpha must be in (0, 2]", + ); + assert!(scale > 0.0, "LevyMutation scale must be > 0"); + for (i, &(lo, hi)) in bounds.iter().enumerate() { + assert!( + lo <= hi, + "LevyMutation bound at index {i} has lo > hi: ({lo}, {hi})", + ); + } + Self { alpha, scale, bounds } + } +} + +impl Variation> for LevyMutation { + fn vary(&mut self, parents: &[Vec], rng: &mut Rng) -> Vec> { + assert!(!parents.is_empty(), "LevyMutation requires at least one parent"); + let alpha = self.alpha; + // Mantegna's algorithm σ for the numerator Normal: + // sigma_u = (Γ(1+α)·sin(π·α/2) / (Γ((1+α)/2)·α·2^((α-1)/2)))^(1/α) + // Denominator Normal has σ = 1. + let sigma_u = mantegna_sigma_u(alpha); + let normal_u = Normal::new(0.0, sigma_u).expect("Normal::new(0, sigma_u)"); + let normal_v = Normal::new(0.0, 1.0).expect("Normal::new(0, 1)"); + let mut child = parents[0].clone(); + for (j, x) in child.iter_mut().enumerate() { + let u: f64 = normal_u.sample(rng); + let v: f64 = normal_v.sample(rng); + let step = u / v.abs().powf(1.0 / alpha); + *x += self.scale * step; + if let Some(&(lo, hi)) = self.bounds.get(j) { + *x = x.clamp(lo, hi); + } + } + vec![child] + } +} + +fn mantegna_sigma_u(alpha: f64) -> f64 { + // Γ-related constants. We compute Γ(z) via libm if the std::f64::gamma + // isn't available; fall back to a small Lanczos approximation. + fn gamma(z: f64) -> f64 { + // Stirling-ish via the standard recursion + Lanczos coefficients. + // For the typical α ∈ [1, 2] range we hit, the expressions Γ(1+α) + // and Γ((1+α)/2) are well-behaved. + let g = 7.0; + let p = [ + 0.999_999_999_999_809_93, + 676.520_368_121_885_1, + -1_259.139_216_722_4023, + 771.323_428_777_653_13, + -176.615_029_162_140_59, + 12.507_343_278_686_905, + -0.138_571_095_265_720_12, + 9.984_369_578_019_571_6e-6, + 1.505_632_735_149_311_6e-7, + ]; + if z < 0.5 { + std::f64::consts::PI / ((std::f64::consts::PI * z).sin() * gamma(1.0 - z)) + } else { + let z = z - 1.0; + let mut x = p[0]; + for (i, &pi) in p.iter().enumerate().skip(1) { + x += pi / (z + i as f64); + } + let t = z + g + 0.5; + (2.0 * std::f64::consts::PI).sqrt() * t.powf(z + 0.5) * (-t).exp() * x + } + } + let num = gamma(1.0 + alpha) * (std::f64::consts::PI * alpha / 2.0).sin(); + let den = gamma((1.0 + alpha) / 2.0) * alpha * 2.0_f64.powf((alpha - 1.0) / 2.0); + (num / den).powf(1.0 / alpha) +} + #[cfg(test)] mod tests { use super::*; @@ -426,6 +531,42 @@ mod tests { let _ = SimulatedBinaryCrossover::new(vec![(0.0, 1.0)], -1.0, 0.5); } + #[test] + fn levy_mutation_returns_one_child_in_bounds() { + let mut m = LevyMutation::new(1.5, 0.1, vec![(-1.0, 1.0); 4]); + let mut rng = rng_from_seed(42); + let parent = vec![0.0_f64; 4]; + for _ in 0..50 { + let children = m.vary(std::slice::from_ref(&parent), &mut rng); + assert_eq!(children.len(), 1); + assert_eq!(children[0].len(), 4); + for &x in &children[0] { + assert!((-1.0..=1.0).contains(&x), "out of bounds: {x}"); + } + } + } + + #[test] + fn levy_mutation_unbounded_works() { + let mut m = LevyMutation::new(1.5, 0.5, Vec::new()); + let mut rng = rng_from_seed(0); + let parent = vec![0.0_f64; 3]; + let children = m.vary(std::slice::from_ref(&parent), &mut rng); + assert_eq!(children[0].len(), 3); + } + + #[test] + #[should_panic(expected = "alpha must be in (0, 2]")] + fn levy_alpha_out_of_range_panics() { + let _ = LevyMutation::new(0.0, 0.1, Vec::new()); + } + + #[test] + #[should_panic(expected = "scale must be > 0")] + fn levy_zero_scale_panics() { + let _ = LevyMutation::new(1.5, 0.0, Vec::new()); + } + #[test] fn polynomial_mutation_keeps_child_in_bounds() { let mut m = PolynomialMutation::new(vec![(-1.0, 1.0); 5], 5.0, 1.0); diff --git a/src/prelude.rs b/src/prelude.rs index 44c7a0e..0c89970 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -18,7 +18,7 @@ pub use crate::pareto::{ pub use crate::operators::{ BitFlipMutation, BoundedGaussianMutation, CompositeVariation, GaussianMutation, - PolynomialMutation, RealBounds, SimulatedBinaryCrossover, SwapMutation, + LevyMutation, PolynomialMutation, RealBounds, SimulatedBinaryCrossover, SwapMutation, }; pub use crate::algorithms::{