Files
heuropt/examples/compare.rs
T
swaits 5728ee14e2 feat(examples): add DTLZ2 (3-obj) and NSGA-III to comparison harness
NSGA-III's value over NSGA-II shows up at 3+ objectives, where
crowding distance loses its diversity signal. Adds a third comparison
section to `examples/compare.rs`:

DTLZ2 (3-objective, 12-D, the textbook benchmark for many-objective
algorithms): unit-sphere-octant Pareto front. Compares RandomSearch,
NSGA-II, SPEA2, and NSGA-III on:

- mean distance from front points to the unit sphere
  (closed-form: |1 - sqrt(f1² + f2² + f3²)|),
- spacing,
- front size,
- wall-clock ms.

NSGA-III config: H=12 reference divisions (91 reference points,
matching the canonical setup from Deb & Jain 2014).

Also wires NSGA-III into the existing ZDT1 (2-objective) section even
though it's not its sweet spot — useful as a regression check that the
algorithm at least keeps up with NSGA-II on bi-objective problems.
2026-05-04 19:57:21 -06:00

590 lines
19 KiB
Rust
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
//! Multi-seed algorithm comparison harness.
//!
//! Runs every applicable optimizer on each test problem across N seeds and
//! prints aggregate quality metrics. Adding a new algorithm to the
//! comparison is a single-line edit to the runner table — see the bottom
//! of this file.
//!
//! ```bash
//! cargo run --release --example compare
//! ```
use std::f64::consts::PI;
use std::time::Instant;
use heuropt::metrics::{hypervolume::hypervolume_2d, spacing::spacing};
use heuropt::prelude::*;
const SEEDS: u64 = 10;
const ZDT1_DIM: usize = 30;
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.
const ZDT1_REFERENCE: [f64; 2] = [11.0, 11.0];
const RASTRIGIN_DIM: usize = 5;
const RASTRIGIN_BUDGET: usize = 50_000;
const DTLZ2_OBJECTIVES: usize = 3;
const DTLZ2_K: usize = 10;
const DTLZ2_DIM: usize = DTLZ2_OBJECTIVES + DTLZ2_K - 1; // 12
const DTLZ2_BUDGET: usize = 30_000;
// -----------------------------------------------------------------------------
// Test problems
// -----------------------------------------------------------------------------
struct Zdt1 {
dim: usize,
}
impl Problem for Zdt1 {
type Decision = Vec<f64>;
fn objectives(&self) -> ObjectiveSpace {
ObjectiveSpace::new(vec![Objective::minimize("f1"), Objective::minimize("f2")])
}
fn evaluate(&self, x: &Vec<f64>) -> 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])
}
}
struct Dtlz2 {
num_objectives: usize,
dim: usize,
}
impl Problem for Dtlz2 {
type Decision = Vec<f64>;
fn objectives(&self) -> ObjectiveSpace {
ObjectiveSpace::new(
(0..self.num_objectives)
.map(|i| Objective::minimize(format!("f{}", i + 1)))
.collect(),
)
}
fn evaluate(&self, x: &Vec<f64>) -> 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;
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)
}
}
struct Rastrigin {
dim: usize,
}
impl Problem for Rastrigin {
type Decision = Vec<f64>;
fn objectives(&self) -> ObjectiveSpace {
ObjectiveSpace::new(vec![Objective::minimize("f")])
}
fn evaluate(&self, x: &Vec<f64>) -> 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::<f64>();
Evaluation::new(vec![value])
}
}
// -----------------------------------------------------------------------------
// Run results + metrics aggregation
// -----------------------------------------------------------------------------
#[derive(Clone)]
struct MoRun {
front: Vec<Candidate<Vec<f64>>>,
wall_ms: u128,
}
#[derive(Clone)]
struct SoRun {
best_value: f64,
wall_ms: u128,
}
fn mean_l2_to_zdt1_front(front: &[Candidate<Vec<f64>>]) -> 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
}
fn mean_std(values: &[f64]) -> (f64, f64) {
let n = values.len() as f64;
let mean = values.iter().sum::<f64>() / n;
let var = values.iter().map(|v| (v - mean).powi(2)).sum::<f64>() / n;
(mean, var.sqrt())
}
// -----------------------------------------------------------------------------
// ZDT1 algorithm runners
// -----------------------------------------------------------------------------
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() }
}
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() }
}
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() }
}
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() }
}
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)
// -----------------------------------------------------------------------------
fn dtlz2_problem() -> Dtlz2 {
Dtlz2 { num_objectives: DTLZ2_OBJECTIVES, dim: DTLZ2_DIM }
}
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() }
}
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() }
}
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() }
}
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|`.
fn mean_distance_to_dtlz2_front(front: &[Candidate<Vec<f64>>]) -> 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::<f64>().sqrt();
(norm - 1.0).abs()
})
.sum();
total / front.len() as f64
}
// -----------------------------------------------------------------------------
// Rastrigin algorithm runners
// -----------------------------------------------------------------------------
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(),
}
}
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(),
}
}
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(),
}
}
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(),
}
}
// -----------------------------------------------------------------------------
// Main
// -----------------------------------------------------------------------------
fn run_zdt1_comparison() {
println!(
"== ZDT1 (dim={ZDT1_DIM}, {ZDT1_BUDGET} evals/run × {SEEDS} seeds) =="
);
println!("metric arrows: hypervolume↑ (higher better), others↓ (lower better)");
println!();
println!(
"{:<14} {:>16} {:>14} {:>14} {:>10} {:>10}",
"algorithm", "hypervolume", "spacing", "mean L2", "front", "ms",
);
println!("{}", "-".repeat(82));
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),
("SPEA2", zdt1_spea2),
("NSGA-II", zdt1_nsga2),
("NSGA-III", zdt1_nsga3),
];
for (name, runner) in runners {
let runs: Vec<MoRun> = (0..SEEDS).map(runner).collect();
let hv: Vec<f64> = runs
.iter()
.map(|r| hypervolume_2d(&r.front, &zdt1_objs, ZDT1_REFERENCE))
.collect();
let sp: Vec<f64> =
runs.iter().map(|r| spacing(&r.front, &zdt1_objs)).collect();
let l2: Vec<f64> =
runs.iter().map(|r| mean_l2_to_zdt1_front(&r.front)).collect();
let fs: Vec<f64> = runs.iter().map(|r| r.front.len() as f64).collect();
let ms: Vec<f64> = 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);
println!(
"{:<14} {:>16} {:>14} {:>14} {:>10} {:>10}",
name,
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}"),
);
}
}
fn run_dtlz2_comparison() {
println!();
println!(
"== DTLZ2 (3-obj, dim={DTLZ2_DIM}, {DTLZ2_BUDGET} evals/run × {SEEDS} seeds) =="
);
println!("Pareto front: unit sphere octant (Σf²=1, all f≥0); 'mean dist' is |‖f‖−1|");
println!();
println!(
"{:<14} {:>16} {:>14} {:>10} {:>10}",
"algorithm", "mean dist↓", "spacing↓", "front", "ms",
);
println!("{}", "-".repeat(70));
let dtlz2 = dtlz2_problem();
let dtlz2_objs = dtlz2.objectives();
type Runner = fn(u64) -> MoRun;
let runners: &[(&str, Runner)] = &[
("RandomSearch", dtlz2_random),
("NSGA-II", dtlz2_nsga2),
("SPEA2", dtlz2_spea2),
("NSGA-III", dtlz2_nsga3),
];
for (name, runner) in runners {
let runs: Vec<MoRun> = (0..SEEDS).map(runner).collect();
let dist: Vec<f64> =
runs.iter().map(|r| mean_distance_to_dtlz2_front(&r.front)).collect();
let sp: Vec<f64> =
runs.iter().map(|r| spacing(&r.front, &dtlz2_objs)).collect();
let fs: Vec<f64> = runs.iter().map(|r| r.front.len() as f64).collect();
let ms: Vec<f64> = 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 run_rastrigin_comparison() {
println!();
println!(
"== Rastrigin (dim={RASTRIGIN_DIM}, {RASTRIGIN_BUDGET} evals/run × {SEEDS} seeds) =="
);
println!("global minimum: f = 0 (lower is better)");
println!();
println!("{:<14} {:>20} {:>10}", "algorithm", "best f", "ms");
println!("{}", "-".repeat(48));
type Runner = fn(u64) -> SoRun;
let runners: &[(&str, Runner)] = &[
("RandomSearch", rastrigin_random),
("PAES", rastrigin_paes),
("NSGA-II", rastrigin_nsga2),
("DE", rastrigin_de),
];
for (name, runner) in runners {
let runs: Vec<SoRun> = (0..SEEDS).map(runner).collect();
let best: Vec<f64> = runs.iter().map(|r| r.best_value).collect();
let ms: Vec<f64> = 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 main() {
run_zdt1_comparison();
run_dtlz2_comparison();
run_rastrigin_comparison();
}