metalforge 0.3.0

//! MOEA/D — Multi-objective Evolutionary Algorithm based on Decomposition
//! (Zhang & Li 2007, IEEE Trans. Evol. Comput. 11(6), 712–731).
//!
//! The third canonical MOEA family, next to dominance-based
//! ([`NsgaII`](super::NsgaII)/[`NsgaIII`](super::NsgaIII)) and
//! indicator-based ([`SmsEmoa`](super::SmsEmoa)): the multi-objective problem
//! is **decomposed** into `N` scalar subproblems — here by the paper's
//! recommended **Tchebycheff** scalarization
//! `g(x | w, z*) = maxⱼ wⱼ · |fⱼ(x) − z*ⱼ|` over Das–Dennis weight vectors —
//! and all subproblems are optimized *simultaneously*, each using only its
//! `T` nearest-neighbor subproblems for mating and replacement. The ideal
//! point `z*` is the running per-objective minimum.
//!
//! This implementation follows the original 2007 algorithm with the two
//! standard robustness refinements from MOEA/D-DE (Li & Zhang 2009), which
//! pymoo's implementation also adopts:
//!
//! - with probability `1 − delta` the mating pool is the whole population
//!   instead of the neighborhood (escapes neighborhood-wide stagnation), and
//! - a replacement cap `nr` limits how many neighbors one offspring may
//!   overwrite (preserves diversity).
//!
//! Variation is SBX + polynomial mutation (shared with the NSGA family).
//! Returns the non-dominated set of the final population as a
//! [`ParetoFront`]. Deterministic for a given seed.

use super::nsga2::{fast_non_dominated_sort, mutate, sbx};
use super::nsga3::das_dennis;
use crate::problem::MultiProblem;
use crate::rng::Rng;
use crate::solution::{MultiSolution, ParetoFront};
use crate::termination::Termination;

/// MOEA/D configuration. The population size is the number of Das–Dennis
/// weight vectors, `C(M + divisions − 1, divisions)` for `M` objectives.
#[derive(Debug, Clone, Copy)]
pub struct Moead {
    /// Das–Dennis divisions for the weight vectors (99 ⇒ 100 subproblems for
    /// 2 objectives; 13 ⇒ 105 for 3).
    pub divisions: usize,
    /// Neighborhood size `T` (clamped to the population size).
    pub neighbors: usize,
    /// Probability of mating within the neighborhood (vs the whole
    /// population).
    pub delta: f64,
    /// Maximum number of neighbor solutions one offspring may replace.
    pub max_replacements: usize,
    /// SBX distribution index `η_c`.
    pub crossover_eta: f64,
    /// Probability of applying SBX to a parent pair.
    pub crossover_prob: f64,
    /// Polynomial-mutation distribution index `η_m`.
    pub mutation_eta: f64,
    /// Per-variable mutation probability; `None` uses `1/dim`.
    pub mutation_prob: Option<f64>,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for Moead {
    fn default() -> Self {
        Moead {
            divisions: 99,
            neighbors: 20,
            delta: 0.9,
            max_replacements: 2,
            crossover_eta: 20.0,
            crossover_prob: 1.0,
            mutation_eta: 20.0,
            mutation_prob: None,
            seed: 42,
        }
    }
}

impl Moead {
    /// Approximates the Pareto front of `problem` within the evaluation
    /// budget of `term` (its `target` is ignored).
    pub fn optimize(&self, problem: &dyn MultiProblem, term: &Termination) -> ParetoFront {
        crate::problem::validate_multi(problem)
            .unwrap_or_else(|e| panic!("Moead: invalid problem: {e}"));
        let bounds = problem.bounds();
        let dim = bounds.len();
        let m = problem.n_objectives();
        let pm = self.mutation_prob.unwrap_or(1.0 / dim as f64);
        let mut rng = Rng::new(self.seed);

        // Weight vectors define the subproblems; zero components are lifted
        // to a small positive value so the Tchebycheff max never ignores an
        // objective entirely (the usual 1e-6 repair).
        let weights: Vec<Vec<f64>> = das_dennis(m, self.divisions.max(1))
            .into_iter()
            .map(|w| w.into_iter().map(|v| v.max(1e-6)).collect())
            .collect();
        let n = weights.len();
        let t = self.neighbors.clamp(2, n);

        // B(i): indices of the T nearest weight vectors (Euclidean).
        let neighborhoods: Vec<Vec<usize>> = (0..n)
            .map(|i| {
                let mut order: Vec<usize> = (0..n).collect();
                order.sort_by(|&a, &b| {
                    let da = dist2(&weights[i], &weights[a]);
                    let db = dist2(&weights[i], &weights[b]);
                    da.partial_cmp(&db).unwrap_or(std::cmp::Ordering::Equal)
                });
                order.truncate(t);
                order
            })
            .collect();

        let eval = |x: &[f64]| -> Vec<f64> {
            let mut o = problem.objectives(x);
            for v in &mut o {
                if !v.is_finite() {
                    *v = f64::INFINITY;
                }
            }
            o
        };

        // Initial population: one solution per subproblem (the budget binds).
        let mut pop: Vec<MultiSolution> = Vec::with_capacity(n);
        let mut ideal = vec![f64::INFINITY; m];
        let mut evaluations = 0usize;
        for _ in 0..n {
            if evaluations >= term.max_evaluations {
                break;
            }
            let x: Vec<f64> = bounds
                .iter()
                .map(|&(lo, hi)| rng.uniform_in(lo, hi))
                .collect();
            let objectives = eval(&x);
            evaluations += 1;
            for (z, &o) in ideal.iter_mut().zip(&objectives) {
                if o.is_finite() {
                    *z = z.min(o);
                }
            }
            pop.push(MultiSolution { x, objectives });
        }
        if pop.is_empty() {
            return ParetoFront {
                solutions: Vec::new(),
                evaluations,
            };
        }
        let full_pop = pop.len() == n;

        while full_pop && evaluations < term.max_evaluations {
            #[allow(clippy::needless_range_loop)] // i indexes weights AND neighborhoods
            for i in 0..n {
                if evaluations >= term.max_evaluations {
                    break;
                }
                // Mating pool: the neighborhood with probability delta, the
                // whole population otherwise (MOEA/D-DE refinement).
                let use_nb = rng.uniform() < self.delta;
                let pick = |rng: &mut Rng| -> usize {
                    if use_nb {
                        neighborhoods[i][rng.index(t)]
                    } else {
                        rng.index(n)
                    }
                };
                let a = pick(&mut rng);
                let b = pick(&mut rng);

                let (mut child, _) = sbx(
                    &pop[a].x,
                    &pop[b].x,
                    bounds,
                    self.crossover_eta,
                    self.crossover_prob,
                    &mut rng,
                );
                mutate(&mut child, bounds, self.mutation_eta, pm, &mut rng);
                let objectives = eval(&child);
                evaluations += 1;
                for (z, &o) in ideal.iter_mut().zip(&objectives) {
                    if o.is_finite() {
                        *z = z.min(o);
                    }
                }
                let child = MultiSolution {
                    x: child,
                    objectives,
                };

                // Replace at most `nr` mating-pool members the child improves
                // (by their own Tchebycheff value), scanning in random order.
                let pool: &[usize] = if use_nb {
                    &neighborhoods[i]
                } else {
                    // Whole-population replacement scans a random sample of
                    // the same size as a neighborhood (keeps cost flat).
                    &[]
                };
                let mut replaced = 0usize;
                if pool.is_empty() {
                    // Sample T random subproblems for replacement.
                    for _ in 0..t {
                        if replaced >= self.max_replacements {
                            break;
                        }
                        let j = rng.index(n);
                        replaced += try_replace(&mut pop[j], &child, &weights[j], &ideal) as usize;
                    }
                } else {
                    // Random scan order over the neighborhood.
                    let mut order: Vec<usize> = pool.to_vec();
                    shuffle(&mut order, &mut rng);
                    for &j in &order {
                        if replaced >= self.max_replacements {
                            break;
                        }
                        replaced += try_replace(&mut pop[j], &child, &weights[j], &ideal) as usize;
                    }
                }
            }
        }

        let fronts = fast_non_dominated_sort(&pop);
        let solutions = fronts[0].iter().map(|&i| pop[i].clone()).collect();
        ParetoFront {
            solutions,
            evaluations,
        }
    }
}

/// Replaces `slot` with `child` when the child has a strictly better
/// Tchebycheff value for the slot's weight vector. Returns whether it did.
fn try_replace(
    slot: &mut MultiSolution,
    child: &MultiSolution,
    weight: &[f64],
    ideal: &[f64],
) -> bool {
    if tchebycheff(&child.objectives, weight, ideal) < tchebycheff(&slot.objectives, weight, ideal)
    {
        *slot = child.clone();
        true
    } else {
        false
    }
}

/// Tchebycheff scalarization `maxⱼ wⱼ · |fⱼ − z*ⱼ|` (Zhang & Li 2007, eq. 5).
fn tchebycheff(objectives: &[f64], weight: &[f64], ideal: &[f64]) -> f64 {
    objectives
        .iter()
        .zip(weight)
        .zip(ideal)
        .map(|((&f, &w), &z)| w * (f - z).abs())
        .fold(f64::NEG_INFINITY, f64::max)
}

fn dist2(a: &[f64], b: &[f64]) -> f64 {
    a.iter().zip(b).map(|(x, y)| (x - y) * (x - y)).sum()
}

/// Fisher–Yates shuffle with the crate RNG.
fn shuffle(v: &mut [usize], rng: &mut Rng) {
    for i in (1..v.len()).rev() {
        let j = rng.index(i + 1);
        v.swap(i, j);
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn tchebycheff_is_weighted_max_distance() {
        let g = tchebycheff(&[2.0, 3.0], &[0.5, 0.5], &[1.0, 1.0]);
        assert!((g - 1.0).abs() < 1e-12); // max(0.5·1, 0.5·2)
    }
}