From ab545e268a627895a0af30ad49faa7a985433f3b Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Wed, 13 May 2026 18:37:47 -0600 Subject: [PATCH] docs(examples): add bi-objective TSP crossover-comparison demo tsp_operators_compare.rs runs NSGA-II four times on the KroAB-25 bi-objective TSP, holding everything constant except the crossover operator. Ranks OX, PMX, CX, ERX by hypervolume against a fixed reference point, plus front size, unique-point count, and runtime. Pedagogical demonstration that the right comparison metric for a Pareto search is hypervolume, not single-objective fitness. --- examples/tsp_operators_compare.rs | 227 ++++++++++++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100644 examples/tsp_operators_compare.rs diff --git a/examples/tsp_operators_compare.rs b/examples/tsp_operators_compare.rs new file mode 100644 index 0000000..3ecad5e --- /dev/null +++ b/examples/tsp_operators_compare.rs @@ -0,0 +1,227 @@ +//! Crossover showdown on the bi-objective TSP from `btsp_kroab.rs`. +//! +//! Runs NSGA-II four times on the same KroAB-25 instance, holding everything +//! constant except the **crossover** operator. The mutation +//! ([`InversionMutation`]), initializer, population, generations, and seed +//! are identical across runs. +//! +//! Operators compared: +//! - [`OrderCrossover`] (OX) +//! - [`PartiallyMappedCrossover`] (PMX) +//! - [`CycleCrossover`] (CX) +//! - [`EdgeRecombinationCrossover`] (ERX) +//! +//! Each run is ranked by **hypervolume** (the standard Pareto-front quality +//! metric), not by single-objective fitness — for a Pareto search, "best +//! length on A" or "best length on B" alone is a misleading scoreboard. +//! +//! Run with: +//! +//! ```bash +//! cargo run --release --example tsp_operators_compare +//! ``` + +use heuropt::metrics::hypervolume_2d; +use heuropt::prelude::*; +use std::time::Instant; + +/// First 25 cities of TSPLIB KroA100 (EUC_2D). +const KROA_25: [(f64, f64); 25] = [ + (1380.0, 939.0), (2848.0, 96.0), (3510.0, 1671.0), (457.0, 334.0), + (3888.0, 666.0), (984.0, 965.0), (2721.0, 1482.0), (1286.0, 525.0), + (2716.0, 1432.0),(738.0, 1325.0), (1251.0, 1832.0), (2728.0, 1698.0), + (3815.0, 169.0), (3683.0, 1533.0),(1247.0, 1945.0), (123.0, 862.0), + (1234.0, 1946.0),(252.0, 1240.0), (611.0, 673.0), (2576.0, 1676.0), + (928.0, 1700.0), (53.0, 857.0), (1807.0, 1711.0), (274.0, 1420.0), + (2574.0, 946.0), +]; + +/// First 25 cities of TSPLIB KroB100 (EUC_2D). +const KROB_25: [(f64, f64); 25] = [ + (3140.0, 1401.0),(556.0, 1056.0), (3675.0, 1522.0), (1182.0, 1853.0), + (3595.0, 1340.0),(1936.0, 953.0), (2722.0, 1311.0), (2839.0, 2055.0), + (2253.0, 1242.0),(3142.0, 1591.0),(627.0, 1336.0), (936.0, 211.0), + (4014.0, 471.0), (1376.0, 1452.0),(3289.0, 593.0), (1453.0, 67.0), + (1014.0, 1944.0),(2811.0, 1080.0),(3010.0, 1290.0), (1817.0, 1517.0), + (510.0, 458.0), (1717.0, 1693.0),(1252.0, 1633.0), (1693.0, 1374.0), + (539.0, 1378.0), +]; + +const N_CITIES: usize = 25; +const REF_POINT: [f64; 2] = [40_000.0, 40_000.0]; + +fn euc2d_matrix(coords: &[(f64, f64)]) -> Vec> { + let n = coords.len(); + let mut d = vec![vec![0.0_f64; n]; n]; + for i in 0..n { + for j in (i + 1)..n { + let dx = coords[i].0 - coords[j].0; + let dy = coords[i].1 - coords[j].1; + let dij = (dx * dx + dy * dy).sqrt().round(); + d[i][j] = dij; + d[j][i] = dij; + } + } + d +} + +struct BTsp { + dist_a: Vec>, + dist_b: Vec>, +} + +impl BTsp { + fn new() -> Self { + Self { + dist_a: euc2d_matrix(&KROA_25), + dist_b: euc2d_matrix(&KROB_25), + } + } + fn tour_length(d: &[Vec], tour: &[usize]) -> f64 { + let n = tour.len(); + let mut total = 0.0; + for i in 0..n { + total += d[tour[i]][tour[(i + 1) % n]]; + } + total + } +} + +impl Problem for BTsp { + type Decision = Vec; + fn objectives(&self) -> ObjectiveSpace { + ObjectiveSpace::new(vec![ + Objective::minimize("length_A"), + Objective::minimize("length_B"), + ]) + } + fn evaluate(&self, tour: &Vec) -> Evaluation { + Evaluation::new(vec![ + Self::tour_length(&self.dist_a, tour), + Self::tour_length(&self.dist_b, tour), + ]) + } +} + +struct RunSummary { + name: &'static str, + front_size: usize, + front_unique: usize, + corner_a: (f64, f64), + corner_b: (f64, f64), + hypervolume: f64, + seconds: f64, +} + +fn run_once(name: &'static str, problem: &BTsp, crossover: C) -> RunSummary +where + C: Variation>, +{ + let mut optimizer = Nsga2::new( + Nsga2Config { + population_size: 200, + generations: 500, + seed: 11, + }, + ShuffledPermutation { n: N_CITIES }, + CompositeVariation { + crossover, + mutation: InversionMutation, + }, + ); + let t0 = Instant::now(); + let result = optimizer.run(problem); + let seconds = t0.elapsed().as_secs_f64(); + + let mut front: Vec<&Candidate>> = result.pareto_front.iter().collect(); + front.sort_by(|a, b| { + a.evaluation.objectives[0] + .partial_cmp(&b.evaluation.objectives[0]) + .unwrap_or(std::cmp::Ordering::Equal) + }); + + let mut seen: Vec<(i64, i64)> = Vec::new(); + for c in &front { + let o = &c.evaluation.objectives; + let k = (o[0] as i64, o[1] as i64); + if !seen.contains(&k) { + seen.push(k); + } + } + + let corner_a = front + .first() + .map(|c| (c.evaluation.objectives[0], c.evaluation.objectives[1])) + .unwrap_or((f64::NAN, f64::NAN)); + let corner_b = front + .last() + .map(|c| (c.evaluation.objectives[0], c.evaluation.objectives[1])) + .unwrap_or((f64::NAN, f64::NAN)); + + let owned: Vec>> = result.pareto_front.to_vec(); + let hv = hypervolume_2d(&owned, &problem.objectives(), REF_POINT); + + RunSummary { + name, + front_size: result.pareto_front.len(), + front_unique: seen.len(), + corner_a, + corner_b, + hypervolume: hv, + seconds, + } +} + +fn main() { + let problem = BTsp::new(); + println!("Bi-objective TSP (KroAB-25): NSGA-II crossover showdown"); + println!("Same population, generations, seed across all runs."); + println!("Mutation held constant at InversionMutation."); + println!("Reference point for hypervolume: ({:.0}, {:.0})", REF_POINT[0], REF_POINT[1]); + println!(); + + let runs = vec![ + run_once("Order (OX)", &problem, OrderCrossover), + run_once("PartiallyMapped (PMX)", &problem, PartiallyMappedCrossover), + run_once("Cycle (CX)", &problem, CycleCrossover), + run_once("EdgeRecomb (ERX)", &problem, EdgeRecombinationCrossover), + ]; + + println!( + " {:<24} | {:>5} {:>5} | {:>17} | {:>17} | {:>14} | {:>6}", + "crossover", + "size", + "uniq", + "A-corner (A, B)", + "B-corner (A, B)", + "hypervolume", + "time" + ); + println!(" {}", "-".repeat(106)); + for r in &runs { + println!( + " {:<24} | {:>5} {:>5} | ({:>6.0},{:>6.0}) | ({:>6.0},{:>6.0}) | {:>14.0} | {:>5.2}s", + r.name, + r.front_size, + r.front_unique, + r.corner_a.0, + r.corner_a.1, + r.corner_b.0, + r.corner_b.1, + r.hypervolume, + r.seconds, + ); + } + println!(); + + // Pick the winner by hypervolume (largest dominated area = best front). + let winner = runs + .iter() + .max_by(|a, b| { + a.hypervolume + .partial_cmp(&b.hypervolume) + .unwrap_or(std::cmp::Ordering::Equal) + }) + .expect("non-empty runs"); + println!("Best by hypervolume: {} ({:.0})", winner.name, winner.hypervolume); +}