test(fuzz): add cargo-fuzz harness for Pareto and operator hot paths

fuzz/fuzz_targets/clamp_to_bounds.rs

@@ -0,0 +1,88 @@
+#![no_main]
+//! Fuzz `ClampToBounds` + `ProjectToSimplex` repair operators for
+//! idempotence and target-set membership.
+
+use arbitrary::Arbitrary;
+use libfuzzer_sys::fuzz_target;
+
+use heuropt::prelude::*;
+
+#[derive(Arbitrary, Debug)]
+struct Input {
+    bounds: Vec<(f64, f64)>,
+    x: Vec<f64>,
+    simplex_total: f64,
+}
+
+fuzz_target!(|input: Input| {
+    if input.bounds.is_empty() || input.bounds.len() > 16 {
+        return;
+    }
+    if input.x.len() != input.bounds.len() {
+        return;
+    }
+    let bounds: Vec<(f64, f64)> = input
+        .bounds
+        .iter()
+        .filter_map(|&(lo, hi)| {
+            if lo.is_finite() && hi.is_finite() && lo < hi {
+                Some((lo, hi))
+            } else {
+                None
+            }
+        })
+        .collect();
+    if bounds.len() != input.bounds.len() {
+        return;
+    }
+    // Restrict to a numerically-reasonable magnitude range for repair
+    // operators — they are invoked downstream of evolutionary search where
+    // candidate magnitudes are bounded.
+    if input.x.iter().any(|v| !v.is_finite() || v.abs() > 1e30) {
+        return;
+    }
+
+    let mut x = input.x.clone();
+    let mut clamp = ClampToBounds::new(bounds.clone());
+    clamp.repair(&mut x);
+    for (j, &v) in x.iter().enumerate() {
+        let (lo, hi) = bounds[j];
+        assert!(v >= lo && v <= hi, "clamp out of bounds");
+    }
+    let after_one = x.clone();
+    clamp.repair(&mut x);
+    assert_eq!(x, after_one, "clamp not idempotent");
+
+    // Simplex projection only meaningful when total > 0 and dim >= 1.
+    // The Duchi/Held-Wolfe projection loses precision when |x| ≫ total
+    // (τ becomes indistinguishable from max(x) in f64). Restrict to inputs
+    // within the algorithm's well-conditioned regime, |x_i| ≤ total · 1e6.
+    let max_abs = input.x.iter().fold(0.0_f64, |a, &b| a.max(b.abs()));
+    if input.simplex_total.is_finite()
+        && input.simplex_total > 1.0
+        && input.simplex_total < 1e9
+        && max_abs <= input.simplex_total * 1e6
+    {
+        let mut y = input.x.clone();
+        let mut proj = ProjectToSimplex::new(input.simplex_total);
+        proj.repair(&mut y);
+        for &v in &y {
+            assert!(v >= 0.0, "project negative entry");
+        }
+        let s: f64 = y.iter().sum();
+        assert!(
+            (s - input.simplex_total).abs() < 1e-6 * input.simplex_total.max(1.0),
+            "project sum {s} != target {}",
+            input.simplex_total,
+        );
+        let after = y.clone();
+        proj.repair(&mut y);
+        for (a, b) in after.iter().zip(y.iter()) {
+            let scale = a.abs().max(b.abs()).max(1.0);
+            assert!(
+                (a - b).abs() < 1e-9 * scale,
+                "project not idempotent: {a} vs {b}",
+            );
+        }
+    }
+});

fuzz/fuzz_targets/crowding_distance.rs

@@ -0,0 +1,50 @@
+#![no_main]
+//! Fuzz `crowding_distance` for shape and non-negativity.
+
+use arbitrary::Arbitrary;
+use libfuzzer_sys::fuzz_target;
+
+use heuropt::core::candidate::Candidate;
+use heuropt::core::evaluation::Evaluation;
+use heuropt::core::objective::{Objective, ObjectiveSpace};
+use heuropt::pareto::crowding::crowding_distance;
+
+#[derive(Arbitrary, Debug)]
+struct Input {
+    points: Vec<(f64, f64)>,
+}
+
+fuzz_target!(|input: Input| {
+    if input.points.len() > 64 {
+        return;
+    }
+    // Bound magnitudes — crowding's `(max - min)` and per-axis gaps can
+    // both overflow to +∞ when points span ±f64::MAX, yielding inf/inf=NaN.
+    if input
+        .points
+        .iter()
+        .any(|&(a, b)| !a.is_finite() || !b.is_finite() || a.abs() > 1e150 || b.abs() > 1e150)
+    {
+        return;
+    }
+    let space = ObjectiveSpace::new(vec![Objective::minimize("f1"), Objective::minimize("f2")]);
+    let pop: Vec<Candidate<()>> = input
+        .points
+        .iter()
+        .map(|&(a, b)| Candidate::new((), Evaluation::new(vec![a, b])))
+        .collect();
+
+    let front: Vec<usize> = (0..pop.len()).collect();
+    let d = crowding_distance(&pop, &front, &space);
+    assert_eq!(d.len(), front.len(), "crowding distance length mismatch");
+    for (i, &v) in d.iter().enumerate() {
+        assert!(v >= 0.0 || v.is_infinite(), "negative crowding[{i}] = {v}");
+        assert!(!v.is_nan(), "NaN crowding[{i}]");
+    }
+    // If size <= 2, every entry is +∞.
+    if pop.len() <= 2 {
+        for (i, &v) in d.iter().enumerate() {
+            assert!(v.is_infinite(), "size<=2 crowding[{i}] not inf: {v}");
+        }
+    }
+});

fuzz/fuzz_targets/hypervolume_2d.rs

@@ -0,0 +1,47 @@
+#![no_main]
+//! Fuzz `hypervolume_2d` for non-negativity and reference-point handling.
+
+use arbitrary::Arbitrary;
+use libfuzzer_sys::fuzz_target;
+
+use heuropt::core::candidate::Candidate;
+use heuropt::core::evaluation::Evaluation;
+use heuropt::core::objective::{Objective, ObjectiveSpace};
+use heuropt::metrics::hypervolume::hypervolume_2d;
+
+#[derive(Arbitrary, Debug)]
+struct Input {
+    points: Vec<(f64, f64)>,
+    ref_point: (f64, f64),
+}
+
+fuzz_target!(|input: Input| {
+    if input.points.len() > 64 {
+        return;
+    }
+    // Non-finite floats are permitted by Evaluation, but HV is undefined
+    // there — restrict to finite for this property.
+    if !input.ref_point.0.is_finite() || !input.ref_point.1.is_finite() {
+        return;
+    }
+    if input
+        .points
+        .iter()
+        .any(|&(a, b)| !a.is_finite() || !b.is_finite())
+    {
+        return;
+    }
+
+    let space = ObjectiveSpace::new(vec![Objective::minimize("f1"), Objective::minimize("f2")]);
+    let pop: Vec<Candidate<()>> = input
+        .points
+        .iter()
+        .map(|&(a, b)| Candidate::new((), Evaluation::new(vec![a, b])))
+        .collect();
+    let hv = hypervolume_2d(&pop, &space, [input.ref_point.0, input.ref_point.1]);
+    // HV can be +∞ when the dominated rectangle area overflows f64 (e.g. a
+    // ref point at f64::MAX with deeply negative front coords). The
+    // contracted invariants are non-negativity and non-NaN.
+    assert!(hv >= 0.0, "HV negative: {hv}");
+    assert!(!hv.is_nan(), "HV is NaN");
+});

fuzz/fuzz_targets/non_dominated_sort.rs

@@ -0,0 +1,62 @@
+#![no_main]
+//! Fuzz `non_dominated_sort` for partition correctness.
+//!
+//! Invariants checked:
+//!   * Every population index appears in exactly one front.
+//!   * Earlier fronts dominate later fronts (no backwards domination).
+//!   * No panics on any vector of finite or non-finite objective values.
+
+use arbitrary::Arbitrary;
+use libfuzzer_sys::fuzz_target;
+
+use heuropt::core::candidate::Candidate;
+use heuropt::core::evaluation::Evaluation;
+use heuropt::core::objective::{Objective, ObjectiveSpace};
+use heuropt::pareto::dominance::{Dominance, pareto_compare};
+use heuropt::pareto::sort::non_dominated_sort;
+
+#[derive(Arbitrary, Debug)]
+struct Input {
+    objectives: Vec<(f64, f64)>,
+}
+
+fuzz_target!(|input: Input| {
+    if input.objectives.is_empty() || input.objectives.len() > 32 {
+        return;
+    }
+    let space = ObjectiveSpace::new(vec![Objective::minimize("f1"), Objective::minimize("f2")]);
+    let pop: Vec<Candidate<()>> = input
+        .objectives
+        .iter()
+        .map(|&(a, b)| Candidate::new((), Evaluation::new(vec![a, b])))
+        .collect();
+
+    let fronts = non_dominated_sort(&pop, &space);
+
+    // Partition: every index appears exactly once.
+    let mut seen = vec![false; pop.len()];
+    for front in &fronts {
+        for &idx in front {
+            assert!(!seen[idx], "index {idx} in multiple fronts");
+            seen[idx] = true;
+        }
+    }
+    for (i, &was) in seen.iter().enumerate() {
+        assert!(was, "index {i} missing from all fronts");
+    }
+
+    // Earlier fronts cannot be dominated by later fronts.
+    for (k, fk) in fronts.iter().enumerate() {
+        for fl in fronts.iter().skip(k + 1) {
+            for &i in fk {
+                for &j in fl {
+                    let r = pareto_compare(&pop[i].evaluation, &pop[j].evaluation, &space);
+                    assert!(
+                        !matches!(r, Dominance::DominatedBy),
+                        "front-{k}/{i} dominated by later front",
+                    );
+                }
+            }
+        }
+    }
+});

fuzz/fuzz_targets/pareto_archive.rs

@@ -0,0 +1,58 @@
+#![no_main]
+//! Fuzz `ParetoArchive` for the non-domination invariant under arbitrary
+//! insertion/truncation sequences.
+
+use arbitrary::Arbitrary;
+use libfuzzer_sys::fuzz_target;
+
+use heuropt::core::candidate::Candidate;
+use heuropt::core::evaluation::Evaluation;
+use heuropt::core::objective::{Objective, ObjectiveSpace};
+use heuropt::pareto::archive::ParetoArchive;
+use heuropt::pareto::dominance::{Dominance, pareto_compare};
+
+#[derive(Arbitrary, Debug)]
+enum Op {
+    Insert(f64, f64),
+    Truncate(u8),
+}
+
+#[derive(Arbitrary, Debug)]
+struct Input {
+    ops: Vec<Op>,
+}
+
+fuzz_target!(|input: Input| {
+    if input.ops.len() > 64 {
+        return;
+    }
+    let space = ObjectiveSpace::new(vec![Objective::minimize("f1"), Objective::minimize("f2")]);
+    let mut archive: ParetoArchive<()> = ParetoArchive::new(space.clone());
+
+    for op in input.ops {
+        match op {
+            Op::Insert(a, b) => {
+                let cand = Candidate::new((), Evaluation::new(vec![a, b]));
+                archive.insert(cand);
+            }
+            Op::Truncate(n) => archive.truncate(n as usize),
+        }
+    }
+
+    // Members must be pairwise non-dominated.
+    let m = archive.members();
+    for i in 0..m.len() {
+        for j in 0..m.len() {
+            if i == j {
+                continue;
+            }
+            let r = pareto_compare(&m[i].evaluation, &m[j].evaluation, &space);
+            assert!(
+                !matches!(r, Dominance::DominatedBy),
+                "archive member {i} dominated by {j}: {:?} vs {:?}",
+                m[i].evaluation.objectives,
+                m[j].evaluation.objectives,
+            );
+        }
+    }
+});

fuzz/fuzz_targets/pareto_compare.rs

@@ -0,0 +1,68 @@
+#![no_main]
+//! Fuzz `pareto_compare` for anti-symmetry and reflexivity.
+//!
+//! Invariants checked:
+//!   * `compare(a, b)` and `compare(b, a)` form an anti-symmetric pair
+//!     (`Dominates ↔ DominatedBy`, `Equal ↔ Equal`, `NonDominated ↔ NonDominated`).
+//!   * `compare(a, a) == Equal`.
+//!   * No panics on any combination of finite/non-finite floats.
+
+use arbitrary::Arbitrary;
+use libfuzzer_sys::fuzz_target;
+
+use heuropt::core::evaluation::Evaluation;
+use heuropt::core::objective::{Objective, ObjectiveSpace};
+use heuropt::pareto::dominance::{Dominance, pareto_compare};
+
+#[derive(Arbitrary, Debug)]
+struct Input {
+    a_objs: Vec<f64>,
+    b_objs: Vec<f64>,
+    a_violation: f64,
+    b_violation: f64,
+    minimize_mask: u8,
+}
+
+fuzz_target!(|input: Input| {
+    if input.a_objs.is_empty() || input.a_objs.len() != input.b_objs.len() {
+        return;
+    }
+    if input.a_objs.len() > 8 {
+        return;
+    }
+    let m = input.a_objs.len();
+    let space = ObjectiveSpace::new(
+        (0..m)
+            .map(|i| {
+                if (input.minimize_mask >> i) & 1 == 0 {
+                    Objective::minimize(format!("f{i}"))
+                } else {
+                    Objective::maximize(format!("f{i}"))
+                }
+            })
+            .collect(),
+    );
+
+    let a = Evaluation::constrained(input.a_objs.clone(), input.a_violation);
+    let b = Evaluation::constrained(input.b_objs.clone(), input.b_violation);
+
+    let ab = pareto_compare(&a, &b, &space);
+    let ba = pareto_compare(&b, &a, &space);
+    let aa = pareto_compare(&a, &a, &space);
+
+    // Anti-symmetry pairs.
+    let antisymmetric = matches!(
+        (ab, ba),
+        (Dominance::Dominates, Dominance::DominatedBy)
+            | (Dominance::DominatedBy, Dominance::Dominates)
+            | (Dominance::Equal, Dominance::Equal)
+            | (Dominance::NonDominated, Dominance::NonDominated),
+    );
+    assert!(antisymmetric, "asymmetric: ab={ab:?}, ba={ba:?}");
+
+    // Reflexivity (when objectives are finite — NaNs make equality
+    // ill-defined, so skip the check there).
+    if input.a_objs.iter().all(|v| v.is_finite()) && input.a_violation.is_finite() {
+        assert_eq!(aa, Dominance::Equal);
+    }
+});

fuzz/fuzz_targets/sbx_polymut.rs

@@ -0,0 +1,99 @@
+#![no_main]
+//! Fuzz SBX + PolynomialMutation: in-bounds parents must produce in-bounds
+//! children for any seed and any (η, per-variable-probability) pair.
+
+use arbitrary::Arbitrary;
+use libfuzzer_sys::fuzz_target;
+
+use heuropt::core::rng::rng_from_seed;
+use heuropt::prelude::*;
+
+#[derive(Arbitrary, Debug)]
+struct Input {
+    bounds: Vec<(f64, f64)>,
+    eta_sbx: f64,
+    eta_pm: f64,
+    pvp_sbx: f64,
+    pvp_pm: f64,
+    a_frac: Vec<f64>,
+    b_frac: Vec<f64>,
+    seed: u64,
+}
+
+fuzz_target!(|input: Input| {
+    let n = input.bounds.len();
+    if n == 0 || n > 8 {
+        return;
+    }
+    if !(input.eta_sbx.is_finite() && input.eta_pm.is_finite()) {
+        return;
+    }
+    if !(input.eta_sbx >= 1.0 && input.eta_sbx <= 100.0) {
+        return;
+    }
+    if !(input.eta_pm >= 1.0 && input.eta_pm <= 100.0) {
+        return;
+    }
+    let pvp_sbx = match input.pvp_sbx {
+        v if v.is_finite() && (0.0..=1.0).contains(&v) => v,
+        _ => return,
+    };
+    let pvp_pm = match input.pvp_pm {
+        v if v.is_finite() && (0.0..=1.0).contains(&v) => v,
+        _ => return,
+    };
+    // Sanitize bounds: lo < hi, finite.
+    let bounds: Vec<(f64, f64)> = input
+        .bounds
+        .iter()
+        .filter_map(|&(lo, hi)| {
+            if lo.is_finite() && hi.is_finite() && hi - lo > 1e-9 {
+                Some((lo, hi))
+            } else {
+                None
+            }
+        })
+        .collect();
+    if bounds.len() != n {
+        return;
+    }
+    if input.a_frac.len() < n || input.b_frac.len() < n {
+        return;
+    }
+
+    let p1: Vec<f64> = bounds
+        .iter()
+        .zip(&input.a_frac)
+        .map(|(&(lo, hi), &f)| {
+            let frac = if f.is_finite() { f.fract().abs() } else { 0.5 };
+            lo + frac * (hi - lo)
+        })
+        .collect();
+    let p2: Vec<f64> = bounds
+        .iter()
+        .zip(&input.b_frac)
+        .map(|(&(lo, hi), &f)| {
+            let frac = if f.is_finite() { f.fract().abs() } else { 0.5 };
+            lo + frac * (hi - lo)
+        })
+        .collect();
+
+    let mut rng = rng_from_seed(input.seed);
+    let mut sbx = SimulatedBinaryCrossover::new(bounds.clone(), input.eta_sbx, pvp_sbx);
+    let kids = sbx.vary(&[p1, p2], &mut rng);
+    assert_eq!(kids.len(), 2);
+    for c in &kids {
+        for (j, &v) in c.iter().enumerate() {
+            let (lo, hi) = bounds[j];
+            assert!(v >= lo && v <= hi, "SBX child[{j}] = {v} out of [{lo}, {hi}]");
+        }
+    }
+
+    let mut pm = PolynomialMutation::new(bounds.clone(), input.eta_pm, pvp_pm);
+    let mutated = pm.vary(std::slice::from_ref(&kids[0]), &mut rng);
+    assert_eq!(mutated.len(), 1);
+    for (j, &v) in mutated[0].iter().enumerate() {
+        let (lo, hi) = bounds[j];
+        assert!(v >= lo && v <= hi, "PM child[{j}] = {v} out of [{lo}, {hi}]");
+    }
+});

fuzz/fuzz_targets/spacing.rs

@@ -0,0 +1,42 @@
+#![no_main]
+//! Fuzz the `spacing` metric for non-negativity.
+
+use arbitrary::Arbitrary;
+use libfuzzer_sys::fuzz_target;
+
+use heuropt::core::candidate::Candidate;
+use heuropt::core::evaluation::Evaluation;
+use heuropt::core::objective::{Objective, ObjectiveSpace};
+use heuropt::metrics::spacing::spacing;
+
+#[derive(Arbitrary, Debug)]
+struct Input {
+    points: Vec<(f64, f64)>,
+}
+
+fuzz_target!(|input: Input| {
+    if input.points.len() > 64 {
+        return;
+    }
+    // Bound magnitudes so distance computations don't overflow to
+    // inf-inf=NaN — `spacing` is documented to operate on values produced
+    // by `as_minimization` of problem evaluations, not arbitrary f64s.
+    if input
+        .points
+        .iter()
+        .any(|&(a, b)| !a.is_finite() || !b.is_finite() || a.abs() > 1e150 || b.abs() > 1e150)
+    {
+        return;
+    }
+    let space = ObjectiveSpace::new(vec![Objective::minimize("f1"), Objective::minimize("f2")]);
+    let pop: Vec<Candidate<()>> = input
+        .points
+        .iter()
+        .map(|&(a, b)| Candidate::new((), Evaluation::new(vec![a, b])))
+        .collect();
+    let s = spacing(&pop, &space);
+    // Spacing can overflow to +∞ when point coordinates straddle ±f64::MAX.
+    // Contract is non-negative + non-NaN.
+    assert!(s >= 0.0, "spacing negative: {s}");
+    assert!(!s.is_nan(), "spacing is NaN");
+});