From 283d7429bba189a5e1881a95a8175f33e93fb484 Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Tue, 5 May 2026 08:23:45 -0600 Subject: [PATCH] feat(algorithms): add Rvea (Reference Vector-guided EA) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Cheng, Jin, Olhofer & Sendhoff 2016 RVEA: many-objective MOEA built around a fixed set of Das–Dennis reference vectors. Each generation: - Generate offspring via random parent selection + variation + evaluation - Combine population + offspring; translate by ideal point z* - Associate every member with the reference vector whose angle to the translated objective vector is smallest - For each occupied vector, keep the member with the smallest Angle-Penalized Distance (APD) score; the rest are dropped - APD = (1 + α(t)·θ_max·γ) · |f − z*| where γ is the angle to the associated reference and α(t) = (t / t_max)^2 anneals the angle penalty over the run This produces well-spread fronts at high objective counts where Pareto-rank methods (NSGA-II, SPEA2) lose discrimination. --- src/algorithms/mod.rs | 2 + src/algorithms/rvea.rs | 328 +++++++++++++++++++++++++++++++++++++++++ src/prelude.rs | 3 +- 3 files changed, 332 insertions(+), 1 deletion(-) create mode 100644 src/algorithms/rvea.rs diff --git a/src/algorithms/mod.rs b/src/algorithms/mod.rs index ff3803a..b2840c1 100644 --- a/src/algorithms/mod.rs +++ b/src/algorithms/mod.rs @@ -15,6 +15,7 @@ pub mod paes; pub(crate) mod parallel_eval; pub mod particle_swarm; pub mod random_search; +pub mod rvea; pub mod simulated_annealing; pub mod sms_emoa; pub mod spea2; @@ -35,6 +36,7 @@ pub use nsga3::*; pub use paes::*; pub use particle_swarm::*; pub use random_search::*; +pub use rvea::*; pub use simulated_annealing::*; pub use sms_emoa::*; pub use spea2::*; diff --git a/src/algorithms/rvea.rs b/src/algorithms/rvea.rs new file mode 100644 index 0000000..7a1b821 --- /dev/null +++ b/src/algorithms/rvea.rs @@ -0,0 +1,328 @@ +//! `Rvea` — Cheng, Jin, Olhofer & Sendhoff 2016 Reference Vector-guided EA. + +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::reference_points::das_dennis; +use crate::traits::{Initializer, Optimizer, Variation}; + +/// Configuration for [`Rvea`]. +#[derive(Debug, Clone)] +pub struct RveaConfig { + /// Constant population size. + pub population_size: usize, + /// Number of generations. + pub generations: usize, + /// Number of divisions `H` for Das–Dennis reference vectors. Pop size + /// should be roughly `binomial(H + M − 1, M − 1)`. + pub reference_divisions: usize, + /// Penalty exponent `α`. The paper recommends 2.0. + pub alpha: f64, + /// Seed for the deterministic RNG. + pub seed: u64, +} + +impl Default for RveaConfig { + fn default() -> Self { + Self { + population_size: 100, + generations: 250, + reference_divisions: 12, + alpha: 2.0, + seed: 42, + } + } +} + +/// Reference Vector-guided Evolutionary Algorithm. +#[derive(Debug, Clone)] +pub struct Rvea { + /// Algorithm configuration. + pub config: RveaConfig, + /// Initial-decision sampler. + pub initializer: I, + /// Offspring-producing variation operator. + pub variation: V, +} + +impl Rvea { + /// Construct an `Rvea`. + pub fn new(config: RveaConfig, initializer: I, variation: V) -> Self { + Self { config, initializer, variation } + } +} + +impl Optimizer

for Rvea +where + P: Problem + Sync, + P::Decision: Send, + I: Initializer, + V: Variation, +{ + fn run(&mut self, problem: &P) -> OptimizationResult { + assert!(self.config.population_size > 0, "Rvea population_size must be > 0"); + let n = self.config.population_size; + let objectives = problem.objectives(); + let m = objectives.len(); + // Reference vectors normalized to unit norm. + let raw_refs = das_dennis(m, self.config.reference_divisions); + let references: Vec> = raw_refs.into_iter().map(unit_normalize).collect(); + assert!(!references.is_empty(), "Rvea: no reference vectors generated"); + + // Smallest angle between any two reference vectors — used to scale + // the APD penalty term. + let theta_max = smallest_neighbor_angle(&references); + 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 gen_idx 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(), "Rvea 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(); + + // Combine + APD-based survival. + let mut combined: Vec> = Vec::with_capacity(2 * n); + combined.extend(population); + combined.extend(offspring); + + // Ideal point z*. + let m_dim = m; + let mut ideal = vec![f64::INFINITY; m_dim]; + for c in &combined { + let oriented = objectives.as_minimization(&c.evaluation.objectives); + for (k, v) in oriented.iter().enumerate() { + if *v < ideal[k] { + ideal[k] = *v; + } + } + } + // Translate. + let translated: Vec> = combined + .iter() + .map(|c| { + let oriented = objectives.as_minimization(&c.evaluation.objectives); + oriented.iter().enumerate().map(|(k, v)| v - ideal[k]).collect() + }) + .collect(); + + // Associate each member with its closest-angle reference vector. + let mut assoc: Vec = vec![0; combined.len()]; + let mut angles: Vec = vec![0.0; combined.len()]; + for (i, t) in translated.iter().enumerate() { + let (best_ref, best_angle) = closest_reference(t, &references); + assoc[i] = best_ref; + angles[i] = best_angle; + } + + // For each occupied reference vector, keep the member with the + // smallest APD score. + let alpha_t = (gen_idx as f64 / (self.config.generations as f64).max(1.0)) + .powf(self.config.alpha); + let mut keep: Vec> = vec![None; references.len()]; + for i in 0..combined.len() { + let r = assoc[i]; + let length: f64 = translated[i].iter().map(|v| v * v).sum::().sqrt(); + let theta_max_safe = theta_max.max(1e-12); + let penalty = 1.0 + (m_dim as f64) * alpha_t * (angles[i] / theta_max_safe); + let apd = penalty * length; + match keep[r] { + None => keep[r] = Some((i, apd)), + Some((_, current)) if apd < current => keep[r] = Some((i, apd)), + _ => {} + } + } + + let mut next: Vec> = + keep.into_iter().flatten().map(|(i, _)| combined[i].clone()).collect(); + // If we ended up with fewer than n (some references unfilled), + // backfill with the lowest-APD remaining candidates. + if next.len() < n { + let mut all_apds: Vec<(usize, f64)> = (0..combined.len()) + .map(|i| { + let length: f64 = + translated[i].iter().map(|v| v * v).sum::().sqrt(); + let theta_max_safe = theta_max.max(1e-12); + let penalty = 1.0 + (m_dim as f64) * alpha_t * (angles[i] / theta_max_safe); + (i, penalty * length) + }) + .collect(); + all_apds + .sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal)); + for (i, _) in all_apds { + if next.len() >= n { + break; + } + if !next + .iter() + .any(|c| std::ptr::eq(c as *const _, &combined[i] as *const _)) + { + next.push(combined[i].clone()); + } + } + } + // If too many (only possible if the reference set has > n + // vectors), truncate by APD. + if next.len() > n { + next.truncate(n); + } + population = next; + } + + let front = pareto_front(&population, &objectives); + let best = best_candidate(&population, &objectives); + OptimizationResult::new( + Population::new(population), + front, + best, + evaluations, + self.config.generations, + ) + } +} + +fn unit_normalize(mut v: Vec) -> Vec { + let n: f64 = v.iter().map(|x| x * x).sum::().sqrt(); + if n > 1e-12 { + for x in v.iter_mut() { + *x /= n; + } + } + v +} + +fn closest_reference(point: &[f64], references: &[Vec]) -> (usize, f64) { + let length: f64 = point.iter().map(|v| v * v).sum::().sqrt().max(1e-12); + let mut best = 0; + let mut best_angle = f64::INFINITY; + for (i, r) in references.iter().enumerate() { + let dot: f64 = point.iter().zip(r.iter()).map(|(a, b)| a * b).sum(); + let cosine = (dot / length).clamp(-1.0, 1.0); + let angle = cosine.acos(); + if angle < best_angle { + best_angle = angle; + best = i; + } + } + (best, best_angle) +} + +fn smallest_neighbor_angle(references: &[Vec]) -> f64 { + let mut min_angle = f64::INFINITY; + for i in 0..references.len() { + for j in (i + 1)..references.len() { + let dot: f64 = references[i] + .iter() + .zip(references[j].iter()) + .map(|(a, b)| a * b) + .sum(); + let angle = dot.clamp(-1.0, 1.0).acos(); + if angle < min_angle { + min_angle = angle; + } + } + } + if !min_angle.is_finite() { std::f64::consts::FRAC_PI_4 } else { min_angle } +} + +#[allow(unused_imports)] +use crate::core::objective::Objective; + +#[cfg(test)] +mod tests { + use super::*; + use crate::operators::{ + CompositeVariation, PolynomialMutation, RealBounds, SimulatedBinaryCrossover, + }; + use crate::tests_support::SchafferN1; + + fn make_optimizer( + seed: u64, + ) -> Rvea> { + 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), + }; + Rvea::new( + RveaConfig { + population_size: 20, + generations: 15, + reference_divisions: 19, + alpha: 2.0, + seed, + }, + initializer, + variation, + ) + } + + #[test] + fn produces_pareto_front() { + let mut opt = make_optimizer(1); + let r = opt.run(&SchafferN1); + 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_population_size_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 = Rvea::new( + RveaConfig { + population_size: 0, + generations: 1, + reference_divisions: 5, + alpha: 2.0, + seed: 0, + }, + initializer, + variation, + ); + let _ = opt.run(&SchafferN1); + } +} diff --git a/src/prelude.rs b/src/prelude.rs index 656cefe..4720687 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -26,7 +26,8 @@ pub use crate::algorithms::{ DifferentialEvolutionConfig, GeneticAlgorithm, GeneticAlgorithmConfig, HillClimber, HillClimberConfig, Hype, HypeConfig, Ibea, IbeaConfig, Moead, MoeadConfig, Mopso, MopsoConfig, Nsga2, Nsga2Config, Nsga3, Nsga3Config, Paes, PaesConfig, ParticleSwarm, - ParticleSwarmConfig, RandomSearch, RandomSearchConfig, SimulatedAnnealing, + ParticleSwarmConfig, RandomSearch, RandomSearchConfig, Rvea, RveaConfig, + SimulatedAnnealing, SimulatedAnnealingConfig, SmsEmoa, SmsEmoaConfig, Spea2, Spea2Config, TabuSearch, TabuSearchConfig, Umda, UmdaConfig,