radiate-alters 1.2.22

use radiate_core::{BoundedGene, Chromosome, FloatGene, Gene, Mutate, Rate, random_provider};
use radiate_utils::{Float, Primitive};

// Use it when:
// 	- You’re evolving floating-point representations (like real-valued neural nets, control parameters, orbital mechanics).
// 	- You want bounded and unbiased mutation behavior.
// 	- You care about the distribution of the mutations. Unlike Gaussian, polynomial gives more control over tail behavior.

// It’s especially powerful in:
// 	- Multi-objective algorithms (like NSGA-II)
// 	- High-precision tuning problems
// 	- Bounded domains where classic Gaussian mutation may overshoot

// The eta parameter in polynomial mutation controls the shape of the mutation distribution.
// In other words, how local or global your mutations are:
// 	- eta is the distribution index (usually denoted as η_m in literature like Deb’s NSGA-II paper).
// 	- It determines the exploration vs. exploitation trade-off:
// 	- Low eta (e.g. 1–5): leads to bigger mutations, promoting exploration.
// 	- High eta (e.g. 20–100): leads to smaller, fine-grained mutations, good for local search.
#[derive(Debug)]
pub struct PolynomialMutator {
    rate: Rate,
    eta: f32,
}

impl PolynomialMutator {
    pub fn new(rate: impl Into<Rate>, eta: f32) -> Self {
        let rate = rate.into();
        PolynomialMutator { rate, eta }
    }

    fn polynomial_mutation(&self, value: f64, min: f64, max: f64, eta: f64) -> f64 {
        let u = random_provider::random::<f64>();

        if (max - min).abs() < f64::EPSILON {
            return value;
        }

        let delta1 = (value - min) / (max - min);
        let delta2 = (max - value) / (max - min);

        let mutq = if u <= 0.5 {
            // Left side of the polynomial
            let term1 = 2.0 * u;
            let term2 = (1.0 - 2.0 * u) * (1.0 - delta1).powf(eta + 1.0);
            (term1 + term2).powf(1.0 / (eta + 1.0))
        } else {
            // Right side of the polynomial
            let term1 = 2.0 * (1.0 - u);
            let term2 = 2.0 * (u - 0.5) * (1.0 - delta2).powf(eta + 1.0);
            1.0 - (term1 + term2).powf(1.0 / (eta + 1.0))
        };

        min + mutq * (max - min)
    }
}

impl<F, C> Mutate<C> for PolynomialMutator
where
    F: Float + Primitive,
    C: Chromosome<Gene = FloatGene<F>>,
{
    fn rate(&self) -> Rate {
        self.rate.clone()
    }

    #[inline]
    fn mutate_gene(&self, gene: &mut C::Gene) -> usize {
        // TODO: Should these be from the bounds?
        let min = gene.min().extract::<f64>().unwrap();
        let max = gene.max().extract::<f64>().unwrap();
        let value = gene.allele().extract::<f64>().unwrap();
        let eta = self.eta as f64;

        let new_value = self.polynomial_mutation(value, min, max, eta);

        let clamped_value = new_value.clamp(min, max);
        *gene.allele_mut() = clamped_value.extract::<F>().unwrap();
        1
    }
}