//! Shared workload for the `compare` example and the `compare_profile` //! gungraun benchmark. //! //! Holds every problem definition, the ~88 algorithm runner functions, and //! the table-printing presentation layer. The example //! (`examples/compare.rs`) is a thin shim that calls [`run_all`]; the //! profiling benchmark (`benches/compare_profile.rs`) reuses the runner //! functions directly so it exercises the exact same workload. //! //! This file is `#[path]`-included as a module by both, and lives under //! `examples/_shared/` (no `main.rs`) so cargo does not auto-discover it //! as an example. use std::f64::consts::PI; use std::time::Instant; use rand::Rng as _; use heuropt::core::rng::Rng; use heuropt::metrics::{hypervolume::hypervolume_2d, spacing::spacing}; use heuropt::prelude::*; pub const SEEDS: u64 = 10; pub const ZDT1_DIM: usize = 30; pub const ZDT1_BUDGET: usize = 25_000; // Standard ZDT1 reference point. Using [11, 11] (rather than the // near-front [1.1, 1.1]) so under-converged algorithms with large `g` // values still register a meaningful — if poor — hypervolume. pub const ZDT1_REFERENCE: [f64; 2] = [11.0, 11.0]; pub const RASTRIGIN_DIM: usize = 5; pub const RASTRIGIN_BUDGET: usize = 50_000; pub const DTLZ2_OBJECTIVES: usize = 3; pub const DTLZ2_K: usize = 10; pub const DTLZ2_DIM: usize = DTLZ2_OBJECTIVES + DTLZ2_K - 1; // 12 pub const DTLZ2_BUDGET: usize = 30_000; pub const ROSENBROCK_DIM: usize = 5; pub const ROSENBROCK_BUDGET: usize = 30_000; pub const ACKLEY_DIM: usize = 5; pub const ACKLEY_BUDGET: usize = 30_000; pub const ZDT3_DIM: usize = 30; pub const ZDT3_BUDGET: usize = 25_000; pub const ZDT3_REFERENCE: [f64; 2] = [11.0, 11.0]; pub const DTLZ1_OBJECTIVES: usize = 3; pub const DTLZ1_K: usize = 5; pub const DTLZ1_DIM: usize = DTLZ1_OBJECTIVES + DTLZ1_K - 1; pub const DTLZ1_BUDGET: usize = 30_000; // ----------------------------------------------------------------------------- // Test problems // ----------------------------------------------------------------------------- pub struct Zdt1 { dim: usize, } impl Problem for Zdt1 { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![Objective::minimize("f1"), Objective::minimize("f2")]) } fn evaluate(&self, x: &Vec) -> Evaluation { let f1 = x[0]; let tail_sum: f64 = x[1..].iter().sum(); let g = 1.0 + 9.0 * tail_sum / (self.dim as f64 - 1.0); let f2 = g * (1.0 - (f1 / g).sqrt()); Evaluation::new(vec![f1, f2]) } } pub struct Dtlz2 { num_objectives: usize, dim: usize, } impl Problem for Dtlz2 { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new( (0..self.num_objectives) .map(|i| Objective::minimize(format!("f{}", i + 1))) .collect(), ) } fn evaluate(&self, x: &Vec) -> Evaluation { let m = self.num_objectives; let g: f64 = x[(m - 1)..self.dim].iter().map(|v| (v - 0.5).powi(2)).sum(); let scale = 1.0 + g; let mut f = vec![0.0_f64; m]; for i in 0..m { let mut prod = scale; #[allow(clippy::needless_range_loop)] // Body indexes `x[j]`. for j in 0..(m - i - 1) { prod *= (x[j] * std::f64::consts::FRAC_PI_2).cos(); } if i > 0 { prod *= (x[m - i - 1] * std::f64::consts::FRAC_PI_2).sin(); } f[i] = prod; } Evaluation::new(f) } } pub struct Rosenbrock { dim: usize, } impl Problem for Rosenbrock { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![Objective::minimize("f")]) } fn evaluate(&self, x: &Vec) -> Evaluation { let f: f64 = (0..self.dim - 1) .map(|i| { let a = 1.0 - x[i]; let b = x[i + 1] - x[i] * x[i]; a * a + 100.0 * b * b }) .sum(); Evaluation::new(vec![f]) } } pub struct Ackley { dim: usize, } impl Problem for Ackley { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![Objective::minimize("f")]) } fn evaluate(&self, x: &Vec) -> Evaluation { let n = self.dim as f64; let sum_sq: f64 = x.iter().map(|v| v * v).sum(); let sum_cos: f64 = x.iter().map(|v| (2.0 * PI * v).cos()).sum(); let f = -20.0 * (-0.2 * (sum_sq / n).sqrt()).exp() - (sum_cos / n).exp() + 20.0 + std::f64::consts::E; Evaluation::new(vec![f]) } } pub struct Zdt3 { dim: usize, } impl Problem for Zdt3 { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![Objective::minimize("f1"), Objective::minimize("f2")]) } fn evaluate(&self, x: &Vec) -> Evaluation { let f1 = x[0]; let tail_sum: f64 = x[1..].iter().sum(); let g = 1.0 + 9.0 * tail_sum / (self.dim as f64 - 1.0); let r = f1 / g; let f2 = g * (1.0 - r.sqrt() - r * (10.0 * PI * f1).sin()); Evaluation::new(vec![f1, f2]) } } pub struct Dtlz1 { num_objectives: usize, dim: usize, } impl Problem for Dtlz1 { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new( (0..self.num_objectives) .map(|i| Objective::minimize(format!("f{}", i + 1))) .collect(), ) } fn evaluate(&self, x: &Vec) -> Evaluation { let m = self.num_objectives; let k = self.dim - (m - 1); let g_term: f64 = x[(m - 1)..self.dim] .iter() .map(|v| (v - 0.5).powi(2) - (20.0 * PI * (v - 0.5)).cos()) .sum(); let g = 100.0 * (k as f64 + g_term); let mut f = vec![0.0_f64; m]; for i in 0..m { let mut prod = 0.5 * (1.0 + g); #[allow(clippy::needless_range_loop)] // body indexes x[j]. for j in 0..(m - i - 1) { prod *= x[j]; } if i > 0 { prod *= 1.0 - x[m - i - 1]; } f[i] = prod; } Evaluation::new(f) } } pub struct Rastrigin { dim: usize, } impl Problem for Rastrigin { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![Objective::minimize("f")]) } fn evaluate(&self, x: &Vec) -> Evaluation { let n = self.dim as f64; let value = 10.0 * n + x.iter() .map(|v| v * v - 10.0 * (2.0 * PI * v).cos()) .sum::(); Evaluation::new(vec![value]) } } // ----------------------------------------------------------------------------- // Run results + metrics aggregation // ----------------------------------------------------------------------------- #[derive(Clone)] pub struct MoRun { front: Vec>>, wall_ms: u128, } #[derive(Clone)] pub struct SoRun { best_value: f64, wall_ms: u128, } pub fn mean_l2_to_zdt1_front(front: &[Candidate>]) -> f64 { if front.is_empty() { return f64::INFINITY; } let samples: Vec<(f64, f64)> = (0..=1000) .map(|i| { let f1 = i as f64 / 1000.0; (f1, 1.0 - f1.sqrt()) }) .collect(); let mut total = 0.0; for c in front { let f1 = c.evaluation.objectives[0]; let f2 = c.evaluation.objectives[1]; let mut best = f64::INFINITY; for &(rf1, rf2) in &samples { let d = ((rf1 - f1).powi(2) + (rf2 - f2).powi(2)).sqrt(); if d < best { best = d; } } total += best; } total / front.len() as f64 } pub fn mean_std(values: &[f64]) -> (f64, f64) { let n = values.len() as f64; let mean = values.iter().sum::() / n; let var = values.iter().map(|v| (v - mean).powi(2)).sum::() / n; (mean, var.sqrt()) } /// Print an aligned text table: column 0 left-justified, the rest /// right-justified. Column widths are derived from the actual cell /// contents (header *and* every row), so the separator and all rows line /// up no matter how the value magnitudes vary. /// /// All cells must be ASCII — width is measured with `str::len`, so a /// multi-byte character (e.g. `±`, `ε`) would silently break alignment on /// terminals that render it at a different column width. Callers format /// `mean +/- std` rather than `mean ± std` for exactly this reason. pub fn print_table(header: &[&str], rows: &[Vec]) { let ncols = header.len(); let mut widths: Vec = header.iter().map(|h| h.len()).collect(); for row in rows { for c in 0..ncols { widths[c] = widths[c].max(row[c].len()); } } let fmt_row = |cells: &[String]| -> String { let mut out = String::new(); for c in 0..ncols { if c > 0 { out.push_str(" "); } let w = widths[c]; if c == 0 { out.push_str(&format!("{:w$}", cells[c])); } } out }; let header_owned: Vec = header.iter().map(|s| s.to_string()).collect(); let header_line = fmt_row(&header_owned); println!("{header_line}"); println!("{}", "-".repeat(header_line.len())); for row in rows { println!("{}", fmt_row(row)); } } // ----------------------------------------------------------------------------- // ZDT1 algorithm runners // ----------------------------------------------------------------------------- pub fn zdt1_random(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let initializer = RealBounds::new(vec![(0.0, 1.0); ZDT1_DIM]); let config = RandomSearchConfig { iterations: ZDT1_BUDGET, batch_size: 1, seed, }; let mut opt = RandomSearch::new(config, initializer); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_paes(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let initializer = RealBounds::new(vec![(0.0, 1.0); ZDT1_DIM]); let variation = BoundedGaussianMutation::new(0.05, vec![(0.0, 1.0); ZDT1_DIM]); let config = PaesConfig { iterations: ZDT1_BUDGET, archive_size: 100, seed, }; let mut opt = Paes::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_spea2(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; let pop = 100; let arc = 100; // SPEA2 evaluates `pop_size` per generation after the initial population. let gens = (ZDT1_BUDGET - pop) / pop; let config = Spea2Config { population_size: pop, archive_size: arc, generations: gens, seed, }; let mut opt = Spea2::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_nsga2(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; let pop = 100; let gens = ZDT1_BUDGET / pop; let config = Nsga2Config { population_size: pop, generations: gens, seed, }; let mut opt = Nsga2::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_sms_emoa(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; // SMS-EMOA is steady-state and computes O(N²) hypervolumes per // iteration, so we run it on a smaller population for a smaller // total budget to keep wall time tractable. let config = SmsEmoaConfig { population_size: 40, generations: 4_000, reference_point: vec![11.0, 11.0], seed, }; let mut opt = SmsEmoa::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_hype(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; let pop = 80; // Each generation does multiple HV-estimation passes (parent selection + // survival truncation), each with mc_samples × N work. Keep budget // smaller than NSGA-II's so wall time is tractable. let gens = 80; let config = HypeConfig { population_size: pop, generations: gens, reference_point: vec![11.0, 11.0], mc_samples: 1_000, seed, }; let mut opt = Hype::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_rvea(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; let pop = 100; let gens = ZDT1_BUDGET / pop; let config = RveaConfig { population_size: pop, generations: gens, reference_divisions: 99, alpha: 2.0, seed, }; let mut opt = Rvea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_pesa2(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; let pop = 50; let gens = (ZDT1_BUDGET - pop) / pop; let config = PesaIIConfig { population_size: pop, archive_size: 100, generations: gens, grid_divisions: 20, seed, }; let mut opt = PesaII::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_epsilon_moea(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; let config = EpsilonMoeaConfig { population_size: 50, evaluations: ZDT1_BUDGET, epsilon: vec![0.01, 0.01], seed, }; let mut opt = EpsilonMoea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_mopso(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = RealBounds::new(vec![(0.0, 1.0); ZDT1_DIM]); let swarm = 100; let gens = (ZDT1_BUDGET - 2 * swarm) / swarm; let config = MopsoConfig { swarm_size: swarm, generations: gens, archive_size: 100, inertia: 0.7, cognitive: 1.5, social: 1.5, seed, }; let mut opt = Mopso::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_ibea(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; let pop = 100; let gens = ZDT1_BUDGET / pop; let config = IbeaConfig { population_size: pop, generations: gens, kappa: 0.05, seed, }; let mut opt = Ibea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_moead(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; // 99 divisions → 100 weights for 2 obj. Each generation evaluates one // child per weight (so `n_weights` evals/gen). let pop = 100; let gens = (ZDT1_BUDGET - pop) / pop; let config = MoeadConfig { generations: gens, reference_divisions: 99, neighborhood_size: 20, seed, }; let mut opt = Moead::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt1_nsga3(seed: u64) -> MoRun { let problem = Zdt1 { dim: ZDT1_DIM }; let bounds = vec![(0.0, 1.0); ZDT1_DIM]; 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 / ZDT1_DIM as f64), }; let pop = 100; let gens = ZDT1_BUDGET / pop; let config = Nsga3Config { population_size: pop, generations: gens, // 99 ref points for 2 objectives — same density as the population. reference_divisions: 99, seed, }; let mut opt = Nsga3::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } // ----------------------------------------------------------------------------- // DTLZ2 algorithm runners (3-objective) // ----------------------------------------------------------------------------- pub fn dtlz2_problem() -> Dtlz2 { Dtlz2 { num_objectives: DTLZ2_OBJECTIVES, dim: DTLZ2_DIM, } } pub fn dtlz2_random(seed: u64) -> MoRun { let problem = dtlz2_problem(); let initializer = RealBounds::new(vec![(0.0, 1.0); DTLZ2_DIM]); let config = RandomSearchConfig { iterations: DTLZ2_BUDGET, batch_size: 1, seed, }; let mut opt = RandomSearch::new(config, initializer); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_nsga2(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; let pop = 92; // close to the 91-ref-point NSGA-III pop, for fairness let gens = DTLZ2_BUDGET / pop; let config = Nsga2Config { population_size: pop, generations: gens, seed, }; let mut opt = Nsga2::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_spea2(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; let pop = 92; let arc = 92; let gens = (DTLZ2_BUDGET - pop) / pop; let config = Spea2Config { population_size: pop, archive_size: arc, generations: gens, seed, }; let mut opt = Spea2::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_sms_emoa(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; // 3-D HV-by-slicing is significantly more expensive than 2-D; smaller // pop and gen budget here to keep the comparison tractable. let pop = 40; let gens = 4_000; let config = SmsEmoaConfig { population_size: pop, generations: gens, reference_point: vec![3.0, 3.0, 3.0], seed, }; let mut opt = SmsEmoa::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_hype(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; let pop = 80; let gens = 80; let config = HypeConfig { population_size: pop, generations: gens, reference_point: vec![3.0, 3.0, 3.0], mc_samples: 1_000, seed, }; let mut opt = Hype::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_rvea(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; let pop = 92; let gens = DTLZ2_BUDGET / pop; let config = RveaConfig { population_size: pop, generations: gens, reference_divisions: 12, alpha: 2.0, seed, }; let mut opt = Rvea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_pesa2(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; let pop = 50; let gens = (DTLZ2_BUDGET - pop) / pop; let config = PesaIIConfig { population_size: pop, archive_size: 100, generations: gens, grid_divisions: 12, seed, }; let mut opt = PesaII::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_epsilon_moea(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; let config = EpsilonMoeaConfig { population_size: 50, evaluations: DTLZ2_BUDGET, epsilon: vec![0.05, 0.05, 0.05], seed, }; let mut opt = EpsilonMoea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_mopso(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = RealBounds::new(vec![(0.0, 1.0); DTLZ2_DIM]); let swarm = 92; let gens = (DTLZ2_BUDGET - 2 * swarm) / swarm; let config = MopsoConfig { swarm_size: swarm, generations: gens, archive_size: 100, inertia: 0.7, cognitive: 1.5, social: 1.5, seed, }; let mut opt = Mopso::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_ibea(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; let pop = 92; let gens = DTLZ2_BUDGET / pop; let config = IbeaConfig { population_size: pop, generations: gens, kappa: 0.05, seed, }; let mut opt = Ibea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_moead(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; // 12 divisions for 3 objectives = 91 weights — same density as NSGA-III. let pop = 91; let gens = (DTLZ2_BUDGET - pop) / pop; let config = MoeadConfig { generations: gens, reference_divisions: 12, neighborhood_size: 20, seed, }; let mut opt = Moead::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz2_nsga3(seed: u64) -> MoRun { let problem = dtlz2_problem(); let bounds = vec![(0.0, 1.0); DTLZ2_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ2_DIM as f64), }; // H=12 → 91 reference points (the canonical NSGA-III 3-objective set). // Population is sized to match: the spec recommends pop ≈ #refs. let pop = 92; let gens = DTLZ2_BUDGET / pop; let config = Nsga3Config { population_size: pop, generations: gens, reference_divisions: 12, seed, }; let mut opt = Nsga3::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } /// DTLZ2's analytical Pareto front is the unit sphere octant in objective /// space (`Σ f_i² = 1`, all `f_i ≥ 0`). The closest-point distance from /// `f` to that surface is `|‖f‖ - 1|`. pub fn mean_distance_to_dtlz2_front(front: &[Candidate>]) -> f64 { if front.is_empty() { return f64::INFINITY; } let total: f64 = front .iter() .map(|c| { let norm: f64 = c .evaluation .objectives .iter() .map(|v| v * v) .sum::() .sqrt(); (norm - 1.0).abs() }) .sum(); total / front.len() as f64 } // ----------------------------------------------------------------------------- // Rastrigin algorithm runners // ----------------------------------------------------------------------------- pub fn rastrigin_random(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let initializer = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let config = RandomSearchConfig { iterations: RASTRIGIN_BUDGET, batch_size: 1, seed, }; let mut opt = RandomSearch::new(config, initializer); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_paes(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let initializer = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let variation = BoundedGaussianMutation::new(0.3, vec![(-5.12, 5.12); RASTRIGIN_DIM]); let config = PaesConfig { iterations: RASTRIGIN_BUDGET, archive_size: 32, seed, }; let mut opt = Paes::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_nsga2(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let bounds = vec![(-5.12, 5.12); RASTRIGIN_DIM]; 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 / RASTRIGIN_DIM as f64), }; let pop = 50; let gens = RASTRIGIN_BUDGET / pop; let config = Nsga2Config { population_size: pop, generations: gens, seed, }; let mut opt = Nsga2::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_de(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let bounds = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let pop = 50; let gens = (RASTRIGIN_BUDGET - pop) / pop; // initial pop also evaluates let config = DifferentialEvolutionConfig { population_size: pop, generations: gens, differential_weight: 0.5, crossover_probability: 0.9, seed, }; let mut opt = DifferentialEvolution::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_hill_climber(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let initializer = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let variation = BoundedGaussianMutation::new(0.3, vec![(-5.12, 5.12); RASTRIGIN_DIM]); let config = HillClimberConfig { iterations: RASTRIGIN_BUDGET, seed, }; let mut opt = HillClimber::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_simulated_annealing(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let initializer = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let variation = BoundedGaussianMutation::new(0.5, vec![(-5.12, 5.12); RASTRIGIN_DIM]); let config = SimulatedAnnealingConfig { iterations: RASTRIGIN_BUDGET, initial_temperature: 5.0, final_temperature: 1e-3, seed, }; let mut opt = SimulatedAnnealing::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_genetic_algorithm(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let bounds = vec![(-5.12, 5.12); RASTRIGIN_DIM]; 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 / RASTRIGIN_DIM as f64), }; let pop = 50; let gens = (RASTRIGIN_BUDGET - pop) / pop; let config = GeneticAlgorithmConfig { population_size: pop, generations: gens, tournament_size: 2, elitism: 2, seed, }; let mut opt = GeneticAlgorithm::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_particle_swarm(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let bounds = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let swarm = 40; // PSO does swarm + swarm·gens + final evaluations. Approximate budget. let gens = (RASTRIGIN_BUDGET - 2 * swarm) / swarm; let config = ParticleSwarmConfig { swarm_size: swarm, generations: gens, inertia: 0.7, cognitive: 1.5, social: 1.5, seed, }; let mut opt = ParticleSwarm::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_cma_es(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let bounds = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let pop = 16; let gens = RASTRIGIN_BUDGET / pop; let config = CmaEsConfig { population_size: pop, generations: gens, initial_sigma: 1.0, eigen_decomposition_period: 1, initial_mean: None, seed, }; let mut opt = CmaEs::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_ipop_cma_es(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let bounds = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let pop = 16; let config = IpopCmaEsConfig { initial_population_size: pop, total_generations: RASTRIGIN_BUDGET / pop, initial_sigma: 1.0, eigen_decomposition_period: 1, stall_generations: None, seed, }; let mut opt = IpopCmaEs::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rastrigin_one_plus_one_es(seed: u64) -> SoRun { let problem = Rastrigin { dim: RASTRIGIN_DIM }; let bounds = RealBounds::new(vec![(-5.12, 5.12); RASTRIGIN_DIM]); let config = OnePlusOneEsConfig { iterations: RASTRIGIN_BUDGET, initial_sigma: 1.0, adaptation_period: 50, step_increase: 1.22, seed, }; let mut opt = OnePlusOneEs::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rosenbrock_nelder_mead(_seed: u64) -> SoRun { let problem = rosenbrock_problem(); let bounds = RealBounds::new(vec![(-5.0, 10.0); ROSENBROCK_DIM]); let config = NelderMeadConfig { iterations: ROSENBROCK_BUDGET / 4, ..NelderMeadConfig::default() }; let mut opt = NelderMead::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rosenbrock_one_plus_one_es(seed: u64) -> SoRun { let problem = rosenbrock_problem(); let bounds = RealBounds::new(vec![(-5.0, 10.0); ROSENBROCK_DIM]); let config = OnePlusOneEsConfig { iterations: ROSENBROCK_BUDGET, initial_sigma: 1.0, adaptation_period: 30, step_increase: 1.22, seed, }; let mut opt = OnePlusOneEs::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn rosenbrock_bo(seed: u64) -> SoRun { let problem = rosenbrock_problem(); let bounds = RealBounds::new(vec![(-5.0, 10.0); ROSENBROCK_DIM]); let config = BayesianOptConfig { initial_samples: 10, iterations: 50, length_scales: None, signal_variance: 1.0, noise_variance: 1e-6, acquisition_samples: 1_000, seed, }; let mut opt = BayesianOpt::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn ackley_bo(seed: u64) -> SoRun { let problem = ackley_problem(); let bounds = RealBounds::new(vec![(-32.768, 32.768); ACKLEY_DIM]); let config = BayesianOptConfig { initial_samples: 10, iterations: 50, length_scales: None, signal_variance: 1.0, noise_variance: 1e-6, acquisition_samples: 1_000, seed, }; let mut opt = BayesianOpt::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } // ----------------------------------------------------------------------------- // Rosenbrock + Ackley runners (a curated SO subset on each) // ----------------------------------------------------------------------------- pub fn rosenbrock_problem() -> Rosenbrock { Rosenbrock { dim: ROSENBROCK_DIM, } } pub fn ackley_problem() -> Ackley { Ackley { dim: ACKLEY_DIM } } macro_rules! so_run_de { ($problem_expr:expr, $dim:expr, $bounds_lo:expr, $bounds_hi:expr, $budget:expr, $seed:expr) => {{ let problem = $problem_expr; let bounds = RealBounds::new(vec![($bounds_lo, $bounds_hi); $dim]); let pop = 50; let gens = ($budget - pop) / pop; let config = DifferentialEvolutionConfig { population_size: pop, generations: gens, differential_weight: 0.5, crossover_probability: 0.9, seed: $seed, }; let mut opt = DifferentialEvolution::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } }}; } macro_rules! so_run_cma { ($problem_expr:expr, $dim:expr, $bounds_lo:expr, $bounds_hi:expr, $budget:expr, $seed:expr) => {{ let problem = $problem_expr; let bounds = RealBounds::new(vec![($bounds_lo, $bounds_hi); $dim]); let pop = 16; let config = CmaEsConfig { population_size: pop, generations: $budget / pop, initial_sigma: 1.0, eigen_decomposition_period: 1, initial_mean: None, seed: $seed, }; let mut opt = CmaEs::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } }}; } macro_rules! so_run_pso { ($problem_expr:expr, $dim:expr, $bounds_lo:expr, $bounds_hi:expr, $budget:expr, $seed:expr) => {{ let problem = $problem_expr; let bounds = RealBounds::new(vec![($bounds_lo, $bounds_hi); $dim]); let swarm = 40; let config = ParticleSwarmConfig { swarm_size: swarm, generations: ($budget - 2 * swarm) / swarm, inertia: 0.7, cognitive: 1.5, social: 1.5, seed: $seed, }; let mut opt = ParticleSwarm::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } }}; } macro_rules! so_run_tlbo { ($problem_expr:expr, $dim:expr, $bounds_lo:expr, $bounds_hi:expr, $budget:expr, $seed:expr) => {{ let problem = $problem_expr; let bounds = RealBounds::new(vec![($bounds_lo, $bounds_hi); $dim]); let pop = 30; // TLBO does ~2N evaluations per generation. let gens = ($budget - pop) / (2 * pop); let config = TlboConfig { population_size: pop, generations: gens, seed: $seed, }; let mut opt = Tlbo::new(config, bounds); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } }}; } pub fn rosenbrock_de(seed: u64) -> SoRun { so_run_de!( rosenbrock_problem(), ROSENBROCK_DIM, -5.0, 10.0, ROSENBROCK_BUDGET, seed ) } pub fn rosenbrock_cma(seed: u64) -> SoRun { so_run_cma!( rosenbrock_problem(), ROSENBROCK_DIM, -5.0, 10.0, ROSENBROCK_BUDGET, seed ) } pub fn rosenbrock_pso(seed: u64) -> SoRun { so_run_pso!( rosenbrock_problem(), ROSENBROCK_DIM, -5.0, 10.0, ROSENBROCK_BUDGET, seed ) } pub fn rosenbrock_tlbo(seed: u64) -> SoRun { so_run_tlbo!( rosenbrock_problem(), ROSENBROCK_DIM, -5.0, 10.0, ROSENBROCK_BUDGET, seed ) } pub fn ackley_de(seed: u64) -> SoRun { so_run_de!( ackley_problem(), ACKLEY_DIM, -32.768, 32.768, ACKLEY_BUDGET, seed ) } pub fn ackley_cma(seed: u64) -> SoRun { so_run_cma!( ackley_problem(), ACKLEY_DIM, -32.768, 32.768, ACKLEY_BUDGET, seed ) } pub fn ackley_pso(seed: u64) -> SoRun { so_run_pso!( ackley_problem(), ACKLEY_DIM, -32.768, 32.768, ACKLEY_BUDGET, seed ) } pub fn ackley_tlbo(seed: u64) -> SoRun { so_run_tlbo!( ackley_problem(), ACKLEY_DIM, -32.768, 32.768, ACKLEY_BUDGET, seed ) } // ----------------------------------------------------------------------------- // ZDT3 runners (curated MO subset) // ----------------------------------------------------------------------------- pub fn zdt3_problem() -> Zdt3 { Zdt3 { dim: ZDT3_DIM } } pub fn zdt3_nsga2(seed: u64) -> MoRun { let problem = zdt3_problem(); let bounds = vec![(0.0, 1.0); ZDT3_DIM]; 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 / ZDT3_DIM as f64), }; let pop = 100; let config = Nsga2Config { population_size: pop, generations: ZDT3_BUDGET / pop, seed, }; let mut opt = Nsga2::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt3_moead(seed: u64) -> MoRun { let problem = zdt3_problem(); let bounds = vec![(0.0, 1.0); ZDT3_DIM]; 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 / ZDT3_DIM as f64), }; let pop = 100; let config = MoeadConfig { generations: (ZDT3_BUDGET - pop) / pop, reference_divisions: 99, neighborhood_size: 20, seed, }; let mut opt = Moead::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt3_ibea(seed: u64) -> MoRun { let problem = zdt3_problem(); let bounds = vec![(0.0, 1.0); ZDT3_DIM]; 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 / ZDT3_DIM as f64), }; let pop = 100; let config = IbeaConfig { population_size: pop, generations: ZDT3_BUDGET / pop, kappa: 0.05, seed, }; let mut opt = Ibea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt3_age_moea(seed: u64) -> MoRun { let problem = zdt3_problem(); let bounds = vec![(0.0, 1.0); ZDT3_DIM]; 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 / ZDT3_DIM as f64), }; let pop = 100; let config = AgeMoeaConfig { population_size: pop, generations: ZDT3_BUDGET / pop, seed, }; let mut opt = AgeMoea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn zdt3_knea(seed: u64) -> MoRun { let problem = zdt3_problem(); let bounds = vec![(0.0, 1.0); ZDT3_DIM]; 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 / ZDT3_DIM as f64), }; let pop = 100; let config = KneaConfig { population_size: pop, generations: ZDT3_BUDGET / pop, seed, }; let mut opt = Knea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } // ----------------------------------------------------------------------------- // DTLZ1 runners (curated many-obj subset) // ----------------------------------------------------------------------------- pub fn dtlz1_problem() -> Dtlz1 { Dtlz1 { num_objectives: DTLZ1_OBJECTIVES, dim: DTLZ1_DIM, } } pub fn dtlz1_nsga3(seed: u64) -> MoRun { let problem = dtlz1_problem(); let bounds = vec![(0.0, 1.0); DTLZ1_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ1_DIM as f64), }; let pop = 92; let config = Nsga3Config { population_size: pop, generations: DTLZ1_BUDGET / pop, reference_divisions: 12, seed, }; let mut opt = Nsga3::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz1_moead(seed: u64) -> MoRun { let problem = dtlz1_problem(); let bounds = vec![(0.0, 1.0); DTLZ1_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ1_DIM as f64), }; let pop = 91; let config = MoeadConfig { generations: (DTLZ1_BUDGET - pop) / pop, reference_divisions: 12, neighborhood_size: 20, seed, }; let mut opt = Moead::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz1_age_moea(seed: u64) -> MoRun { let problem = dtlz1_problem(); let bounds = vec![(0.0, 1.0); DTLZ1_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ1_DIM as f64), }; let pop = 92; let config = AgeMoeaConfig { population_size: pop, generations: DTLZ1_BUDGET / pop, seed, }; let mut opt = AgeMoea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn dtlz1_grea(seed: u64) -> MoRun { let problem = dtlz1_problem(); let bounds = vec![(0.0, 1.0); DTLZ1_DIM]; let initializer = RealBounds::new(bounds.clone()); let variation = CompositeVariation { crossover: SimulatedBinaryCrossover::new(bounds.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(bounds, 20.0, 1.0 / DTLZ1_DIM as f64), }; let pop = 92; let config = GreaConfig { population_size: pop, generations: DTLZ1_BUDGET / pop, grid_divisions: 8, seed, }; let mut opt = Grea::new(config, initializer, variation); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } /// Mean L2 distance from each front point to the analytical DTLZ1 front /// (`Σf_i = 0.5`, all `f_i ≥ 0`). Closed-form: signed distance from the /// hyperplane projected to non-negative. pub fn mean_distance_to_dtlz1_front(front: &[Candidate>]) -> f64 { if front.is_empty() { return f64::INFINITY; } let total: f64 = front .iter() .map(|c| { let s: f64 = c.evaluation.objectives.iter().sum(); (s - 0.5).abs() }) .sum(); total / front.len() as f64 } // ----------------------------------------------------------------------------- // Main // ----------------------------------------------------------------------------- pub fn run_zdt1_comparison() { println!("== ZDT1 (dim={ZDT1_DIM}, {ZDT1_BUDGET} evals/run × {SEEDS} seeds) =="); println!("Zitzler-Deb-Thiele 2-objective benchmark: {ZDT1_DIM} real variables, one smooth"); println!("convex Pareto front f2 = 1 - sqrt(f1). Hard because 29 of 30 variables must"); println!("collapse to 0 before the front is even reachable, and only then can the"); println!("population spread along it. Optimum: mean L2 -> 0 (the front is known exactly)."); println!("sorted best-first by hypervolume (higher better; spacing / mean L2 lower better)"); println!(); let zdt1 = Zdt1 { dim: ZDT1_DIM }; let zdt1_objs = zdt1.objectives(); type Runner = fn(u64) -> MoRun; let runners: &[(&str, Runner)] = &[ ("RandomSearch", zdt1_random), ("PAES", zdt1_paes), ("MOPSO", zdt1_mopso), ("SPEA2", zdt1_spea2), ("PESA-II", zdt1_pesa2), ("eps-MOEA", zdt1_epsilon_moea), ("IBEA", zdt1_ibea), ("HypE", zdt1_hype), ("SMS-EMOA", zdt1_sms_emoa), ("RVEA", zdt1_rvea), ("NSGA-II", zdt1_nsga2), ("NSGA-III", zdt1_nsga3), ("MOEA/D", zdt1_moead), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let hv: Vec = runs .iter() .map(|r| hypervolume_2d(&r.front, &zdt1_objs, ZDT1_REFERENCE)) .collect(); let sp: Vec = runs.iter().map(|r| spacing(&r.front, &zdt1_objs)).collect(); let l2: Vec = runs .iter() .map(|r| mean_l2_to_zdt1_front(&r.front)) .collect(); let fs: Vec = runs.iter().map(|r| r.front.len() as f64).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (hv_m, hv_s) = mean_std(&hv); let (sp_m, sp_s) = mean_std(&sp); let (l2_m, l2_s) = mean_std(&l2); let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); rows.push(( hv_m, vec![ name.to_string(), format!("{hv_m:.4}+/-{hv_s:.4}"), format!("{sp_m:.4}+/-{sp_s:.4}"), format!("{l2_m:.4}+/-{l2_s:.4}"), format!("{fs_m:.0}"), format!("{ms_m:.0}"), ], )); } // Higher hypervolume is better. rows.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table( &[ "algorithm", "hypervolume", "spacing", "mean L2", "front", "ms", ], &table, ); } pub fn run_dtlz2_comparison() { println!(); println!("== DTLZ2 (3-obj, dim={DTLZ2_DIM}, {DTLZ2_BUDGET} evals/run × {SEEDS} seeds) =="); println!("Deb-Thiele-Laumanns-Zitzler 3-objective; the Pareto front is the unit-sphere"); println!("octant (Σf² = 1, all f >= 0) -- a curved 2-D surface embedded in 3-D objective"); println!("space. Hard because crowding/spacing must work in a higher dimension."); println!("'mean dist' = |‖f‖ - 1|, so 0 means perfectly on the sphere (the known optimum)."); println!("sorted best-first by mean dist (lower is better)"); println!(); let dtlz2 = dtlz2_problem(); let dtlz2_objs = dtlz2.objectives(); type Runner = fn(u64) -> MoRun; let runners: &[(&str, Runner)] = &[ ("RandomSearch", dtlz2_random), ("MOPSO", dtlz2_mopso), ("NSGA-II", dtlz2_nsga2), ("SPEA2", dtlz2_spea2), ("PESA-II", dtlz2_pesa2), ("eps-MOEA", dtlz2_epsilon_moea), ("IBEA", dtlz2_ibea), ("HypE", dtlz2_hype), ("SMS-EMOA", dtlz2_sms_emoa), ("RVEA", dtlz2_rvea), ("NSGA-III", dtlz2_nsga3), ("MOEA/D", dtlz2_moead), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let dist: Vec = runs .iter() .map(|r| mean_distance_to_dtlz2_front(&r.front)) .collect(); let sp: Vec = runs .iter() .map(|r| spacing(&r.front, &dtlz2_objs)) .collect(); let fs: Vec = runs.iter().map(|r| r.front.len() as f64).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (d_m, d_s) = mean_std(&dist); let (sp_m, sp_s) = mean_std(&sp); let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); rows.push(( d_m, vec![ name.to_string(), format!("{d_m:.4}+/-{d_s:.4}"), format!("{sp_m:.4}+/-{sp_s:.4}"), format!("{fs_m:.0}"), format!("{ms_m:.0}"), ], )); } // Lower mean distance to the true front is better. rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table( &["algorithm", "mean dist", "spacing", "front", "ms"], &table, ); } pub fn run_rastrigin_comparison() { println!(); println!("== Rastrigin (dim={RASTRIGIN_DIM}, {RASTRIGIN_BUDGET} evals/run × {SEEDS} seeds) =="); println!( "Highly multimodal trap: f = 10n + Σ(x_i² - 10·cos(2π·x_i)), {RASTRIGIN_DIM} dimensions." ); println!( "Hard because a near-quadratic global bowl is overlaid with ~10^{RASTRIGIN_DIM} regularly" ); println!("spaced local minima -- any greedy step lands in the nearest dimple."); println!("Global optimum: f = 0 at the origin."); println!("sorted best-first (lower is better)"); println!(); type Runner = fn(u64) -> SoRun; let runners: &[(&str, Runner)] = &[ ("RandomSearch", rastrigin_random), ("HillClimber", rastrigin_hill_climber), ("(1+1)-ES", rastrigin_one_plus_one_es), ("SimulatedAnneal", rastrigin_simulated_annealing), ("PAES", rastrigin_paes), ("GA", rastrigin_genetic_algorithm), ("PSO", rastrigin_particle_swarm), ("NSGA-II", rastrigin_nsga2), ("DE", rastrigin_de), ("CMA-ES", rastrigin_cma_es), ("IPOP-CMA-ES", rastrigin_ipop_cma_es), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let best: Vec = runs.iter().map(|r| r.best_value).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (b_m, b_s) = mean_std(&best); let (ms_m, _) = mean_std(&ms); rows.push(( b_m, vec![ name.to_string(), format!("{b_m:.4e} +/- {b_s:.2e}"), format!("{ms_m:.0}"), ], )); } // Lower best f is better. rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table(&["algorithm", "best f", "ms"], &table); } pub fn run_rosenbrock_comparison() { println!(); println!("== Rosenbrock (dim={ROSENBROCK_DIM}, {ROSENBROCK_BUDGET} evals × {SEEDS} seeds) =="); println!( "Rosenbrock's banana valley: f = Σ(100·(x_{{i+1}} - x_i²)² + (1 - x_i)²), {ROSENBROCK_DIM} dims." ); println!("Hard because the minimum sits in a long, bent, near-flat valley -- easy to enter,"); println!("very slow to crawl along to the tip. Global optimum: f = 0 at the all-ones point."); println!("sorted best-first (lower is better)"); println!(); type Runner = fn(u64) -> SoRun; let runners: &[(&str, Runner)] = &[ ("DE", rosenbrock_de), ("PSO", rosenbrock_pso), ("CMA-ES", rosenbrock_cma), ("TLBO", rosenbrock_tlbo), ("(1+1)-ES", rosenbrock_one_plus_one_es), ("Nelder-Mead", rosenbrock_nelder_mead), ("BO (60 evals)", rosenbrock_bo), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let best: Vec = runs.iter().map(|r| r.best_value).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (b_m, b_s) = mean_std(&best); let (ms_m, _) = mean_std(&ms); rows.push(( b_m, vec![ name.to_string(), format!("{b_m:.4e} +/- {b_s:.2e}"), format!("{ms_m:.0}"), ], )); } rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table(&["algorithm", "best f", "ms"], &table); } pub fn run_ackley_comparison() { println!(); println!("== Ackley (dim={ACKLEY_DIM}, {ACKLEY_BUDGET} evals × {SEEDS} seeds) =="); println!("Ackley's function: a near-flat outer plateau with shallow ripples surrounding a"); println!( "single deep, narrow global basin, {ACKLEY_DIM} dimensions. Hard because the gradient is" ); println!("almost zero far from the optimum, giving local search little to follow."); println!("Global optimum: f = 0 at the origin."); println!("sorted best-first (lower is better)"); println!(); type Runner = fn(u64) -> SoRun; let runners: &[(&str, Runner)] = &[ ("DE", ackley_de), ("PSO", ackley_pso), ("CMA-ES", ackley_cma), ("TLBO", ackley_tlbo), ("BO (60 evals)", ackley_bo), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let best: Vec = runs.iter().map(|r| r.best_value).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (b_m, b_s) = mean_std(&best); let (ms_m, _) = mean_std(&ms); rows.push(( b_m, vec![ name.to_string(), format!("{b_m:.4e} +/- {b_s:.2e}"), format!("{ms_m:.0}"), ], )); } rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table(&["algorithm", "best f", "ms"], &table); } pub fn run_zdt3_comparison() { println!(); println!("== ZDT3 (dim={ZDT3_DIM}, {ZDT3_BUDGET} evals × {SEEDS} seeds) =="); println!("Zitzler-Deb-Thiele 2-objective with a DISCONNECTED front: five separate arcs"); println!("rather than one curve. Hard because an algorithm has to discover and populate"); println!("every arc while not stranding solutions in the dominated gaps between them."); println!("Optimum: cover all five arcs; scored by hypervolume vs [11, 11]."); println!("sorted best-first by hypervolume (higher better; spacing lower better)"); println!(); let problem = zdt3_problem(); let objs = problem.objectives(); type Runner = fn(u64) -> MoRun; let runners: &[(&str, Runner)] = &[ ("NSGA-II", zdt3_nsga2), ("MOEA/D", zdt3_moead), ("IBEA", zdt3_ibea), ("AGE-MOEA", zdt3_age_moea), ("KnEA", zdt3_knea), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let hv: Vec = runs .iter() .map(|r| hypervolume_2d(&r.front, &objs, ZDT3_REFERENCE)) .collect(); let sp: Vec = runs.iter().map(|r| spacing(&r.front, &objs)).collect(); let fs: Vec = runs.iter().map(|r| r.front.len() as f64).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (hv_m, hv_s) = mean_std(&hv); let (sp_m, sp_s) = mean_std(&sp); let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); rows.push(( hv_m, vec![ name.to_string(), format!("{hv_m:.4}+/-{hv_s:.4}"), format!("{sp_m:.4}+/-{sp_s:.4}"), format!("{fs_m:.0}"), format!("{ms_m:.0}"), ], )); } rows.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table( &["algorithm", "hypervolume", "spacing", "front", "ms"], &table, ); } pub fn run_dtlz1_comparison() { println!(); println!("== DTLZ1 (3-obj, dim={DTLZ1_DIM}, {DTLZ1_BUDGET} evals × {SEEDS} seeds) =="); println!("Deb-Thiele-Laumanns-Zitzler 3-objective; the Pareto front is the linear simplex"); println!("Σf = 0.5 in the positive octant. Hard because a deceptive multimodal 'g' term"); println!("riddles the approach with a huge number of local fronts -- only fully-converged"); println!("runs land on the simplex. Optimum: mean dist -> 0."); println!("sorted best-first by mean dist (lower is better)"); println!(); let problem = dtlz1_problem(); let objs = problem.objectives(); type Runner = fn(u64) -> MoRun; let runners: &[(&str, Runner)] = &[ ("NSGA-III", dtlz1_nsga3), ("MOEA/D", dtlz1_moead), ("AGE-MOEA", dtlz1_age_moea), ("GrEA", dtlz1_grea), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let dist: Vec = runs .iter() .map(|r| mean_distance_to_dtlz1_front(&r.front)) .collect(); let sp: Vec = runs.iter().map(|r| spacing(&r.front, &objs)).collect(); let fs: Vec = runs.iter().map(|r| r.front.len() as f64).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (d_m, d_s) = mean_std(&dist); let (sp_m, sp_s) = mean_std(&sp); let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); rows.push(( d_m, vec![ name.to_string(), format!("{d_m:.4}+/-{d_s:.4}"), format!("{sp_m:.4}+/-{sp_s:.4}"), format!("{fs_m:.0}"), format!("{ms_m:.0}"), ], )); } rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table( &["algorithm", "mean dist", "spacing", "front", "ms"], &table, ); } // ============================================================================= // Combinatorial / sequencing problems // // These use a different decision encoding — `Vec` permutations and // `Vec` bitstrings — than the continuous problems above, so the // applicable-algorithm roster differs: the real-vector methods (CMA-ES, DE, // PSO, …) cannot run here, while permutation-native methods (Ant Colony, // Tabu Search) and the EDA-style ones can. // ============================================================================= pub const TSP_CITIES: usize = 15; pub const TSP_BUDGET: usize = 8_000; pub const JSS_BUDGET: usize = 8_000; pub const KNAPSACK_BUDGET: usize = 20_000; // ---- Single-objective TSP: a regular polygon on the unit circle ------------- /// `TSP_CITIES` points equally spaced on the unit circle. The optimal tour /// just visits them in angular order; its length is the regular-polygon /// perimeter `2·n·sin(π/n)`, which gives an exact known-best baseline. pub struct RingTsp { distances: Vec>, } impl RingTsp { fn new() -> Self { let pts: Vec<(f64, f64)> = (0..TSP_CITIES) .map(|i| { let a = 2.0 * PI * (i as f64) / (TSP_CITIES as f64); (a.cos(), a.sin()) }) .collect(); let distances = (0..TSP_CITIES) .map(|i| { (0..TSP_CITIES) .map(|j| { let (xi, yi) = pts[i]; let (xj, yj) = pts[j]; ((xi - xj).powi(2) + (yi - yj).powi(2)).sqrt() }) .collect() }) .collect(); Self { distances } } /// Known optimal tour length — the regular-polygon perimeter. fn optimal_length() -> f64 { 2.0 * TSP_CITIES as f64 * (PI / TSP_CITIES as f64).sin() } } impl Problem for RingTsp { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![Objective::minimize("tour_length")]) } fn evaluate(&self, tour: &Vec) -> Evaluation { let n = tour.len(); let mut total = 0.0; for i in 0..n { total += self.distances[tour[i]][tour[(i + 1) % n]]; } Evaluation::new(vec![total]) } } pub fn tsp_random(seed: u64) -> SoRun { let problem = RingTsp::new(); let mut opt = RandomSearch::new( RandomSearchConfig { iterations: TSP_BUDGET, batch_size: 1, seed, }, ShuffledPermutation { n: TSP_CITIES }, ); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn tsp_hill_climber(seed: u64) -> SoRun { let problem = RingTsp::new(); let mut opt = HillClimber::new( HillClimberConfig { iterations: TSP_BUDGET, seed, }, ShuffledPermutation { n: TSP_CITIES }, InversionMutation, ); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn tsp_simulated_annealing(seed: u64) -> SoRun { let problem = RingTsp::new(); let mut opt = SimulatedAnnealing::new( SimulatedAnnealingConfig { iterations: TSP_BUDGET, initial_temperature: 1.0, final_temperature: 1e-3, seed, }, ShuffledPermutation { n: TSP_CITIES }, InversionMutation, ); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn tsp_tabu_search(seed: u64) -> SoRun { let problem = RingTsp::new(); // Each step considers 16 inversion (2-opt-style) neighbours. let neighbors = |tour: &Vec, rng: &mut Rng| { let mut m = InversionMutation; (0..16) .map(|_| m.vary(std::slice::from_ref(tour), rng).pop().unwrap()) .collect() }; let mut opt = TabuSearch::new( TabuSearchConfig { iterations: TSP_BUDGET / 16, tabu_tenure: 20, seed, }, ShuffledPermutation { n: TSP_CITIES }, neighbors, ); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn tsp_genetic_algorithm(seed: u64) -> SoRun { let problem = RingTsp::new(); let mut opt = GeneticAlgorithm::new( GeneticAlgorithmConfig { population_size: 80, generations: TSP_BUDGET / 80, tournament_size: 3, elitism: 2, seed, }, ShuffledPermutation { n: TSP_CITIES }, CompositeVariation { crossover: OrderCrossover, mutation: InversionMutation, }, ); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn tsp_ant_colony(seed: u64) -> SoRun { let problem = RingTsp::new(); let distances = problem.distances.clone(); let mut opt = AntColonyTsp::new( AntColonyTspConfig { ants: 20, generations: TSP_BUDGET / 20, alpha: 1.0, beta: 3.0, evaporation: 0.5, deposit: 1.0, initial_pheromone: 1.0, seed, }, distances, ); let t0 = Instant::now(); let result = opt.run(&problem); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } // ---- Single-objective job-shop scheduling: Fisher & Thompson FT06 ----------- pub const JSS_JOBS: usize = 6; pub const JSS_MACHINES: usize = 6; /// FT06 routing — machine id of the k-th operation of job j. pub const FT06_MACHINE: [[usize; JSS_MACHINES]; JSS_JOBS] = [ [2, 0, 1, 3, 5, 4], [1, 2, 4, 5, 0, 3], [2, 3, 5, 0, 1, 4], [1, 0, 2, 3, 4, 5], [2, 1, 4, 5, 0, 3], [1, 3, 5, 0, 4, 2], ]; /// FT06 processing times — duration of the k-th operation of job j. pub const FT06_TIME: [[f64; JSS_MACHINES]; JSS_JOBS] = [ [1.0, 3.0, 6.0, 7.0, 3.0, 6.0], [8.0, 5.0, 10.0, 10.0, 10.0, 4.0], [5.0, 4.0, 8.0, 9.0, 1.0, 7.0], [5.0, 5.0, 5.0, 3.0, 8.0, 9.0], [9.0, 3.0, 5.0, 4.0, 3.0, 1.0], [3.0, 3.0, 9.0, 10.0, 4.0, 1.0], ]; /// FT06's optimal makespan is a long-settled benchmark value. pub const FT06_OPTIMAL_MAKESPAN: f64 = 55.0; /// FT06 job-shop, makespan objective. Operation-string encoding: a length-36 /// multiset where job id `j` appears `JSS_MACHINES` times; the k-th /// occurrence of `j` is its k-th operation. pub struct Ft06Makespan; impl Problem for Ft06Makespan { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![Objective::minimize("makespan")]) } fn evaluate(&self, schedule: &Vec) -> Evaluation { let mut job_next = [0_usize; JSS_JOBS]; let mut job_clock = [0.0_f64; JSS_JOBS]; let mut machine_clock = [0.0_f64; JSS_MACHINES]; for &job in schedule { let k = job_next[job]; let m = FT06_MACHINE[job][k]; let t = FT06_TIME[job][k]; let start = job_clock[job].max(machine_clock[m]); let end = start + t; job_clock[job] = end; machine_clock[m] = end; job_next[job] = k + 1; } let makespan = machine_clock.iter().cloned().fold(0.0_f64, f64::max); Evaluation::new(vec![makespan]) } } /// Precedence-Order Crossover — multiset-preserving recombination for the /// operation-string encoding (the strict-permutation crossovers would /// corrupt the multiset). #[derive(Debug, Clone, Copy, Default)] pub struct PrecedenceOrderCrossover; impl Variation> for PrecedenceOrderCrossover { fn vary(&mut self, parents: &[Vec], rng: &mut Rng) -> Vec> { assert!(parents.len() >= 2, "POX requires 2 parents"); let (p1, p2) = (&parents[0], &parents[1]); let mut in_j1 = [false; JSS_JOBS]; loop { for slot in &mut in_j1 { *slot = rng.random_bool(0.5); } let c = in_j1.iter().filter(|&&b| b).count(); if c > 0 && c < JSS_JOBS { break; } } vec![pox_child(p1, p2, &in_j1), pox_child(p2, p1, &in_j1)] } } pub fn pox_child(donor: &[usize], filler: &[usize], in_donor_set: &[bool]) -> Vec { let n = donor.len(); let mut child = vec![usize::MAX; n]; for k in 0..n { if in_donor_set[donor[k]] { child[k] = donor[k]; } } let mut fill_idx = 0; for &v in filler { if !in_donor_set[v] { while fill_idx < n && child[fill_idx] != usize::MAX { fill_idx += 1; } child[fill_idx] = v; fill_idx += 1; } } child } pub fn jss_initializer() -> ShuffledMultisetPermutation { ShuffledMultisetPermutation::new(vec![JSS_MACHINES; JSS_JOBS]) } pub fn jss_random(seed: u64) -> SoRun { let mut opt = RandomSearch::new( RandomSearchConfig { iterations: JSS_BUDGET, batch_size: 1, seed, }, jss_initializer(), ); let t0 = Instant::now(); let result = opt.run(&Ft06Makespan); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn jss_hill_climber(seed: u64) -> SoRun { let mut opt = HillClimber::new( HillClimberConfig { iterations: JSS_BUDGET, seed, }, jss_initializer(), InsertionMutation, ); let t0 = Instant::now(); let result = opt.run(&Ft06Makespan); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn jss_simulated_annealing(seed: u64) -> SoRun { let mut opt = SimulatedAnnealing::new( SimulatedAnnealingConfig { iterations: JSS_BUDGET, initial_temperature: 5.0, final_temperature: 1e-2, seed, }, jss_initializer(), InsertionMutation, ); let t0 = Instant::now(); let result = opt.run(&Ft06Makespan); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn jss_tabu_search(seed: u64) -> SoRun { let neighbors = |schedule: &Vec, rng: &mut Rng| { let mut m = InsertionMutation; (0..16) .map(|_| m.vary(std::slice::from_ref(schedule), rng).pop().unwrap()) .collect() }; let mut opt = TabuSearch::new( TabuSearchConfig { iterations: JSS_BUDGET / 16, tabu_tenure: 20, seed, }, jss_initializer(), neighbors, ); let t0 = Instant::now(); let result = opt.run(&Ft06Makespan); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } pub fn jss_genetic_algorithm(seed: u64) -> SoRun { let mut opt = GeneticAlgorithm::new( GeneticAlgorithmConfig { population_size: 80, generations: JSS_BUDGET / 80, tournament_size: 3, elitism: 2, seed, }, jss_initializer(), CompositeVariation { crossover: PrecedenceOrderCrossover, mutation: InsertionMutation, }, ); let t0 = Instant::now(); let result = opt.run(&Ft06Makespan); SoRun { best_value: result.best.unwrap().evaluation.objectives[0], wall_ms: t0.elapsed().as_millis(), } } // ---- Bi-objective 0/1 knapsack: Zitzler & Thiele style ---------------------- pub const KNAPSACK_N: usize = 30; pub const KP_PROFIT_A: [f64; KNAPSACK_N] = [ 61.0, 17.0, 92.0, 49.0, 73.0, 28.0, 84.0, 36.0, 55.0, 78.0, 23.0, 91.0, 12.0, 67.0, 45.0, 58.0, 33.0, 71.0, 14.0, 26.0, 87.0, 42.0, 19.0, 65.0, 30.0, 51.0, 79.0, 22.0, 47.0, 88.0, ]; pub const KP_PROFIT_B: [f64; KNAPSACK_N] = [ 24.0, 81.0, 16.0, 67.0, 29.0, 73.0, 41.0, 60.0, 52.0, 19.0, 77.0, 34.0, 95.0, 22.0, 71.0, 88.0, 56.0, 27.0, 64.0, 90.0, 18.0, 43.0, 79.0, 31.0, 85.0, 25.0, 38.0, 92.0, 70.0, 13.0, ]; pub const KP_WEIGHT: [f64; KNAPSACK_N] = [ 35.0, 58.0, 22.0, 71.0, 14.0, 86.0, 31.0, 53.0, 78.0, 19.0, 44.0, 16.0, 67.0, 88.0, 25.0, 51.0, 33.0, 74.0, 12.0, 47.0, 63.0, 28.0, 91.0, 36.0, 55.0, 17.0, 82.0, 41.0, 24.0, 68.0, ]; /// Hypervolume reference point for the knapsack front. Both objectives are /// maximized, so in the minimization-oriented frame any positive profit is /// below 0 — `[0, 0]` is dominated by every feasible solution. pub const KNAPSACK_REFERENCE: [f64; 2] = [0.0, 0.0]; /// Bi-objective 0/1 knapsack: two profit vectors, one shared capacity. /// Weight overruns are penalized in both objectives so the Pareto front is /// composed of feasible solutions. pub struct BiKnapsack { capacity: f64, } impl BiKnapsack { fn new() -> Self { Self { capacity: 0.5 * KP_WEIGHT.iter().sum::(), } } } impl Problem for BiKnapsack { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![ Objective::maximize("profit_a"), Objective::maximize("profit_b"), ]) } fn evaluate(&self, take: &Vec) -> Evaluation { let (mut pa, mut pb, mut w) = (0.0, 0.0, 0.0); for (i, &t) in take.iter().enumerate() { if t { pa += KP_PROFIT_A[i]; pb += KP_PROFIT_B[i]; w += KP_WEIGHT[i]; } } let penalty = 1000.0 * (w - self.capacity).max(0.0); Evaluation::new(vec![pa - penalty, pb - penalty]) } } /// Random binary initializer — each bit 50/50 independently. #[derive(Debug, Clone, Copy)] pub struct RandomBinary { n: usize, } impl Initializer> for RandomBinary { fn initialize(&mut self, size: usize, rng: &mut Rng) -> Vec> { (0..size) .map(|_| (0..self.n).map(|_| rng.random_bool(0.5)).collect()) .collect() } } /// One-point crossover for binary chromosomes. #[derive(Debug, Clone, Copy, Default)] pub struct OnePointCrossoverBool; impl Variation> for OnePointCrossoverBool { fn vary(&mut self, parents: &[Vec], rng: &mut Rng) -> Vec> { assert!( parents.len() >= 2, "OnePointCrossoverBool requires 2 parents" ); let (p1, p2) = (&parents[0], &parents[1]); let n = p1.len(); if n < 2 { return vec![p1.clone(), p2.clone()]; } let cut = rng.random_range(1..n); let mut c1 = Vec::with_capacity(n); let mut c2 = Vec::with_capacity(n); c1.extend_from_slice(&p1[..cut]); c1.extend_from_slice(&p2[cut..]); c2.extend_from_slice(&p2[..cut]); c2.extend_from_slice(&p1[cut..]); vec![c1, c2] } } pub fn knap_binary_variation() -> CompositeVariation { CompositeVariation { crossover: OnePointCrossoverBool, mutation: BitFlipMutation { probability: 1.0 / KNAPSACK_N as f64, }, } } /// One run's knapsack metrics. The bi-objective front lives over `Vec` /// decisions, so `MoRun` (which fixes `Vec`) does not fit — we keep /// only the aggregate quality numbers. #[derive(Clone)] pub struct KnapRun { hypervolume: f64, front_size: usize, wall_ms: u128, } pub fn knap_run(front: &[Candidate>], problem: &BiKnapsack, wall_ms: u128) -> KnapRun { KnapRun { hypervolume: hypervolume_2d(front, &problem.objectives(), KNAPSACK_REFERENCE), front_size: front.len(), wall_ms, } } pub fn knapsack_random(seed: u64) -> KnapRun { let problem = BiKnapsack::new(); let mut opt = RandomSearch::new( RandomSearchConfig { iterations: KNAPSACK_BUDGET, batch_size: 1, seed, }, RandomBinary { n: KNAPSACK_N }, ); let t0 = Instant::now(); let result = opt.run(&problem); let ms = t0.elapsed().as_millis(); knap_run(&result.pareto_front, &problem, ms) } pub fn knapsack_nsga2(seed: u64) -> KnapRun { let problem = BiKnapsack::new(); let mut opt = Nsga2::new( Nsga2Config { population_size: 100, generations: KNAPSACK_BUDGET / 100, seed, }, RandomBinary { n: KNAPSACK_N }, knap_binary_variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); let ms = t0.elapsed().as_millis(); knap_run(&result.pareto_front, &problem, ms) } pub fn knapsack_spea2(seed: u64) -> KnapRun { let problem = BiKnapsack::new(); let mut opt = Spea2::new( Spea2Config { population_size: 100, archive_size: 100, generations: KNAPSACK_BUDGET / 100, seed, }, RandomBinary { n: KNAPSACK_N }, knap_binary_variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); let ms = t0.elapsed().as_millis(); knap_run(&result.pareto_front, &problem, ms) } pub fn knapsack_nsga3(seed: u64) -> KnapRun { let problem = BiKnapsack::new(); let mut opt = Nsga3::new( Nsga3Config { population_size: 100, generations: KNAPSACK_BUDGET / 100, reference_divisions: 20, seed, }, RandomBinary { n: KNAPSACK_N }, knap_binary_variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); let ms = t0.elapsed().as_millis(); knap_run(&result.pareto_front, &problem, ms) } pub fn knapsack_ibea(seed: u64) -> KnapRun { let problem = BiKnapsack::new(); let mut opt = Ibea::new( IbeaConfig { population_size: 100, generations: KNAPSACK_BUDGET / 100, kappa: 0.05, seed, }, RandomBinary { n: KNAPSACK_N }, knap_binary_variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); let ms = t0.elapsed().as_millis(); knap_run(&result.pareto_front, &problem, ms) } // ---- Combinatorial comparison runners -------------------------------------- pub fn run_tsp_comparison() { println!(); println!("== TSP ring-{TSP_CITIES} ({TSP_BUDGET} evals/run × {SEEDS} seeds) =="); println!( "{TSP_CITIES} equally-spaced cities on the unit circle; minimize the closed tour length." ); println!( "The space is ({TSP_CITIES}-1)!/2 distinct tours, but cities in convex position have no" ); println!("2-opt local optima -- so this instance cleanly separates methods with good"); println!("neighbourhood moves (inversion = 2-opt) from blind recombination / sampling."); println!( "Known optimum (the polygon perimeter): {:.4}", RingTsp::optimal_length() ); println!("sorted best-first by tour length (lower is better)"); println!(); type Runner = fn(u64) -> SoRun; let runners: &[(&str, Runner)] = &[ ("RandomSearch", tsp_random), ("HillClimber", tsp_hill_climber), ("SimulatedAnneal", tsp_simulated_annealing), ("TabuSearch", tsp_tabu_search), ("GA", tsp_genetic_algorithm), ("AntColony", tsp_ant_colony), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let best: Vec = runs.iter().map(|r| r.best_value).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (b_m, b_s) = mean_std(&best); let (ms_m, _) = mean_std(&ms); rows.push(( b_m, vec![ name.to_string(), format!("{b_m:.4}+/-{b_s:.4}"), format!("{ms_m:.0}"), ], )); } rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table(&["algorithm", "tour length", "ms"], &table); } pub fn run_jss_comparison() { println!(); println!("== JSS FT06 ({JSS_BUDGET} evals/run × {SEEDS} seeds) =="); println!("Fisher & Thompson 1963 6-job × 6-machine job-shop; minimize makespan."); println!("Hard because every job has a fixed machine order, so swapping two"); println!( "operations can ripple delays across the whole schedule. Known optimum: {FT06_OPTIMAL_MAKESPAN:.0}" ); println!("sorted best-first by makespan (lower is better)"); println!(); type Runner = fn(u64) -> SoRun; let runners: &[(&str, Runner)] = &[ ("RandomSearch", jss_random), ("HillClimber", jss_hill_climber), ("SimulatedAnneal", jss_simulated_annealing), ("TabuSearch", jss_tabu_search), ("GA", jss_genetic_algorithm), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let best: Vec = runs.iter().map(|r| r.best_value).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (b_m, b_s) = mean_std(&best); let (ms_m, _) = mean_std(&ms); rows.push(( b_m, vec![ name.to_string(), format!("{b_m:.4}+/-{b_s:.4}"), format!("{ms_m:.0}"), ], )); } rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table(&["algorithm", "makespan", "ms"], &table); } pub fn run_knapsack_comparison() { println!(); println!( "== Knapsack ({KNAPSACK_N} items, bi-objective, {KNAPSACK_BUDGET} evals/run × {SEEDS} seeds) ==" ); println!("Zitzler-Thiele style 0/1 knapsack: two profit vectors, one capacity"); println!("(half the total weight). Hard because the two profit objectives"); println!("conflict and the capacity constraint carves feasible regions out of"); println!("the 2^{KNAPSACK_N} bitstrings. No closed-form optimum; scored by hypervolume"); println!("vs reference {KNAPSACK_REFERENCE:?} (higher is better)."); println!("sorted best-first by hypervolume"); println!(); type Runner = fn(u64) -> KnapRun; let runners: &[(&str, Runner)] = &[ ("RandomSearch", knapsack_random), ("NSGA-II", knapsack_nsga2), ("SPEA2", knapsack_spea2), ("NSGA-III", knapsack_nsga3), ("IBEA", knapsack_ibea), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let hv: Vec = runs.iter().map(|r| r.hypervolume).collect(); let fs: Vec = runs.iter().map(|r| r.front_size as f64).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (hv_m, hv_s) = mean_std(&hv); let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); rows.push(( hv_m, vec![ name.to_string(), format!("{hv_m:.1}+/-{hv_s:.1}"), format!("{fs_m:.0}"), format!("{ms_m:.0}"), ], )); } rows.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table(&["algorithm", "hypervolume", "front", "ms"], &table); } // ============================================================================= // Many-objective problems (4+ objectives) // // The "curse of dimensionality" for multi-objective optimizers: as the // objective count climbs, the fraction of mutually non-dominated solution // pairs rushes toward 1, so Pareto rank alone stops discriminating. // NSGA-II's whole population collapses into front 0 and only crowding // distance is left to steer. Reference-point (NSGA-III), decomposition // (MOEA/D), reference-vector (RVEA), grid (GrEA) and indicator (IBEA, // HypE) methods are built to survive this regime. // // Both DTLZ structs above are already parameterised by objective count, // and the DTLZ1/DTLZ2 distance metrics generalise to any M, so a single // `ManySpec` + generic runners cover every objective count. // ============================================================================= /// A DTLZ instance at an arbitrary objective count, plus the budget and /// the reference-set sizing the many-objective algorithms need. #[derive(Clone, Copy)] pub struct ManySpec { objectives: usize, dim: usize, budget: usize, /// Population for the fixed-population algorithms; also the /// Das-Dennis weight count MOEA/D derives from `reference_divisions`. population: usize, /// Das-Dennis divisions for NSGA-III / RVEA / MOEA/D reference sets. reference_divisions: usize, /// `true` = DTLZ1 (deceptive multimodal, linear-simplex front); /// `false` = DTLZ2 (unit-hypersphere-octant front). is_dtlz1: bool, /// Per-axis HypE hypervolume reference coordinate. hype_ref: f64, } /// One DTLZ problem type so the generic runners have a single `Problem` /// to hand to `run` regardless of which front geometry is in play. pub enum ManyDtlz { D1(Dtlz1), D2(Dtlz2), } impl Problem for ManyDtlz { type Decision = Vec; fn objectives(&self) -> ObjectiveSpace { match self { ManyDtlz::D1(p) => p.objectives(), ManyDtlz::D2(p) => p.objectives(), } } fn evaluate(&self, x: &Vec) -> Evaluation { match self { ManyDtlz::D1(p) => p.evaluate(x), ManyDtlz::D2(p) => p.evaluate(x), } } } impl ManySpec { fn problem(&self) -> ManyDtlz { if self.is_dtlz1 { ManyDtlz::D1(Dtlz1 { num_objectives: self.objectives, dim: self.dim, }) } else { ManyDtlz::D2(Dtlz2 { num_objectives: self.objectives, dim: self.dim, }) } } fn bounds(&self) -> Vec<(f64, f64)> { vec![(0.0, 1.0); self.dim] } fn variation(&self) -> CompositeVariation { let b = self.bounds(); CompositeVariation { crossover: SimulatedBinaryCrossover::new(b.clone(), 30.0, 1.0), mutation: PolynomialMutation::new(b, 20.0, 1.0 / self.dim as f64), } } fn mean_dist(&self, front: &[Candidate>]) -> f64 { if self.is_dtlz1 { mean_distance_to_dtlz1_front(front) } else { mean_distance_to_dtlz2_front(front) } } } pub fn many_random(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = RandomSearch::new( RandomSearchConfig { iterations: spec.budget, batch_size: 1, seed, }, RealBounds::new(spec.bounds()), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn many_nsga2(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = Nsga2::new( Nsga2Config { population_size: spec.population, generations: spec.budget / spec.population, seed, }, RealBounds::new(spec.bounds()), spec.variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn many_nsga3(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = Nsga3::new( Nsga3Config { population_size: spec.population, generations: spec.budget / spec.population, reference_divisions: spec.reference_divisions, seed, }, RealBounds::new(spec.bounds()), spec.variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn many_moead(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = Moead::new( MoeadConfig { generations: spec.budget / spec.population, reference_divisions: spec.reference_divisions, neighborhood_size: 20, seed, }, RealBounds::new(spec.bounds()), spec.variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn many_rvea(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = Rvea::new( RveaConfig { population_size: spec.population, generations: spec.budget / spec.population, reference_divisions: spec.reference_divisions, alpha: 2.0, seed, }, RealBounds::new(spec.bounds()), spec.variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn many_grea(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = Grea::new( GreaConfig { population_size: spec.population, generations: spec.budget / spec.population, grid_divisions: 10, seed, }, RealBounds::new(spec.bounds()), spec.variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn many_ibea(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = Ibea::new( IbeaConfig { population_size: spec.population, generations: spec.budget / spec.population, kappa: 0.05, seed, }, RealBounds::new(spec.bounds()), spec.variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn many_hype(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = Hype::new( HypeConfig { population_size: spec.population, generations: spec.budget / spec.population, reference_point: vec![spec.hype_ref; spec.objectives], mc_samples: 1_000, seed, }, RealBounds::new(spec.bounds()), spec.variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn many_age_moea(spec: ManySpec, seed: u64) -> MoRun { let problem = spec.problem(); let mut opt = AgeMoea::new( AgeMoeaConfig { population_size: spec.population, generations: spec.budget / spec.population, seed, }, RealBounds::new(spec.bounds()), spec.variation(), ); let t0 = Instant::now(); let result = opt.run(&problem); MoRun { front: result.pareto_front, wall_ms: t0.elapsed().as_millis(), } } pub fn run_many_objective_comparison(title: &str, blurb: &[&str], spec: ManySpec) { println!(); println!("== {title} =="); for line in blurb { println!("{line}"); } println!("sorted best-first by mean dist to the true front (lower is better)"); println!(); type Runner = fn(ManySpec, u64) -> MoRun; let runners: &[(&str, Runner)] = &[ ("RandomSearch", many_random), ("NSGA-II", many_nsga2), ("NSGA-III", many_nsga3), ("MOEA/D", many_moead), ("RVEA", many_rvea), ("GrEA", many_grea), ("IBEA", many_ibea), ("HypE", many_hype), ("AGE-MOEA", many_age_moea), ]; let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(|s| runner(spec, s)).collect(); let dist: Vec = runs.iter().map(|r| spec.mean_dist(&r.front)).collect(); let fs: Vec = runs.iter().map(|r| r.front.len() as f64).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (d_m, d_s) = mean_std(&dist); let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); rows.push(( d_m, vec![ name.to_string(), format!("{d_m:.4}+/-{d_s:.4}"), format!("{fs_m:.0}"), format!("{ms_m:.0}"), ], )); } rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); print_table(&["algorithm", "mean dist", "front", "ms"], &table); } /// Run every algorithm runner exactly once at seed 0 — the full `compare` /// workload at a single seed, with no multi-seed aggregation and no /// printing. This is the unit of work the `compare_profile` gungraun /// benchmark profiles under callgrind. /// /// Returns a checksum folded from every run's result so the optimizer /// cannot elide any runner call. /// /// Must enumerate *every* runner function — `benches/compare_profile.rs` /// allow-lists the presentation layer, so a runner missing here would not /// be caught by the bench's dead-code lint (it is caught by the example /// build, where every runner is reachable through `run_all`). #[allow(dead_code)] // entry point for the compare_profile bench; unused by the example pub fn profile_workload() -> u64 { let mut acc: u64 = 0; macro_rules! mo { ($run:expr) => { acc = acc.wrapping_add($run.front.len() as u64); }; } macro_rules! so { ($run:expr) => { acc ^= $run.best_value.to_bits(); }; } macro_rules! kn { ($run:expr) => {{ let r = $run; acc = acc.wrapping_add(r.front_size as u64) ^ r.hypervolume.to_bits(); }}; } // ZDT1 — 2-objective continuous mo!(zdt1_random(0)); mo!(zdt1_paes(0)); mo!(zdt1_spea2(0)); mo!(zdt1_nsga2(0)); mo!(zdt1_sms_emoa(0)); mo!(zdt1_hype(0)); mo!(zdt1_rvea(0)); mo!(zdt1_pesa2(0)); mo!(zdt1_epsilon_moea(0)); mo!(zdt1_mopso(0)); mo!(zdt1_ibea(0)); mo!(zdt1_moead(0)); mo!(zdt1_nsga3(0)); // DTLZ2 — 3-objective continuous mo!(dtlz2_random(0)); mo!(dtlz2_nsga2(0)); mo!(dtlz2_spea2(0)); mo!(dtlz2_sms_emoa(0)); mo!(dtlz2_hype(0)); mo!(dtlz2_rvea(0)); mo!(dtlz2_pesa2(0)); mo!(dtlz2_epsilon_moea(0)); mo!(dtlz2_mopso(0)); mo!(dtlz2_ibea(0)); mo!(dtlz2_moead(0)); mo!(dtlz2_nsga3(0)); // Rastrigin — single-objective multimodal so!(rastrigin_random(0)); so!(rastrigin_paes(0)); so!(rastrigin_nsga2(0)); so!(rastrigin_de(0)); so!(rastrigin_hill_climber(0)); so!(rastrigin_simulated_annealing(0)); so!(rastrigin_genetic_algorithm(0)); so!(rastrigin_particle_swarm(0)); so!(rastrigin_cma_es(0)); so!(rastrigin_ipop_cma_es(0)); so!(rastrigin_one_plus_one_es(0)); // Rosenbrock + Ackley — single-objective smooth / multimodal so!(rosenbrock_nelder_mead(0)); so!(rosenbrock_one_plus_one_es(0)); so!(rosenbrock_bo(0)); so!(ackley_bo(0)); so!(rosenbrock_de(0)); so!(rosenbrock_cma(0)); so!(rosenbrock_pso(0)); so!(rosenbrock_tlbo(0)); so!(ackley_de(0)); so!(ackley_cma(0)); so!(ackley_pso(0)); so!(ackley_tlbo(0)); // ZDT3 — 2-objective, disconnected front mo!(zdt3_nsga2(0)); mo!(zdt3_moead(0)); mo!(zdt3_ibea(0)); mo!(zdt3_age_moea(0)); mo!(zdt3_knea(0)); // DTLZ1 — 3-objective, linear simplex mo!(dtlz1_nsga3(0)); mo!(dtlz1_moead(0)); mo!(dtlz1_age_moea(0)); mo!(dtlz1_grea(0)); // TSP ring — single-objective permutation so!(tsp_random(0)); so!(tsp_hill_climber(0)); so!(tsp_simulated_annealing(0)); so!(tsp_tabu_search(0)); so!(tsp_genetic_algorithm(0)); so!(tsp_ant_colony(0)); // JSS FT06 — single-objective sequencing so!(jss_random(0)); so!(jss_hill_climber(0)); so!(jss_simulated_annealing(0)); so!(jss_tabu_search(0)); so!(jss_genetic_algorithm(0)); // Knapsack — bi-objective binary kn!(knapsack_random(0)); kn!(knapsack_nsga2(0)); kn!(knapsack_spea2(0)); kn!(knapsack_nsga3(0)); kn!(knapsack_ibea(0)); // Many-objective — every runner against each of the three DTLZ specs let many_specs = [ ManySpec { objectives: 4, dim: 13, budget: 40_000, population: 56, reference_divisions: 5, is_dtlz1: false, hype_ref: 3.0, }, ManySpec { objectives: 10, dim: 19, budget: 40_000, population: 55, reference_divisions: 2, is_dtlz1: false, hype_ref: 3.0, }, ManySpec { objectives: 8, dim: 12, budget: 40_000, population: 120, reference_divisions: 3, is_dtlz1: true, hype_ref: 1.0, }, ]; for spec in many_specs { mo!(many_random(spec, 0)); mo!(many_nsga2(spec, 0)); mo!(many_nsga3(spec, 0)); mo!(many_moead(spec, 0)); mo!(many_rvea(spec, 0)); mo!(many_grea(spec, 0)); mo!(many_ibea(spec, 0)); mo!(many_hype(spec, 0)); mo!(many_age_moea(spec, 0)); } acc } pub fn run_all() { run_zdt1_comparison(); run_zdt3_comparison(); run_dtlz2_comparison(); run_dtlz1_comparison(); run_rastrigin_comparison(); run_rosenbrock_comparison(); run_ackley_comparison(); run_tsp_comparison(); run_jss_comparison(); run_knapsack_comparison(); // ---- Many-objective (4+) ---- run_many_objective_comparison( "DTLZ2 4-objective (dim=13, 40000 evals/run × 10 seeds)", &[ "DTLZ2 scaled to 4 objectives -- the entry point to many-objective.", "Front is still the unit-hypersphere octant (Σf² = 1). Already hard:", "with 4 objectives most random pairs of solutions are mutually", "non-dominated, so Pareto rank alone barely discriminates. Optimum:", "mean dist -> 0.", ], ManySpec { objectives: 4, dim: 13, budget: 40_000, population: 56, reference_divisions: 5, is_dtlz1: false, hype_ref: 3.0, }, ); run_many_objective_comparison( "DTLZ2 10-objective (dim=19, 40000 evals/run × 10 seeds)", &[ "DTLZ2 scaled to 10 objectives -- the curse of dimensionality in full.", "In 10-D objective space almost EVERY pair of solutions is mutually", "non-dominated, so NSGA-II's whole population collapses into front 0", "and crowding distance is the only signal left. Reference-point,", "decomposition and indicator methods are built for exactly this.", "Watch the 'front' column: it pins to the population size because", "nothing dominates anything. Optimum: mean dist -> 0.", ], ManySpec { objectives: 10, dim: 19, budget: 40_000, population: 55, reference_divisions: 2, is_dtlz1: false, hype_ref: 3.0, }, ); run_many_objective_comparison( "DTLZ1 8-objective (dim=12, 40000 evals/run × 10 seeds)", &[ "DTLZ1 scaled to 8 objectives -- the brutal one. Stacks the", "many-objective dominance collapse on top of DTLZ1's deceptive", "multimodal g-term (a huge number of local fronts). The true front", "is the linear simplex Σf = 0.5; reaching it at all is the", "achievement. Expect large mean-dist values and wide spreads -- this", "is near the edge of what the catalogue does at this budget.", ], ManySpec { objectives: 8, dim: 12, budget: 40_000, population: 120, reference_divisions: 3, is_dtlz1: true, hype_ref: 1.0, }, ); }