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

@@ -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<Vec<f64>> {
+    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<Vec<f64>>,
+    dist_b: Vec<Vec<f64>>,
+}
+
+impl BTsp {
+    fn new() -> Self {
+        Self {
+            dist_a: euc2d_matrix(&KROA_25),
+            dist_b: euc2d_matrix(&KROB_25),
+        }
+    }
+    fn tour_length(d: &[Vec<f64>], 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<usize>;
+    fn objectives(&self) -> ObjectiveSpace {
+        ObjectiveSpace::new(vec![
+            Objective::minimize("length_A"),
+            Objective::minimize("length_B"),
+        ])
+    }
+    fn evaluate(&self, tour: &Vec<usize>) -> 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<C>(name: &'static str, problem: &BTsp, crossover: C) -> RunSummary
+where
+    C: Variation<Vec<usize>>,
+{
+    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<Vec<usize>>> = 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<Candidate<Vec<usize>>> = 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);
+}