diff --git a/examples/compare.rs b/examples/compare.rs index 772326f..25e072a 100644 --- a/examples/compare.rs +++ b/examples/compare.rs @@ -32,6 +32,21 @@ const DTLZ2_K: usize = 10; const DTLZ2_DIM: usize = DTLZ2_OBJECTIVES + DTLZ2_K - 1; // 12 const DTLZ2_BUDGET: usize = 30_000; +const ROSENBROCK_DIM: usize = 5; +const ROSENBROCK_BUDGET: usize = 30_000; + +const ACKLEY_DIM: usize = 5; +const ACKLEY_BUDGET: usize = 30_000; + +const ZDT3_DIM: usize = 30; +const ZDT3_BUDGET: usize = 25_000; +const ZDT3_REFERENCE: [f64; 2] = [11.0, 11.0]; + +const DTLZ1_OBJECTIVES: usize = 3; +const DTLZ1_K: usize = 5; +const DTLZ1_DIM: usize = DTLZ1_OBJECTIVES + DTLZ1_K - 1; +const DTLZ1_BUDGET: usize = 30_000; + // ----------------------------------------------------------------------------- // Test problems // ----------------------------------------------------------------------------- @@ -92,6 +107,113 @@ impl Problem for Dtlz2 { } } +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]) + } +} + +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]) + } +} + +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]) + } +} + +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) + } +} + struct Rastrigin { dim: usize, } @@ -892,6 +1014,292 @@ fn rastrigin_cma_es(seed: u64) -> SoRun { } } +// ----------------------------------------------------------------------------- +// Rosenbrock + Ackley runners (a curated SO subset on each) +// ----------------------------------------------------------------------------- + +fn rosenbrock_problem() -> Rosenbrock { + Rosenbrock { dim: ROSENBROCK_DIM } +} +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, + 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(), + } + }}; +} + +fn rosenbrock_de(seed: u64) -> SoRun { so_run_de!(rosenbrock_problem(), ROSENBROCK_DIM, -5.0, 10.0, ROSENBROCK_BUDGET, seed) } +fn rosenbrock_cma(seed: u64) -> SoRun { so_run_cma!(rosenbrock_problem(), ROSENBROCK_DIM, -5.0, 10.0, ROSENBROCK_BUDGET, seed) } +fn rosenbrock_pso(seed: u64) -> SoRun { so_run_pso!(rosenbrock_problem(), ROSENBROCK_DIM, -5.0, 10.0, ROSENBROCK_BUDGET, seed) } +fn rosenbrock_tlbo(seed: u64) -> SoRun { so_run_tlbo!(rosenbrock_problem(), ROSENBROCK_DIM, -5.0, 10.0, ROSENBROCK_BUDGET, seed) } + +fn ackley_de(seed: u64) -> SoRun { so_run_de!(ackley_problem(), ACKLEY_DIM, -32.768, 32.768, ACKLEY_BUDGET, seed) } +fn ackley_cma(seed: u64) -> SoRun { so_run_cma!(ackley_problem(), ACKLEY_DIM, -32.768, 32.768, ACKLEY_BUDGET, seed) } +fn ackley_pso(seed: u64) -> SoRun { so_run_pso!(ackley_problem(), ACKLEY_DIM, -32.768, 32.768, ACKLEY_BUDGET, seed) } +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) +// ----------------------------------------------------------------------------- + +fn zdt3_problem() -> Zdt3 { Zdt3 { dim: ZDT3_DIM } } + +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() } +} + +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() } +} + +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() } +} + +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() } +} + +// ----------------------------------------------------------------------------- +// DTLZ1 runners (curated many-obj subset) +// ----------------------------------------------------------------------------- + +fn dtlz1_problem() -> Dtlz1 { + Dtlz1 { num_objectives: DTLZ1_OBJECTIVES, dim: DTLZ1_DIM } +} + +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() } +} + +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() } +} + +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() } +} + +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. +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 // ----------------------------------------------------------------------------- @@ -1056,8 +1464,150 @@ fn run_rastrigin_comparison() { } } +fn run_rosenbrock_comparison() { + println!(); + println!("== Rosenbrock (dim={ROSENBROCK_DIM}, {ROSENBROCK_BUDGET} evals × {SEEDS} seeds) =="); + println!("smooth non-convex valley; global minimum f = 0 at all-ones"); + println!(); + println!("{:<14} {:>20} {:>10}", "algorithm", "best f", "ms"); + println!("{}", "-".repeat(48)); + type Runner = fn(u64) -> SoRun; + let runners: &[(&str, Runner)] = &[ + ("DE", rosenbrock_de), + ("PSO", rosenbrock_pso), + ("CMA-ES", rosenbrock_cma), + ("TLBO", rosenbrock_tlbo), + ]; + 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); + println!( + "{:<14} {:>20} {:>10}", + name, + format!("{b_m:.4e} ± {b_s:.2e}"), + format!("{ms_m:.0}"), + ); + } +} + +fn run_ackley_comparison() { + println!(); + println!("== Ackley (dim={ACKLEY_DIM}, {ACKLEY_BUDGET} evals × {SEEDS} seeds) =="); + println!("smoother multimodal landscape than Rastrigin; global minimum f = 0 at origin"); + println!(); + println!("{:<14} {:>20} {:>10}", "algorithm", "best f", "ms"); + println!("{}", "-".repeat(48)); + type Runner = fn(u64) -> SoRun; + let runners: &[(&str, Runner)] = &[ + ("DE", ackley_de), + ("PSO", ackley_pso), + ("CMA-ES", ackley_cma), + ("TLBO", ackley_tlbo), + ]; + 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); + println!( + "{:<14} {:>20} {:>10}", + name, + format!("{b_m:.4e} ± {b_s:.2e}"), + format!("{ms_m:.0}"), + ); + } +} + +fn run_zdt3_comparison() { + println!(); + println!("== ZDT3 (dim={ZDT3_DIM}, {ZDT3_BUDGET} evals × {SEEDS} seeds) =="); + println!("disconnected Pareto front (not contiguous); spread across gaps matters"); + println!(); + println!( + "{:<14} {:>16} {:>14} {:>10} {:>10}", + "algorithm", "hypervolume↑", "spacing↓", "front", "ms", + ); + println!("{}", "-".repeat(70)); + 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), + ]; + 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); + println!( + "{:<14} {:>16} {:>14} {:>10} {:>10}", + name, + format!("{hv_m:.4}±{hv_s:.4}"), + format!("{sp_m:.4}±{sp_s:.4}"), + format!("{fs_m:.0}"), + format!("{ms_m:.0}"), + ); + } +} + +fn run_dtlz1_comparison() { + println!(); + println!("== DTLZ1 (3-obj, dim={DTLZ1_DIM}, {DTLZ1_BUDGET} evals × {SEEDS} seeds) =="); + println!("Pareto front: linear simplex Σf=0.5 in the positive octant"); + println!(); + println!( + "{:<14} {:>16} {:>14} {:>10} {:>10}", + "algorithm", "mean dist↓", "spacing↓", "front", "ms", + ); + println!("{}", "-".repeat(70)); + 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), + ]; + 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); + println!( + "{:<14} {:>16} {:>14} {:>10} {:>10}", + name, + format!("{d_m:.4}±{d_s:.4}"), + format!("{sp_m:.4}±{sp_s:.4}"), + format!("{fs_m:.0}"), + format!("{ms_m:.0}"), + ); + } +} + fn main() { run_zdt1_comparison(); + run_zdt3_comparison(); run_dtlz2_comparison(); + run_dtlz1_comparison(); run_rastrigin_comparison(); + run_rosenbrock_comparison(); + run_ackley_comparison(); }