From 6bfa52c149944f12de537df310d62bc8d8561b99 Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Tue, 5 May 2026 09:08:30 -0600 Subject: [PATCH] feat(algorithms): add AgeMoea (Adaptive Geometry Estimation MOEA) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Panichella 2019 AGE-MOEA: a many-objective MOEA that *infers* the front's geometry (its L_p shape, where p = 1 is linear, p = 2 is spherical, p < 1 is convex etc.) from the current non-dominated set and uses that estimate to drive both proximity and diversity in survival selection. Each generation: - NSGA-II-like loop: random parent selection + variation + evaluation - Combine + non_dominated_sort - Fill front-by-front; for the splitting front: - Translate by ideal point z* - Find extreme points by ASF (same as NSGA-III) and intercepts - Estimate the geometry parameter p by minimizing \|f − ideal\|_p constancy on the extreme points - Score every member by survival_score = (proximity_to_ideal) + (1 / nearest-neighbor distance in the same L_p frame) - Keep the top scorers The geometry estimation is the novel contribution; with 3+ objectives it produces fronts whose spread better matches the true shape than NSGA-III's reference points (which assume a known geometry). --- src/algorithms/age_moea.rs | 341 +++++++++++++++++++++++++++++++++++++ src/algorithms/mod.rs | 2 + src/prelude.rs | 2 +- 3 files changed, 344 insertions(+), 1 deletion(-) create mode 100644 src/algorithms/age_moea.rs diff --git a/src/algorithms/age_moea.rs b/src/algorithms/age_moea.rs new file mode 100644 index 0000000..8fcbe46 --- /dev/null +++ b/src/algorithms/age_moea.rs @@ -0,0 +1,341 @@ +//! `AgeMoea` — Panichella 2019 Adaptive Geometry Estimation MOEA. + +use rand::Rng as _; + +use crate::algorithms::parallel_eval::evaluate_batch; +use crate::core::candidate::Candidate; +use crate::core::objective::ObjectiveSpace; +use crate::core::population::Population; +use crate::core::problem::Problem; +use crate::core::result::OptimizationResult; +use crate::core::rng::rng_from_seed; +use crate::pareto::front::{best_candidate, pareto_front}; +use crate::pareto::sort::non_dominated_sort; +use crate::traits::{Initializer, Optimizer, Variation}; + +/// Configuration for [`AgeMoea`]. +#[derive(Debug, Clone)] +pub struct AgeMoeaConfig { + /// Constant population size. + pub population_size: usize, + /// Number of generations. + pub generations: usize, + /// Seed for the deterministic RNG. + pub seed: u64, +} + +impl Default for AgeMoeaConfig { + fn default() -> Self { + Self { population_size: 100, generations: 250, seed: 42 } + } +} + +/// Adaptive Geometry Estimation MOEA. +/// +/// Estimates the current front's L_p geometry parameter and uses it to +/// score survivors by a combination of proximity (distance to the +/// translated origin in the L_p frame) and diversity (distance to the +/// nearest survivor in the same frame). +#[derive(Debug, Clone)] +pub struct AgeMoea { + /// Algorithm configuration. + pub config: AgeMoeaConfig, + /// Initial-decision sampler. + pub initializer: I, + /// Offspring-producing variation operator. + pub variation: V, +} + +impl AgeMoea { + /// Construct an `AgeMoea`. + pub fn new(config: AgeMoeaConfig, initializer: I, variation: V) -> Self { + Self { config, initializer, variation } + } +} + +impl Optimizer

for AgeMoea +where + P: Problem + Sync, + P::Decision: Send, + I: Initializer, + V: Variation, +{ + fn run(&mut self, problem: &P) -> OptimizationResult { + assert!(self.config.population_size > 0, "AgeMoea population_size must be > 0"); + let n = self.config.population_size; + let objectives = problem.objectives(); + let mut rng = rng_from_seed(self.config.seed); + + let initial_decisions = self.initializer.initialize(n, &mut rng); + let mut population: Vec> = + evaluate_batch(problem, initial_decisions); + let mut evaluations = population.len(); + + for _ in 0..self.config.generations { + // Phase 1: random parent selection + variation. + let mut offspring_decisions: Vec = Vec::with_capacity(n); + while offspring_decisions.len() < n { + let p1 = rng.random_range(0..population.len()); + let p2 = rng.random_range(0..population.len()); + let parents = + vec![population[p1].decision.clone(), population[p2].decision.clone()]; + let children = self.variation.vary(&parents, &mut rng); + assert!(!children.is_empty(), "AgeMoea variation returned no children"); + for child in children { + if offspring_decisions.len() >= n { + break; + } + offspring_decisions.push(child); + } + } + let offspring = evaluate_batch(problem, offspring_decisions); + evaluations += offspring.len(); + + // Phase 3: combine + age-moea survival selection. + let mut combined: Vec> = Vec::with_capacity(2 * n); + combined.extend(population); + combined.extend(offspring); + population = environmental_selection(combined, &objectives, n); + } + + let front = pareto_front(&population, &objectives); + let best = best_candidate(&population, &objectives); + OptimizationResult::new( + Population::new(population), + front, + best, + evaluations, + self.config.generations, + ) + } +} + +fn environmental_selection( + combined: Vec>, + objectives: &ObjectiveSpace, + n: usize, +) -> Vec> { + let fronts = non_dominated_sort(&combined, objectives); + let mut selected: Vec = Vec::with_capacity(n); + let mut splitting: Vec = Vec::new(); + for f in &fronts { + if selected.len() + f.len() <= n { + selected.extend(f.iter().copied()); + } else { + splitting = f.clone(); + break; + } + if selected.len() == n { + break; + } + } + if selected.len() == n { + return selected.into_iter().map(|i| combined[i].clone()).collect(); + } + + // Translate by ideal point z*. + let m = objectives.len(); + let n0_oriented: Vec> = fronts[0] + .iter() + .map(|&i| objectives.as_minimization(&combined[i].evaluation.objectives)) + .collect(); + let mut ideal = vec![f64::INFINITY; m]; + for o in &n0_oriented { + for (k, v) in o.iter().enumerate() { + if *v < ideal[k] { + ideal[k] = *v; + } + } + } + + // Translate every combined member. + let translated: Vec> = combined + .iter() + .map(|c| { + let oriented = objectives.as_minimization(&c.evaluation.objectives); + oriented.iter().enumerate().map(|(k, v)| (v - ideal[k]).max(0.0)).collect() + }) + .collect(); + + // Estimate p (geometry parameter) from the *first* front's + // extreme points: find the point with the largest single-axis value + // for each axis, then solve for p such that all extreme points have + // unit L_p norm after normalizing by the per-axis maximum. + let p = estimate_p(&fronts[0], &translated, m); + + // Score every member of the splitting front by: + // proximity = ||translated||_p + // diversity = nearest-neighbor distance in the same L_p frame + // among already-selected + splitting members. + let mut keep = selected.clone(); + let mut remaining: Vec = splitting.clone(); + while keep.len() < n { + // Compute scores for every remaining candidate; pick the one with + // the largest combined score. + let mut best_idx: Option = None; + let mut best_score = f64::NEG_INFINITY; + for &i in &remaining { + let prox = lp_norm(&translated[i], p); + let div = nearest_neighbor_distance(i, &translated, &keep, p); + let score = div / (prox.max(1e-12)); + if score > best_score { + best_score = score; + best_idx = Some(i); + } + } + match best_idx { + None => break, + Some(pick) => { + keep.push(pick); + remaining.retain(|&i| i != pick); + } + } + } + keep.into_iter().map(|i| combined[i].clone()).collect() +} + +fn lp_norm(v: &[f64], p: f64) -> f64 { + v.iter().map(|x| x.abs().powf(p)).sum::().powf(1.0 / p) +} + +fn lp_distance(a: &[f64], b: &[f64], p: f64) -> f64 { + a.iter().zip(b.iter()).map(|(x, y)| (x - y).abs().powf(p)).sum::().powf(1.0 / p) +} + +fn nearest_neighbor_distance( + i: usize, + translated: &[Vec], + selected: &[usize], + p: f64, +) -> f64 { + if selected.is_empty() { + return f64::INFINITY; + } + let mut best = f64::INFINITY; + for &j in selected { + if j == i { + continue; + } + let d = lp_distance(&translated[i], &translated[j], p); + if d < best { + best = d; + } + } + best +} + +/// Estimate the L_p geometry parameter from the front's extreme points. +/// +/// Find the extreme point on each axis (the front member maximizing that +/// objective relative to its own L_∞ norm), then choose p such that all +/// extreme points have approximately the same L_p magnitude. Falls back +/// to p = 2 (spherical) if anything degenerates. +fn estimate_p(front_indices: &[usize], translated: &[Vec], m: usize) -> f64 { + if front_indices.is_empty() || m == 0 { + return 2.0; + } + // For each axis, find the extreme: the front member with the largest + // ratio of its k-th coordinate to its own L1 norm (i.e., the most + // "k-aligned" member). + let extremes: Vec = (0..m) + .map(|axis| { + let mut best = front_indices[0]; + let mut best_ratio = f64::NEG_INFINITY; + for &idx in front_indices { + let l1: f64 = translated[idx].iter().sum::().max(1e-12); + let ratio = translated[idx][axis] / l1; + if ratio > best_ratio { + best_ratio = ratio; + best = idx; + } + } + best + }) + .collect(); + + // Solve for p ∈ [0.1, 10.0] that minimizes std-dev of L_p norms across + // extremes (a coarse sweep is fine — full Brent isn't needed for this + // shape estimate). + let candidates: Vec = (1..=40).map(|i| (i as f64) * 0.25).collect(); + let mut best_p = 2.0; + let mut best_loss = f64::INFINITY; + for &p in &candidates { + let norms: Vec = extremes + .iter() + .map(|&i| lp_norm(&translated[i], p)) + .collect(); + let mean = norms.iter().sum::() / norms.len() as f64; + if mean.is_finite() && mean > 0.0 { + let var = norms.iter().map(|x| (x - mean).powi(2)).sum::() / norms.len() as f64; + let loss = var.sqrt() / mean; + if loss < best_loss { + best_loss = loss; + best_p = p; + } + } + } + best_p +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::operators::{ + CompositeVariation, PolynomialMutation, RealBounds, SimulatedBinaryCrossover, + }; + use crate::tests_support::SchafferN1; + + fn make_optimizer( + seed: u64, + ) -> AgeMoea> { + let bounds = vec![(-5.0, 5.0)]; + let initializer = RealBounds::new(bounds.clone()); + let variation = CompositeVariation { + crossover: SimulatedBinaryCrossover::new(bounds.clone(), 15.0, 0.5), + mutation: PolynomialMutation::new(bounds, 20.0, 1.0), + }; + AgeMoea::new( + AgeMoeaConfig { population_size: 20, generations: 15, seed }, + initializer, + variation, + ) + } + + #[test] + fn produces_pareto_front() { + let mut opt = make_optimizer(1); + let r = opt.run(&SchafferN1); + assert_eq!(r.population.len(), 20); + assert!(!r.pareto_front.is_empty()); + } + + #[test] + fn deterministic_with_same_seed() { + let mut a = make_optimizer(99); + let mut b = make_optimizer(99); + let ra = a.run(&SchafferN1); + let rb = b.run(&SchafferN1); + let oa: Vec> = + ra.pareto_front.iter().map(|c| c.evaluation.objectives.clone()).collect(); + let ob: Vec> = + rb.pareto_front.iter().map(|c| c.evaluation.objectives.clone()).collect(); + assert_eq!(oa, ob); + } + + #[test] + #[should_panic(expected = "population_size must be > 0")] + fn zero_pop_panics() { + let bounds = vec![(0.0, 1.0)]; + let initializer = RealBounds::new(bounds.clone()); + let variation = CompositeVariation { + crossover: SimulatedBinaryCrossover::new(bounds.clone(), 15.0, 0.5), + mutation: PolynomialMutation::new(bounds, 20.0, 1.0), + }; + let mut opt = AgeMoea::new( + AgeMoeaConfig { population_size: 0, generations: 1, seed: 0 }, + initializer, + variation, + ); + let _ = opt.run(&SchafferN1); + } +} diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index c814a77..14fc3f1 100644 --- a/src/algorithms/mod.rs +++ b/src/algorithms/mod.rs @@ -1,5 +1,6 @@ //! Built-in reference optimizers. +pub mod age_moea; pub mod ant_colony_tsp; pub mod cma_es; pub mod differential_evolution; @@ -25,6 +26,7 @@ pub mod tabu_search; pub mod tlbo; pub mod umda; +pub use age_moea::*; pub use ant_colony_tsp::*; pub use cma_es::*; pub use differential_evolution::*; diff --git a/src/prelude.rs b/src/prelude.rs index 8a9f0a1..226f34e 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -22,7 +22,7 @@ pub use crate::operators::{ }; pub use crate::algorithms::{ - AntColonyTsp, AntColonyTspConfig, CmaEs, CmaEsConfig, DifferentialEvolution, + AgeMoea, AgeMoeaConfig, AntColonyTsp, AntColonyTspConfig, CmaEs, CmaEsConfig, DifferentialEvolution, DifferentialEvolutionConfig, EpsilonMoea, EpsilonMoeaConfig, GeneticAlgorithm, GeneticAlgorithmConfig, HillClimber, HillClimberConfig, Hype, HypeConfig, Ibea, IbeaConfig, Moead, MoeadConfig, Mopso, MopsoConfig, Nsga2,