Introduction

heuropt is a practical Rust toolkit for heuristic optimization — the art of searching for good answers when the problem is too gnarly to solve analytically.

The kinds of problems heuropt is built for:

• Single-objective: "find the parameters that minimize the loss of this model." Hyperparameter tuning. Curve fitting. Calibration.
• Multi-objective: "find the trade-off curve between cost and accuracy." Engineering design. Portfolio optimization. Fleet scheduling.
• Many-objective (4+): the same idea but with enough objectives that classical Pareto methods break down. Power-grid planning. Airfoil design. Multi-criteria recommendation.

If your problem is differentiable and convex, you don't need this crate — use a gradient solver. heuropt is for the messy problems: landscapes with lots of local minima, decisions that aren't continuous (permutations, bit vectors), or evaluations that are noisy / expensive / black-box.

Why heuropt

There are other Rust optimization crates and many more in Python (pymoo, hyperopt, optuna, DEAP). heuropt's design priorities:

1. Approachable code. No trait objects in the public API. No GATs, HRTBs, generic-RNG plumbing. A junior Rust engineer should be able to read Random Search and write a new optimizer by implementing only the Optimizer<P> trait.
2. One concrete RNG type. Seeded determinism is a property tested across the crate; identical inputs always produce identical outputs.
3. Algorithms that work. Every algorithm is benchmarked against the canonical test problems (ZDT, DTLZ, Rastrigin, Rosenbrock, Ackley) and the results are checked into examples/compare-results.md so you can see what each algorithm's strengths actually are.
4. Testing as a first-class concern. 316+ unit / integration / property tests, eight cargo-fuzz targets in CI, gungraun instruction-count benchmarks. The fuzzers find real bugs and the property tests check actual invariants.

What's in the box

heuropt v0.10 ships 33 algorithms spanning:

• Single-objective continuous: Random Search, Hill Climber, (1+1)-ES, Simulated Annealing, GA, PSO, Differential Evolution, TLBO, CMA-ES, IPOP-CMA-ES, sNES, Nelder-Mead.
• Single-objective other types: UMDA (binary), Tabu Search (any), Ant Colony (permutation).
• Multi-objective (2–3): PAES, NSGA-II, SPEA2, MOPSO, IBEA, SMS-EMOA, HypE, ε-MOEA, PESA-II, AGE-MOEA, KnEA, MOEA/D.
• Many-objective (4+): NSGA-III, RVEA, GrEA.
• Sample-efficient / multi-fidelity: Bayesian Optimization, TPE, Hyperband.

Plus the operators (SBX, PolynomialMutation, BoundedGaussianMutation, LevyMutation, BitFlipMutation, SwapMutation, ClampToBounds, ProjectToSimplex), the metrics (hypervolume, spacing), and the Pareto utilities (dominance, fronts, crowding distance, Das–Dennis reference points, the ParetoArchive) that you'd expect.

Async evaluation (since v0.8, behind the async feature flag): when your evaluate function is IO-bound — calling an HTTP service, an RPC, or a subprocess — implement AsyncProblem and use run_async(&problem, concurrency).await on any algorithm in the catalog. heuropt is the only mainstream optimization library with first-class async support across its entire algorithm set.

How to use this guide

If you're new to heuropt, read it linearly:

1. Five-minute walkthrough — install, define a problem, run an optimizer, look at the result.
2. Defining a problem — the Problem trait in depth: single- vs multi-objective, constraints, custom decision types.
3. Choosing an algorithm — the decision tree, expanded with the reasoning behind each branch.

If you're already up and running, jump into the cookbook for recipes, or comparison for how heuropt stacks up against other libraries.

Five-minute walkthrough

The shortest path from a fresh project to a working optimizer.

1. Add heuropt to your Cargo.toml

[dependencies]
heuropt = "0.10"

The default feature set is small. Optional features:

• parallel — rayon-backed parallel population evaluation.
• serde — Serialize / Deserialize derives on the core data types, plus the heuropt::explorer JSON export module for the heuropt-explorer webapp.
• async — AsyncProblem trait + per-algorithm run_async for IO-bound evaluations.

heuropt = { version = "0.10", features = ["parallel"] }

2. Define a problem and run an optimizer

A problem is a struct that implements the Problem trait. You tell heuropt what kind of decision your problem takes (Vec<f64>, Vec<bool>, …), what objectives it has (minimize or maximize), and how to score one decision.

We'll fit a straight line to a handful of (x, y) data points by finding the slope and intercept that minimize the sum of squared errors — same objective as least-squares regression. For a smooth single-objective continuous problem like this, CMA-ES is a strong default.

use heuropt::prelude::*;

struct LineFit {
    points: Vec<(f64, f64)>,
}

impl Problem for LineFit {
    type Decision = Vec<f64>; // [slope, intercept]

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("sum_squared_error")])
    }

    fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
        let (slope, intercept) = (x[0], x[1]);
        let sse: f64 = self
            .points
            .iter()
            .map(|(px, py)| (py - (slope * px + intercept)).powi(2))
            .sum();
        Evaluation::new(vec![sse])
    }
}

fn main() {
    // Five noisy points roughly on the line y = 2x + 1.
    let problem = LineFit {
        points: vec![(0.0, 1.1), (1.0, 2.9), (2.0, 5.1), (3.0, 6.8), (4.0, 9.2)],
    };

    // Search box: slope and intercept each in [-10, 10].
    let bounds = RealBounds::new(vec![(-10.0, 10.0); 2]);

    let mut opt = CmaEs::new(
        CmaEsConfig {
            population_size: 12,
            generations: 80,
            initial_sigma: 1.0,
            eigen_decomposition_period: 1,
            initial_mean: None,
            seed: 42,
        },
        bounds,
    );

    let result = opt.run(&problem);
    let best = result.best.expect("at least one feasible candidate");
    let (slope, intercept) = (best.decision[0], best.decision[1]);
    println!(
        "best fit: y = {:.4} x + {:.4}   (sse = {:.4e}, evaluations = {})",
        slope, intercept, best.evaluation.objectives[0], result.evaluations,
    );

    println!();
    println!("predictions vs actual:");
    for (px, py) in &problem.points {
        let pred = slope * px + intercept;
        println!(
            "  x = {:.1}   actual = {:.2}   predicted = {:.4}   residual = {:+.4}",
            px, py, pred, py - pred,
        );
    }
}

Run with cargo run --release — heuristic optimization is allergic to debug builds. The actual output:

best fit: y = 2.0100 x + 1.0000   (sse = 1.0700e-1, evaluations = 960)

predictions vs actual:
  x = 0.0   actual = 1.10   predicted = 1.0000   residual = +0.1000
  x = 1.0   actual = 2.90   predicted = 3.0100   residual = -0.1100
  x = 2.0   actual = 5.10   predicted = 5.0200   residual = +0.0800
  x = 3.0   actual = 6.80   predicted = 7.0300   residual = -0.2300
  x = 4.0   actual = 9.20   predicted = 9.0400   residual = +0.1600

Reading the result

CMA-ES recovered slope ≈ 2.01, intercept ≈ 1.00 — within hundredths of the underlying line y = 2x + 1 that the data was sampled from. The residuals are evenly distributed in sign (3 positive, 2 negative) and small in magnitude (the largest is 0.23 at x = 3), which means the fit is balancing the noise rather than chasing any single point.

The total sum of squared errors is 0.107 — that is the value the optimizer was actually minimizing, and it matches the answer you'd get from running numpy.polyfit or solving the normal equations directly. CMA-ES is overkill for a two-parameter problem (closed-form least-squares does it in one step), but the same code shape scales straight up to nonlinear models, robust loss functions, or constrained variants where there is no closed form.

It used 960 evaluations to get there. That's population_size × generations = 12 × 80 = 960, and CMA-ES converges to machine epsilon on problems this clean in well under that budget.

4. What just happened

• Problem is the what you're optimizing.
• CMA-ES (or any other optimizer) is the how.
• CmaEsConfig is a plain public-field struct: there are no builders, no chained setters, just public fields you set directly.
• Optimizer::run returns an OptimizationResult containing the full final population, the pareto_front (just the best for single-objective), the best candidate, the total evaluations, and the number of generations.

5. Where to go next

• Multi-objective: see Defining a problem for how to express two or more objectives, and Choosing an algorithm for which optimizer fits.
• Want to know which algorithm to pick: read the README's decision tree, or jump straight to the choosing-an-algorithm chapter for the long form.
• Production patterns: the cookbook has recipes for parallelism, expensive evaluations, comparing algorithms, and more.

Defining a problem

Everything in heuropt starts with the Problem trait. This chapter walks through every shape it can take.

The trait

pub trait Problem {
    type Decision: Clone;
    fn objectives(&self) -> ObjectiveSpace;
    fn evaluate(&self, decision: &Self::Decision) -> Evaluation;
}

Three things you decide:

1. Decision — the type of the thing you're optimizing. Vec<f64> is by far the most common; Vec<bool> for binary search, Vec<usize> for permutations, your own struct for anything else.
2. objectives — how many objectives you have, what they're called, and whether each is minimized or maximized. Returned as an ObjectiveSpace.
3. evaluate — given one decision, score it. Returns an Evaluation with a vector of objective values (and optionally a constraint-violation scalar).

evaluate takes &self, so caches and lookup tables are easy. It is called many thousands of times during a typical run, so keep it fast.

Single-objective continuous

The Rosenbrock banana — minimize a smooth non-convex valley.

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

struct Rosenbrock;

impl Problem for Rosenbrock {
    type Decision = Vec<f64>;

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("f")])
    }

    fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
        let f: f64 = x.windows(2)
            .map(|w| 100.0 * (w[1] - w[0].powi(2)).powi(2) + (1.0 - w[0]).powi(2))
            .sum();
        Evaluation::new(vec![f])
    }
}
}

Multi-objective

ZDT1 — two objectives that conflict. The Pareto front is the set of non-dominated trade-offs.

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

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 n = x.len() as f64;
        let f1 = x[0];
        let g = 1.0 + 9.0 * x[1..].iter().sum::<f64>() / (n - 1.0);
        let h = 1.0 - (f1 / g).sqrt();
        let f2 = g * h;
        Evaluation::new(vec![f1, f2])
    }
}
}

For multi-objective problems, pick a Pareto-aware optimizer: NSGA-II is the canonical default; MOPSO often wins on smooth-front 2-objective problems; IBEA often wins on disconnected fronts. See choosing-an-algorithm.

Maximizing instead of minimizing

heuropt's internals normalize everything to minimization, but you declare your objective with the orientation that's natural for your problem. A scoring problem might want to maximize:

#![allow(unused)]
fn main() {
use heuropt::prelude::*;
let space = ObjectiveSpace::new(vec![
    Objective::minimize("cost"),
    Objective::maximize("accuracy"),
]);
}

Objective::maximize is a convenience for Direction::Maximize. Mix freely; the Pareto-comparison machinery handles the orientation.

Constraints

heuropt models constraints as a single non-negative scalar constraint_violation on each Evaluation. The convention:

• 0.0 (or negative) means feasible.
• Any positive value means infeasible, and bigger numbers are worse violations.

Pareto-comparison and tournament-selection helpers prefer feasible candidates and break ties on the violation magnitude — so the rule "feasibility comes first" is enforced automatically.

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

struct Constrained;
impl Problem for Constrained {
    type Decision = Vec<f64>;

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("f")])
    }

    fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
        let f: f64 = x.iter().map(|v| v * v).sum();

        // Constraint: x[0] + x[1] >= 1. Violation = how much we miss it by.
        let g1 = (1.0 - (x[0] + x[1])).max(0.0);
        let total_violation: f64 = g1; // sum of max(0, gᵢ) for each constraint

        Evaluation::constrained(vec![f], total_violation)
    }
}
}

If your constraints are very tight and the search keeps hitting them, see Constrain your search with Repair.

Decision types beyond Vec<f64>

Binary (Vec<bool>)

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

struct OneMax { bits: usize }
impl Problem for OneMax {
    type Decision = Vec<bool>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::maximize("ones")])
    }
    fn evaluate(&self, x: &Vec<bool>) -> Evaluation {
        Evaluation::new(vec![x.iter().filter(|b| **b).count() as f64])
    }
}
}

For Vec<bool> problems, UMDA is a parameter-free EDA; GA with BitFlipMutation is the GA route.

Permutations (Vec<usize>)

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

struct Tsp { distances: Vec<Vec<f64>> }
impl Problem for Tsp {
    type Decision = Vec<usize>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("length")])
    }
    fn evaluate(&self, tour: &Vec<usize>) -> Evaluation {
        let mut len = 0.0;
        for w in tour.windows(2) {
            len += self.distances[w[0]][w[1]];
        }
        len += self.distances[*tour.last().unwrap()][tour[0]];
        Evaluation::new(vec![len])
    }
}
}

For permutations, Ant Colony specializes on TSP-style problems; Tabu Search takes a user-supplied neighbor function for arbitrary discrete neighborhoods; Simulated Annealing with SwapMutation is the simplest baseline.

Custom decision types

Any Clone type works. If you have a struct, just implement Clone and you can use it. You'll need to write your own Variation impl to mutate it; see Write your own algorithm.

What Evaluation carries

pub struct Evaluation {
    pub objectives: Vec<f64>,        // one entry per objective
    pub constraint_violation: f64,    // 0.0 = feasible
}

That's it. Construct with Evaluation::new for unconstrained problems or Evaluation::constrained when you have a violation.

Summary

• Implement Problem with your decision type.
• Declare objectives via ObjectiveSpace (mix minimize/maximize freely).
• Return an Evaluation from evaluate.
• For constraints, set constraint_violation > 0 for infeasible decisions; heuropt's selection helpers prefer feasibles automatically.

Next: Choosing an algorithm walks through the decision tree.

Choosing an algorithm

The README has a compact decision tree. This chapter expands it with the reasoning behind each branch.

Step 0: How expensive is one evaluation?

This is the first fork because it changes everything that comes after it.

For the cheap-eval branch, you have the run of the catalog. For the expensive branch, classical evolutionary methods waste your evaluation budget — go to Bayesian Optimization or TPE. For the very expensive branch where each eval has a tunable budget (epochs, MC samples, sim steps), Hyperband over the PartialProblem trait is the move.

Step 1: How many objectives?

The biggest fork.

• One — there's a single best answer. Pick from the single-objective branch.
• Two or three — a Pareto front. Pick from the multi-objective branch.
• Four or more — a many-objective Pareto front; classical multi-objective methods break down here because almost every pair of points is non-dominated. Pick from the many-objective branch.

Pareto front: the set of decisions where you cannot improve any objective without sacrificing another. In a 2-objective minimize problem, plot every solution; the Pareto front is the lower-left envelope.

If you found yourself staring at a single composite score that's a weighted sum of conflicting goals, you probably have a multi-objective problem in disguise. A weighted sum bakes in your preferences before you've seen the trade-off; running a multi-objective optimizer first and picking off the front later is almost always a better workflow (see Pick one answer off a Pareto front).

Step 2 — single-objective continuous

These all take Vec<f64> decisions.

Smooth, low-to-moderate dimension

CMA-ES is the strong default. It adapts the search distribution's covariance to the local landscape. On the comparison harness it hits machine epsilon on Rosenbrock at 30 000 evaluations.

For very low-dimensional smooth problems (≤ 5 dim), Nelder-Mead is deterministic and converges to f = 0 exactly on Rosenbrock.

High dimension, smooth

sNES uses a diagonal covariance — cheaper per step than CMA-ES at the cost of being unable to model rotated landscapes. Worth trying when CMA-ES's O(d²) per-step cost hurts.

Multimodal landscapes

Multimodal = many local minima that aren't the global one. Rastrigin and Ackley are classic traps.

IPOP-CMA-ES is CMA-ES with an increasing-population restart strategy specifically designed for this. On the harness it drops vanilla CMA-ES's Rastrigin score from f = 2.35 to f = 0.13.

Differential Evolution is rarely beaten on cheap multimodal continuous problems. On Rastrigin it ties with (1+1)-ES at f = 0.

Simulated Annealing is a cheap, generic baseline that escapes local optima via temperature decay.

Want parameter-free

TLBO (Teaching-Learning-Based Optimization) has no F, CR, w, or σ to tune. Often a respectable middle-of-the-pack performer.

Smallest possible self-adapting baseline

(1+1)-ES — Rechenberg's 1973 (1+1)-ES with the one-fifth success rule. On the harness it hits f = 0 on Rastrigin in 50 000 evaluations.

Just want a baseline

Random Search. Useful as a sanity check: if your fancy optimizer can't beat random search, something is wrong (with the fancy optimizer or with the problem).

Step 2 — single-objective other types

heuropt's permutation operator toolkit covers four crossovers (OrderCrossover, PartiallyMappedCrossover, CycleCrossover, EdgeRecombinationCrossover) and four mutations (SwapMutation, InversionMutation, InsertionMutation, ScrambleMutation), plus two initializers for strict and multiset permutations. See Optimize a permutation for the full picker.

Step 2 — multi-objective (2 or 3)

Strong default

MOEA/D is the most consistent performer on the harness. It decomposes the problem into many scalar sub-problems (Tchebycheff or weighted sum) and solves them in parallel — fast per generation, and robust: it finishes top-3 on every multi- and many-objective table (convex, disconnected, spherical and linear fronts; 2 through 10 objectives) and is consistently the fastest or near-fastest. It rarely wins a table outright — a specialist usually does — but it never lands badly. One caveat from the literature: MOEA/D's spread depends on the weight-vector distribution and the scalarizing function, so it can leave gaps on highly irregular or degenerate fronts; the DTLZ/ZDT suite here doesn't stress that.

NSGA-II is the other safe default — the canonical Pareto-based EA: fast, well-understood, diversity-preserving via crowding distance. On the harness it's edged out by MOEA/D on every multi-objective table and degrades past ~4 objectives (see the many-objective section), but it stays a solid 2–3-objective pick and is the established choice for combinatorial encodings: drop in ShuffledPermutation + a permutation crossover and it solves bi-objective TSP; drop in a binary initializer and BitFlipMutation and it solves bi-objective knapsack. See Multi-objective combinatorial problems.

Real-valued, smooth front, want best convergence

MOPSO (multi-objective PSO with archive). On ZDT1 it wins hypervolume outright and converges 100× tighter than the dominance-based methods.

Better front quality than the default

IBEA (indicator-based) is consistently the best of the dominance-based methods on the harness — wins ZDT3 hypervolume and DTLZ2 mean distance by 24×. It uses an additive ε-indicator for selection rather than dominance + crowding.

SPEA2 (strength + density) — solid alternative; explicit external archive separate from the population.

SMS-EMOA uses exact hypervolume contribution for selection. Elegant in theory; in practice on the harness budgets here it underperforms NSGA-II. Worth the higher per-step cost only when exact HV contribution is the right discriminator.

Disconnected or non-convex front

A disconnected front (separate arcs, like ZDT3) and a non-convex but contiguous front are different problems — don't conflate them.

For a disconnected front, IBEA is the clear pick: on the harness it wins ZDT3 — the disconnected-front benchmark — outright on hypervolume, with MOEA/D and NSGA-II close behind. Counter-intuitively the geometry-aware methods below trail here: estimating a single front geometry or chasing knee points doesn't help when the front is in pieces (on ZDT3, AGE-MOEA and KnEA finish last).

For a non-convex but contiguous front:

AGE-MOEA estimates the front geometry adaptively (the L_p parameter p is fit from data each generation).

KnEA favors knee points — the regions of the front where small gains in one objective cost large losses in another.

Region-based diversity

PESA-II uses grid hyperboxes to drive selection — divide the objective space into a grid, pick from the least-crowded boxes.

ε-MOEA uses an ε-grid archive that auto-limits its size.

Just one starting decision (no population budget)

PAES — (1+1)-ES with a Pareto archive. Cheap, simple, useful when your evaluations are expensive enough that you can't afford a population.

Step 2 — many-objective (4+)

Strong default

MOEA/D again. Decomposition sidesteps the dominance resistance that breaks Pareto-based methods at high objective count — each scalar sub-problem still has a clear best, even when almost every pair of solutions is mutually non-dominated. On the harness it is #2 on every many-objective table (DTLZ2 at 4 and 10 objectives, DTLZ1 at 8), and fast every time. NSGA-II is the cautionary tale: on DTLZ2 at 10 objectives it finishes last — behind random search — because its crowding distance has no dominance signal left to refine.

Linear / simplex-shaped front (e.g., DTLZ1)

GrEA — grid coords drive ranking. On 3-objective DTLZ1 it beats NSGA-III by 3× and AGE-MOEA by 2.5×, and it wins the 8-objective DTLZ1 table outright.

MOEA/D — also #2 on both DTLZ1 tables.

Curved / unknown front geometry

NSGA-III — reference-point niching; the canonical many-objective method by reputation, though on the harness MOEA/D outperforms it on every table. Reach for it when you specifically want reference-point niching.

AGE-MOEA — estimates L_p geometry per generation.

RVEA — reference vectors with adaptive penalty.

Indicator-based selection

IBEA — additive ε-indicator; doesn't degrade at high obj count.

HypE — Monte Carlo hypervolume estimation; scales to arbitrary objective count where exact HV is too expensive.

Step 3: Are there hard constraints?

heuropt models constraints as a single scalar constraint_violation on each Evaluation. Three escalations when the feasibility region is hard to find:

1. Penalty-only. Just set constraint_violation > 0 for infeasible decisions. The default tournament/Pareto comparisons prefer feasibles automatically.
2. Repair. Implement Repair<D> (or use the provided ClampToBounds / ProjectToSimplex) to project infeasible decisions back into the feasible region. Pair with a Variation in a CompositeVariation for bounds-aware variants.
3. Stochastic ranking. Use stochastic_ranking_select instead of tournament_select_single_objective. It probabilistically explores near-feasibility instead of strict feasibility-first ordering, which helps when feasible regions are narrow.

See Constrain your search with Repair for worked examples.

Step 4: Should you parallelize?

Enable the parallel feature flag if your evaluate takes more than ~50 µs. Population-based algorithms (Random Search, NSGA-II, Differential Evolution, SPEA2, IBEA, MOPSO, …) batch- evaluate via rayon when the feature is on. Seeded runs stay bit-identical to serial mode.

heuropt = { version = "0.10", features = ["parallel"] }

If your evaluation is IO-bound (HTTP request, RPC, subprocess) rather than CPU-bound, use the async feature instead — it gives you AsyncProblem and a run_async(&problem, concurrency).await method on every algorithm in the catalog. See the Async evaluation cookbook recipe.

TL;DR table

Cookbook

Short, focused recipes for the patterns that come up in practice. Each recipe is self-contained and small enough to copy into your own project.

Recipes

• Parallelize evaluation with rayon — when your evaluate is non-trivial CPU work, the parallel feature pays for itself almost immediately.
• Async evaluation — when your evaluate is IO-bound (HTTP / RPC / subprocess), the async feature lets the optimizer await many evaluations concurrently. The differentiating feature vs other optimization libraries.
• Tune a model with expensive evaluations — Bayesian Optimization, TPE, and Hyperband for the 50–500-eval regime.
• Compare two algorithms on your problem — multi-seed harness pattern straight from examples/compare.rs.
• Optimize a permutation (TSP-style) — the permutation operator toolkit (OX / PMX / CX / ERX + Inversion / Insertion / Scramble), plus Ant Colony for distance-matrix TSP.
• Multi-objective combinatorial problems — bi-objective TSP, bi-objective knapsack (Vec<bool>), and 3-objective JSS via NSGA-II / NSGA-III.
• Constrain your search with Repair — bounds, simplex projection, custom repair.
• Pick one answer off a Pareto front — the a-posteriori weighted-decision pattern from the jiggly_tuning example.
• Explore your results in a webapp — export an OptimizationResult to JSON and browse it interactively at heuropt-explorer — parallel coordinates, scatter, range filters, weighted ranking.
• Write your own algorithm — implement Optimizer<P> from scratch, à la the examples/custom_optimizer.rs walkthrough.

Parallelize evaluation with rayon

If a single call to your evaluate takes more than ~50 µs, enabling the parallel feature usually pays for itself immediately on population-based algorithms. Each generation evaluates an entire population, and rayon parallelizes that batch.

Enable the feature

[dependencies]
heuropt = { version = "0.10", features = ["parallel"] }

There's nothing else to opt into in your code. The population-evaluation helper is feature-gated; with parallel on it uses rayon::into_par_iter internally, with parallel off it falls back to plain into_iter.

Determinism still holds

Seeded runs are bit-identical between the serial and parallel modes. The trick is that population members are evaluated in parallel but assembled back into the same order. Variation, selection, and the RNG are all driven by the main thread, so seed-stability tests still pass.

Which algorithms benefit

Algorithms with a per-generation evaluate_batch:

• Random Search, NSGA-II, NSGA-III, SPEA2, MOEA/D, MOPSO, IBEA, SMS-EMOA, HypE, PESA-II, ε-MOEA, AGE-MOEA, KnEA, GrEA, RVEA.
• Differential Evolution and GA benefit on the initial population and offspring batches.

Steady-state algorithms (PAES, Simulated Annealing, Hill Climber, (1+1)-ES) only evaluate one or a few candidates per iteration, so the parallel feature gives them nothing — leave it off if those are your primary optimizers.

Worked example

The Sphere problem is too cheap to actually benefit from parallelism — this example just shows the shape. In real workloads evaluate is the expensive bit (a simulation, a model fit, an HTTP call).

use heuropt::prelude::*;

struct ExpensiveSphere;
impl Problem for ExpensiveSphere {
    type Decision = Vec<f64>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("f")])
    }
    fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
        // Pretend this is a 5 ms simulation.
        std::thread::sleep(std::time::Duration::from_millis(5));
        Evaluation::new(vec![x.iter().map(|v| v * v).sum::<f64>()])
    }
}

fn main() {
    let bounds = vec![(-1.0_f64, 1.0_f64); 5];
    let mut opt = DifferentialEvolution::new(
        DifferentialEvolutionConfig {
            population_size: 16,
            generations: 50,
            differential_weight: 0.5,
            crossover_probability: 0.9,
            seed: 42,
        },
        RealBounds::new(bounds),
    );
    let r = opt.run(&ExpensiveSphere);
    println!("best f = {}", r.best.unwrap().evaluation.objectives[0]);
}

With the parallel feature on, each generation's 16 evaluations run across rayon's worker threads. On a 16-core machine the wall-clock cost per generation drops from 16 × 5 ms = 80 ms to roughly 5 ms + scheduling overhead.

Sizing your thread pool

heuropt uses rayon's global thread pool. Override the size with:

rayon::ThreadPoolBuilder::new().num_threads(8).build_global().unwrap();

Run this before any heuropt call, or use rayon's install API to scope it.

When parallelism doesn't help

• Your evaluate is sub-microsecond (Sphere, Rastrigin, Ackley unweighted) — the rayon scheduling overhead exceeds the work.
• You're already running multiple seeds in parallel at the harness level (see Compare two algorithms). Stacking parallelism rarely helps.
• The algorithm is steady-state (PAES, SA, hill climber).

parallel vs async

Both can be on at once if your evaluation does both substantial CPU work and IO. The two features are independent.

Async evaluation

When your evaluate does IO — calls an HTTP service, sends an RPC, spawns a subprocess — await-ing it from the optimizer is much more efficient than blocking a thread per evaluation. heuropt ships first-class async support behind the async feature flag.

This is the differentiating capability vs pymoo / hyperopt / optuna / DEAP / MOEA Framework — none of those have a native async evaluation path.

Enable the feature

[dependencies]
heuropt = { version = "0.10", features = ["async"] }

# Pick whatever async runtime you want; heuropt itself depends only on
# `futures`. The example below uses tokio.
tokio = { version = "1", features = ["rt-multi-thread", "macros", "time"] }

Implement AsyncProblem

It mirrors the regular Problem trait one-for-one — same Decision type, same objectives(), but evaluate is replaced with evaluate_async returning a future.

#![allow(unused)]
fn main() {
use heuropt::core::async_problem::AsyncProblem;
use heuropt::prelude::*;

struct RemoteService;

impl AsyncProblem for RemoteService {
    type Decision = Vec<f64>;

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("loss")])
    }

    async fn evaluate_async(&self, x: &Vec<f64>) -> Evaluation {
        // Real workload: HTTP call to a model-scoring service, an RPC,
        // a subprocess. Here we just sleep to model 20 ms latency.
        tokio::time::sleep(std::time::Duration::from_millis(20)).await;
        let loss: f64 = x.iter().map(|v| v * v).sum();
        Evaluation::new(vec![loss])
    }
}
}

Run the optimizer with run_async

run_async(&problem, concurrency).await is provided by every algorithm in the catalog as of v0.8. concurrency caps how many evaluations are in-flight at once.

use heuropt::core::async_problem::AsyncProblem;
use heuropt::prelude::*;
struct RemoteService;
impl AsyncProblem for RemoteService {
    type Decision = Vec<f64>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("loss")])
    }
    async fn evaluate_async(&self, x: &Vec<f64>) -> Evaluation {
        Evaluation::new(vec![x.iter().map(|v| v * v).sum::<f64>()])
    }
}
#[tokio::main]
async fn main() {
    let bounds = vec![(-1.0_f64, 1.0_f64); 4];
    let mut opt = DifferentialEvolution::new(
        DifferentialEvolutionConfig {
            population_size: 16,
            generations: 50,
            differential_weight: 0.5,
            crossover_probability: 0.9,
            seed: 42,
        },
        RealBounds::new(bounds),
    );
    let r = opt.run_async(&RemoteService, /* concurrency */ 8).await;
    println!("best: {}", r.best.unwrap().evaluation.objectives[0]);
}

Picking concurrency

Concurrency is the maximum in-flight evaluation count. Tradeoffs:

The bigger you go, the more memory the in-flight futures hold and the more load you put on the downstream service. A reasonable starting point is min(pop_size, 16) and increase only if the downstream service is comfortable.

Determinism

Same seed produces the same final result whether you use run or run_async, provided your async evaluate_async is itself deterministic. heuropt drives the RNG and selection on the main task; only the evaluations are concurrent, and the evaluate_batch_async helper preserves input order before feeding results back to the algorithm.

What the worked example shows

examples/async_eval.rs runs Random Search (200 evaluations × 20 ms each) at concurrency = 1, 4, 16 and Differential Evolution at concurrency = 8. On a recent machine:

RandomSearch with 200 evaluations (20 ms each)

concurrency =  1  elapsed ≈ 4250 ms     (sequential 200 × 20 ms)
concurrency =  4  elapsed ≈ 2100 ms     (2× speedup, batch_size=2 caps it)
concurrency = 16  elapsed ≈ 2100 ms     (same — batch_size dominates)

DifferentialEvolution at concurrency=8
elapsed ≈ 230 ms     (8 ants run in parallel each generation)

Run it yourself: cargo run --release --features async --example async_eval.

Which algorithms support run_async?

All 33 algorithms in the catalog. The shape of the async path depends on the algorithm:

• Population-based / batch-evaluating — NSGA-II, NSGA-III, SPEA2, MOEA/D, IBEA, SMS-EMOA, HypE, ε-MOEA, PESA-II, AGE-MOEA, KnEA, GrEA, RVEA, MOPSO, GA, DE, PSO, CMA-ES, IPOP-CMA-ES, sNES, TLBO, UMDA, Ant Colony, Random Search. Each generation's offspring evaluations are fanned out concurrently up to concurrency.
• Steady-state (one-eval-per-step) — Hill Climber, Simulated Annealing, (1+1)-ES, PAES, Nelder-Mead. The concurrency parameter is accepted for API uniformity but evaluation order is inherently sequential.
• Tabu Search — fans out the K-neighbor batch each step.
• Surrogate (BO, TPE) — fans out the initial design batch, then awaits per-iteration acquisitions sequentially (the surrogate must update before the next point is chosen).
• Hyperband — uses the separate AsyncPartialProblem trait (multi-fidelity); each Successive-Halving rung's evaluations fan out concurrently.

Async vs parallel

Both can be on at once if your evaluation does both substantial CPU work and IO. The two features are independent.

Tune a model with expensive evaluations

Population-based EAs throw thousands of evaluations at a problem. If each evaluation costs a minute (a model training run, a CFD solve, a real-world measurement) you can't afford that. heuropt has three algorithms aimed at this regime.

When each is right

• Black-box, fixed cost per eval, smooth-ish landscape → BO.
• Black-box, fixed cost per eval, no time to tune the surrogate → TPE.
• Each eval has a tunable fidelity → Hyperband.

Bayesian Optimization

A worked example with a synthetic 5-D problem and a 60-evaluation budget — same configuration the compare harness uses.

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

struct Rosenbrock5D;
impl Problem for Rosenbrock5D {
    type Decision = Vec<f64>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("f")])
    }
    fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
        let f: f64 = x.windows(2).map(|w|
            100.0 * (w[1] - w[0].powi(2)).powi(2) + (1.0 - w[0]).powi(2)
        ).sum();
        Evaluation::new(vec![f])
    }
}

let bounds = vec![(-2.048_f64, 2.048_f64); 5];
let mut opt = BayesianOpt::new(
    BayesianOptConfig {
        evaluations: 60,
        initial_samples: 10,
        length_scale: 1.0,
        signal_variance: 1.0,
        noise_variance: 1e-6,
        seed: 42,
    },
    RealBounds::new(bounds),
);
let r = opt.run(&Rosenbrock5D);
println!("best f after 60 evals: {}", r.best.unwrap().evaluation.objectives[0]);
}

Honest disclosure. On the comparison harness this default configuration produces f ≈ 3170 ± 2920 on Rosenbrock 5-D — well below what a tuned BO can do. The default RBF kernel without per-problem hyperparameter tuning is the limitation. For real workloads, consider:

• More evaluations (200+ instead of 60).
• Tuning length_scale to a known scale of your problem (lower for high-frequency landscapes, higher for smooth ones).
• TPE instead of BO if you don't want to tune the kernel.

Tree-structured Parzen Estimator

TPE keeps two density estimates — l(x) over historical good points and g(x) over the rest — and picks new candidates that maximize the ratio. Cheaper per step than a GP and famously robust without hand-tuning.

#![allow(unused)]
fn main() {
use heuropt::prelude::*;
struct Rosenbrock5D;
impl Problem for Rosenbrock5D {
    type Decision = Vec<f64>;
    fn objectives(&self) -> ObjectiveSpace { ObjectiveSpace::new(vec![Objective::minimize("f")]) }
    fn evaluate(&self, _x: &Vec<f64>) -> Evaluation { Evaluation::new(vec![0.0]) }
}
let bounds = vec![(-2.048_f64, 2.048_f64); 5];
let mut opt = Tpe::new(
    TpeConfig {
        evaluations: 60,
        initial_samples: 10,
        gamma: 0.25,
        candidates_per_step: 24,
        bandwidth_factor: 1.06,
        seed: 42,
    },
    RealBounds::new(bounds),
);
let _r = opt.run(&Rosenbrock5D);
}

gamma is the fraction of best points used as l(x); 0.25 is the canonical Bergstra value.

Hyperband

Hyperband needs your problem to implement PartialProblem — that is, you can evaluate at a tunable fidelity (e.g. number of training epochs). The algorithm schedules many cheap-fidelity runs and promotes only the survivors to higher fidelity.

#![allow(unused)]
fn main() {
use heuropt::prelude::*;
use heuropt::core::partial_problem::PartialProblem;

struct ModelTuning;
impl Problem for ModelTuning {
    type Decision = Vec<f64>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("val_loss")])
    }
    fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
        // Full-fidelity eval = train at max_epochs.
        self.evaluate_at_budget(x, 100.0)
    }
}
impl PartialProblem for ModelTuning {
    fn evaluate_at_budget(&self, x: &Vec<f64>, budget: f64) -> Evaluation {
        // Replace with: train your model for `budget` epochs, return val_loss.
        // For demo, pretend more budget = lower noisy loss.
        let lr = x[0];
        let wd = x[1];
        let loss = (lr - 0.001).powi(2) + (wd - 1e-4).powi(2)
                 + 1.0 / (budget + 1.0);
        Evaluation::new(vec![loss])
    }
}

let bounds = vec![(1e-5_f64, 1e-1), (1e-6_f64, 1e-2)];
let mut hyperband = Hyperband::new(
    HyperbandConfig {
        max_budget: 100.0,
        eta: 3.0,
        seed: 42,
    },
    RealBounds::new(bounds),
);
let _r = hyperband.run(&ModelTuning);
}

max_budget is the most epochs (or whatever your fidelity unit is) you'd ever spend on a single config. eta controls how aggressive the elimination is — 3.0 is the classic value; higher means more aggressive culling.

Strategy: combining surrogate + multi-fidelity

The state of the art (BOHB) combines BO with Hyperband: TPE picks the configurations Hyperband then evaluates at increasing fidelity. heuropt doesn't ship a unified BOHB but the building blocks are there — wrap your PartialProblem with a TPE-driven sampler and feed the picks into Hyperband. PRs welcome.

Compare two algorithms on your problem

The harness in examples/compare.rs runs every applicable algorithm against every test problem with N seeds and reports mean ± std. You can lift the same pattern for your own problem in ~30 lines.

The pattern

1. Wrap your problem in a struct that implements Problem.
2. Pick a few candidate algorithms.
3. For each algorithm × seed, run and record the metric you care about.
4. Print mean ± std.

Worked example

use heuropt::prelude::*;
use std::time::Instant;

struct MyProblem;
impl Problem for MyProblem {
    type Decision = Vec<f64>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("f")])
    }
    fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
        // your problem here
        Evaluation::new(vec![x.iter().map(|v| v * v).sum::<f64>()])
    }
}

const SEEDS: u64 = 10;
const DIM: usize = 5;
const BUDGET: usize = 30_000;

fn main() {
    let bounds: Vec<(f64, f64)> = vec![(-5.0, 5.0); DIM];

    let mut best_de = vec![];
    let mut best_cmaes = vec![];
    let mut best_ipop = vec![];
    let mut t_de = vec![];
    let mut t_cmaes = vec![];
    let mut t_ipop = vec![];

    for seed in 0..SEEDS {
        // Differential Evolution
        let t = Instant::now();
        let mut de = DifferentialEvolution::new(
            DifferentialEvolutionConfig {
                population_size: 30,
                generations: BUDGET / 30,
                differential_weight: 0.5,
                crossover_probability: 0.9,
                seed,
            },
            RealBounds::new(bounds.clone()),
        );
        let r = de.run(&MyProblem);
        t_de.push(t.elapsed().as_millis() as f64);
        best_de.push(r.best.unwrap().evaluation.objectives[0]);

        // CMA-ES
        let t = Instant::now();
        let mut cma = CmaEs::new(
            CmaEsConfig {
                population_size: 12,
                generations: BUDGET / 12,
                initial_sigma: 1.0,
                eigen_decomposition_period: 1,
                initial_mean: None,
                seed,
            },
            RealBounds::new(bounds.clone()),
        );
        let r = cma.run(&MyProblem);
        t_cmaes.push(t.elapsed().as_millis() as f64);
        best_cmaes.push(r.best.unwrap().evaluation.objectives[0]);

        // IPOP-CMA-ES
        let t = Instant::now();
        let mut ipop = IpopCmaEs::new(
            IpopCmaEsConfig {
                base: CmaEsConfig {
                    population_size: 12,
                    generations: BUDGET / 12 / 4,
                    initial_sigma: 1.0,
                    eigen_decomposition_period: 1,
                    initial_mean: None,
                    seed,
                },
                max_restarts: 3,
                population_factor: 2.0,
                seed,
            },
            RealBounds::new(bounds.clone()),
        );
        let r = ipop.run(&MyProblem);
        t_ipop.push(t.elapsed().as_millis() as f64);
        best_ipop.push(r.best.unwrap().evaluation.objectives[0]);
    }

    println!("{:<12} {:>14} {:>10}", "algorithm", "best f (mean±std)", "ms");
    print_row("DE",          &best_de,    &t_de);
    print_row("CMA-ES",      &best_cmaes, &t_cmaes);
    print_row("IPOP-CMA-ES", &best_ipop,  &t_ipop);
}

fn print_row(name: &str, values: &[f64], times: &[f64]) {
    let (m, s) = mean_std(values);
    let (t, _) = mean_std(times);
    println!("{:<12} {:>10.3e} ± {:>5.2e} {:>6.0}", name, m, s, t);
}

fn mean_std(xs: &[f64]) -> (f64, f64) {
    let n = xs.len() as f64;
    let m = xs.iter().sum::<f64>() / n;
    let v = xs.iter().map(|x| (x - m).powi(2)).sum::<f64>() / n;
    (m, v.sqrt())
}

What to record

• best.evaluation.objectives[0] for single-objective.
• hypervolume_2d(&result.pareto_front, &space, ref_point) for 2-objective.
• spacing(&result.pareto_front, &space) for front uniformity.
• result.evaluations to cross-check that every algorithm got the same evaluation budget.
• Wall-clock Instant::now() deltas for runtime comparison.

Pitfalls

• Population size matters. Different algorithms have very different sweet spots. Don't just give them all the same population — the README's algorithm pages note typical defaults.
• Different algorithms count "generations" differently. What matters is the total evaluations count. Set generations = BUDGET / population_size to match across algorithms (with caveats for steady-state algorithms like SMS-EMOA that evaluate one offspring per generation).
• One seed is not a comparison. Always run ≥ 5 seeds; ≥ 10 is better. Single-seed comparisons are noise.
• The harness in examples/compare.rs is the canonical version. When in doubt, copy from there.

Optimize a permutation (TSP-style)

When your decision is "an ordering" — visiting cities, scheduling jobs, routing — the natural representation is Vec<usize>. heuropt ships three reasonable starting points:

• A purpose-built algorithm — Ant Colony for TSP-shaped problems with a distance matrix.
• A genetic algorithm with the permutation operator toolkit — the most general option, and the right choice when you want to bring your own evaluator without a pheromone metaphor.
• A trajectory method — Simulated Annealing + SwapMutation for a tiny baseline, or Tabu Search when you have a custom neighbor function.

This recipe walks through all three, with the bulk of the page on the GA toolkit, since it's the most flexible. For the multi-objective versions (bi-objective TSP, bi-objective JSS, Pareto fronts) see Multi-objective combinatorial problems.

The permutation operator toolkit

heuropt ships a complete set of permutation operators in the prelude. You compose them with CompositeVariation into a crossover-plus-mutation pipeline and feed them to any GA-shaped algorithm.

Initializers

Crossovers

All four take two parents and return two children. They assume strict permutations — applying them to multiset encodings (like JSS) will break the multiset.

For TSP specifically, ERX usually wins on Pareto-front quality at the cost of being ~70% slower per crossover. See Multi-objective combinatorial problems for a head-to-head comparison.

Mutations

All five take one parent and return one child. All four below preserve both strict permutations and multiset encodings, so they're safe for JSS too.

Quick "what should I use?" guide

Single-objective TSP with a Genetic Algorithm

This is the toolkit's headline pattern. It mirrors the examples/tsp_ulysses16.rs benchmark, which converges to the known TSPLIB optimum for the 16-city Ulysses instance.

use heuropt::prelude::*;

struct Tsp {
    distances: Vec<Vec<f64>>,
}

impl Problem for Tsp {
    type Decision = Vec<usize>;

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("tour_length")])
    }

    fn evaluate(&self, tour: &Vec<usize>) -> Evaluation {
        let n = tour.len();
        let mut len = 0.0;
        for i in 0..n {
            len += self.distances[tour[i]][tour[(i + 1) % n]];
        }
        Evaluation::new(vec![len])
    }
}

fn main() {
    // Replace with your actual distance matrix.
    let n: usize = 16;
    let distances = vec![vec![0.0_f64; n]; n];
    let problem = Tsp { distances };

    let mut optimizer = GeneticAlgorithm::new(
        GeneticAlgorithmConfig {
            population_size: 150,
            generations: 1500,
            tournament_size: 3,
            elitism: 4,
            seed: 42,
        },
        ShuffledPermutation { n },
        CompositeVariation {
            crossover: OrderCrossover,
            mutation: InversionMutation,
        },
    );
    let r = optimizer.run(&problem);
    let best = r.best.unwrap();
    println!("best tour length: {:.0}", best.evaluation.objectives[0]);
    println!("tour: {:?}", best.decision);
}

A few things to notice:

• Decision type is Vec<usize>. Every operator in the toolkit is generic over the decision type via Variation<Vec<usize>>, so the whole pipeline composes naturally.
• CompositeVariation is the wiring. It runs the crossover first, then runs the mutation on each child. For a single-parent operator pair (e.g., two mutations stacked), it still works — the "crossover" slot just becomes a first-stage mutation.
• Tournament size 3 and elitism 4 are slightly stronger than the defaults; small permutation GAs benefit from a touch more selection pressure.

Job-shop scheduling — multiset encodings

JSS problems use a different encoding: a string of length n_jobs × n_machines where each job id appears n_machines times. The k-th occurrence of job j represents the k-th operation of job j. This is a multiset permutation, not a strict permutation, and the crossovers above (OX, PMX, CX, ERX) will break it because they assume each value appears exactly once.

Use ShuffledMultisetPermutation for the initializer:

use heuropt::prelude::*;

// 6 jobs × 6 machines (FT06 layout): each job id 0..6 appears 6 times.
let initializer = ShuffledMultisetPermutation::new(vec![6; 6]);

For variation you have two options:

1. Mutation only. The four mutations above all preserve the multiset, so you can drive a GA with just InversionMutation or SwapMutation and skip crossover. This works on small JSS instances; on larger ones search becomes slow.

2. Add a JSS-aware crossover. The standard choice is POX (Precedence-preserving Order-based Crossover, Bierwirth 1996). It's not in the library because every JSS instance specifies its own number of distinct ids and POX needs that constant; defining it locally per example keeps the type clean:

   #![allow(unused)]
   fn main() {
   use heuropt::prelude::*;
   use rand::Rng as _;
   
   const N_JOBS: usize = 6;
   
   /// POX — partition job ids into two sets J1/J2; child takes positions
   /// of J1-jobs from parent A and fills the remaining positions with
   /// J2-jobs from parent B in B's order. Preserves the multiset.
   #[derive(Default)]
   struct PrecedenceOrderCrossover;
   
   impl Variation<Vec<usize>> for PrecedenceOrderCrossover {
       fn vary(&mut self, parents: &[Vec<usize>], rng: &mut Rng) -> Vec<Vec<usize>> {
           let p1 = &parents[0];
           let p2 = &parents[1];
           let mut in_j1 = [false; N_JOBS];
           loop {
               for slot in &mut in_j1 {
                   *slot = rng.random_bool(0.5);
               }
               let count = in_j1.iter().filter(|&&b| b).count();
               if count > 0 && count < N_JOBS { break; }
           }
           vec![pox_child(p1, p2, &in_j1), pox_child(p2, p1, &in_j1)]
       }
   }
   
   fn pox_child(donor: &[usize], filler: &[usize], in_donor_set: &[bool]) -> Vec<usize> {
       let n = donor.len();
       let mut child = vec![usize::MAX; n];
       for k in 0..n {
           if in_donor_set[donor[k]] {
               child[k] = donor[k];
           }
       }
       let mut idx = 0;
       for &v in filler {
           if !in_donor_set[v] {
               while idx < n && child[idx] != usize::MAX { idx += 1; }
               child[idx] = v;
               idx += 1;
           }
       }
       child
   }
   }

   See examples/jss_ft06_bi.rs and examples/mo_jss_la01.rs for the complete worked examples.

A full JSS evaluator walks the schedule string left-to-right, tracking per-job operation counters and per-machine clocks:

fn evaluate(&self, schedule: &Vec<usize>) -> Evaluation {
    let mut job_next = [0_usize; N_JOBS];
    let mut job_clock = [0.0_f64; N_JOBS];
    let mut machine_clock = [0.0_f64; N_MACHINES];
    for &job in schedule {
        let k = job_next[job];
        let m = ROUTING[job][k];
        let t = PROCESSING_TIME[job][k];
        let start = job_clock[job].max(machine_clock[m]);
        let end = start + t;
        job_clock[job] = end;
        machine_clock[m] = end;
        job_next[job] = k + 1;
    }
    let makespan = machine_clock.iter().cloned().fold(0.0_f64, f64::max);
    Evaluation::new(vec![makespan])
}

Comparing crossover operators

Tuning the right operator combo matters more than tuning population size. The pattern is: hold everything constant, swap the operator, record the metric:

use heuropt::prelude::*;
use heuropt::metrics::hypervolume_2d;

fn run_with<C: Variation<Vec<usize>>>(crossover: C) -> f64 {
    let mut opt = GeneticAlgorithm::new(
        GeneticAlgorithmConfig { /* identical config */ ..Default::default() },
        ShuffledPermutation { n: 25 },
        CompositeVariation { crossover, mutation: InversionMutation },
    );
    opt.run(&problem).best.unwrap().evaluation.objectives[0]
}

println!("OX:  {:.0}", run_with(OrderCrossover));
println!("PMX: {:.0}", run_with(PartiallyMappedCrossover));
println!("CX:  {:.0}", run_with(CycleCrossover));
println!("ERX: {:.0}", run_with(EdgeRecombinationCrossover));

For Pareto-front problems use hypervolume, not single-objective fitness, as the comparison metric — see examples/tsp_operators_compare.rs for the bi-objective version.

TSP with Ant Colony

When your problem is genuinely TSP-shaped — symmetric distance matrix, visit-every-node — Ant Colony is purpose-built and worth a look. It doesn't use crossover or mutation; instead it deposits pheromone trails that bias future ants toward good edges.

use heuropt::prelude::*;

struct Tsp {
    distances: Vec<Vec<f64>>,
}

impl Problem for Tsp {
    type Decision = Vec<usize>;

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("length")])
    }

    fn evaluate(&self, tour: &Vec<usize>) -> Evaluation {
        let mut len = 0.0;
        for w in tour.windows(2) {
            len += self.distances[w[0]][w[1]];
        }
        len += self.distances[*tour.last().unwrap()][tour[0]];
        Evaluation::new(vec![len])
    }
}

fn main() {
    let cities = vec![(0.0, 0.0), (1.0, 5.0), (5.0, 2.0), (6.0, 6.0), (8.0, 3.0)];
    let n = cities.len();
    let mut distances = vec![vec![0.0; n]; n];
    for i in 0..n {
        for j in 0..n {
            let dx = cities[i].0 - cities[j].0;
            let dy = cities[i].1 - cities[j].1;
            distances[i][j] = (dx * dx + dy * dy).sqrt();
        }
    }
    let problem = Tsp { distances: distances.clone() };

    let mut opt = AntColonyTsp::new(AntColonyTspConfig {
        ants: 20,
        iterations: 200,
        alpha: 1.0,
        beta: 5.0,
        evaporation: 0.5,
        deposit: 1.0,
        distances,
        seed: 42,
    });

    let r = opt.run(&problem);
    let best = r.best.unwrap();
    println!("best tour length: {:.3}", best.evaluation.objectives[0]);
}

alpha weights pheromone influence and beta weights the heuristic (1 / distance). evaporation is the per-iteration pheromone decay. The classic Dorigo paper uses alpha = 1, beta = 2..5, evaporation = 0.1..0.5.

Tiny baseline: SA + SwapMutation

The smallest possible permutation optimizer — one starting decision, no population, one mutation operator. Good as a sanity-check baseline.

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

struct JobShop {
    process_times: Vec<f64>,
}
impl Problem for JobShop {
    type Decision = Vec<usize>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("weighted_completion")])
    }
    fn evaluate(&self, schedule: &Vec<usize>) -> Evaluation {
        let cost: f64 = schedule.iter().enumerate()
            .map(|(i, &job)| (i as f64 + 1.0) * self.process_times[job])
            .sum();
        Evaluation::new(vec![cost])
    }
}

let times = vec![3.0, 1.5, 4.2, 2.7, 5.1];
let n = times.len();
let problem = JobShop { process_times: times };

// SimulatedAnnealing expects exactly one initial decision.
struct OneShuffle { n: usize }
impl Initializer<Vec<usize>> for OneShuffle {
    fn initialize(&mut self, _size: usize, rng: &mut Rng) -> Vec<Vec<usize>> {
        use rand::seq::SliceRandom;
        let mut p: Vec<usize> = (0..self.n).collect();
        p.shuffle(rng);
        vec![p]
    }
}

let mut opt = SimulatedAnnealing::new(
    SimulatedAnnealingConfig {
        iterations: 2000,
        initial_temperature: 5.0,
        final_temperature: 1e-3,
        seed: 7,
    },
    OneShuffle { n },
    SwapMutation,
);
let r = opt.run(&problem);
let best = r.best.unwrap();
println!("best cost: {:.3}", best.evaluation.objectives[0]);
}

Custom neighborhoods: Tabu Search

When you want full control of the move set (e.g., systematic 2-opt for TSP, or insert-and-shift for scheduling), Tabu Search takes your own neighbor function.

use heuropt::prelude::*;
let neighbors = |x: &Vec<usize>, _rng: &mut Rng| -> Vec<Vec<usize>> {
    // All 2-opt neighbors of x.
    let mut out = Vec::new();
    for i in 0..x.len() {
        for j in (i + 2)..x.len() {
            let mut child = x.clone();
            child[i + 1..=j].reverse();
            out.push(child);
        }
    }
    out
};
// Pass `neighbors` to TabuSearch::new(...).

When to use which approach

Multi-objective combinatorial problems

Real combinatorial problems usually have more than one cost. A TSP where every edge has both distance and time; a job-shop where you care about makespan, flow time, and tardiness; a knapsack with two profit metrics and a single weight budget. The decision type is still combinatorial — a permutation, a bitstring — but the objective is a vector, and the answer is a Pareto front rather than a single best.

heuropt's NSGA-II and NSGA-III are fully generic over the decision type. You don't need a separate "combinatorial NSGA" — just plug in the right initializer and variation operators for your encoding.

This recipe walks through three patterns:

• Bi-objective TSP with NSGA-II (Pareto front of two distance matrices over the same cities)
• Bi-objective 0/1 knapsack with NSGA-II (binary encoding)
• 3-objective JSS with NSGA-III (the many-objective successor)

For the single-objective permutation toolkit it builds on, see Optimize a permutation.

Bi-objective TSP

This is the canonical multi-objective combinatorial benchmark (Lust–Teghem 2010). Two TSP instances on the same city set define two distance matrices A and B; the search trades off length under A versus length under B.

use heuropt::prelude::*;
use heuropt::metrics::hypervolume_2d;

struct BiObjectiveTsp {
    dist_a: Vec<Vec<f64>>,
    dist_b: Vec<Vec<f64>>,
}

impl BiObjectiveTsp {
    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 BiObjectiveTsp {
    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),
        ])
    }
}

fn main() {
    let n: usize = 25;
    let dist_a = vec![vec![0.0_f64; n]; n]; // your matrix A
    let dist_b = vec![vec![0.0_f64; n]; n]; // your matrix B
    let problem = BiObjectiveTsp { dist_a, dist_b };

    let mut optimizer = Nsga2::new(
        Nsga2Config {
            population_size: 200,
            generations: 600,
            seed: 11,
        },
        ShuffledPermutation { n },
        CompositeVariation {
            crossover: EdgeRecombinationCrossover,
            mutation: InversionMutation,
        },
    );
    let result = optimizer.run(&problem);

    println!("Pareto-front size: {}", result.pareto_front.len());

    // Hypervolume against a generous reference point (larger than any
    // length you'd reasonably see). Use this as the single-number
    // quality metric for the run.
    let ref_point = [40_000.0, 40_000.0];
    let hv = hypervolume_2d(&result.pareto_front, &problem.objectives(), ref_point);
    println!("Hypervolume vs. {:?}: {:.0}", ref_point, hv);
}

EdgeRecombinationCrossover (ERX) is the standout crossover for TSP. On a 25-city bi-objective instance it produces about twice the front diversity of OX, PMX, or CX — see examples/tsp_operators_compare.rs for a head-to-head benchmark.

Bi-objective 0/1 knapsack — Vec<bool> decisions

NSGA-II works over Vec<bool> the same way. The Zitzler–Thiele bi-objective knapsack is the textbook benchmark: each item has two profit values and a single weight; you maximize both profits under one capacity constraint.

use heuropt::prelude::*;
use rand::Rng as _;

const N_ITEMS: usize = 30;

struct BiKnapsack {
    profits_a: Vec<f64>,
    profits_b: Vec<f64>,
    weights: Vec<f64>,
    capacity: f64,
}

impl Problem for BiKnapsack {
    type Decision = Vec<bool>;

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![
            Objective::maximize("profit_A"),
            Objective::maximize("profit_B"),
        ])
    }

    fn evaluate(&self, take: &Vec<bool>) -> Evaluation {
        let (pa, pb, w) = take.iter().enumerate().fold(
            (0.0_f64, 0.0_f64, 0.0_f64),
            |(pa, pb, w), (i, &t)| {
                if t {
                    (pa + self.profits_a[i], pb + self.profits_b[i], w + self.weights[i])
                } else {
                    (pa, pb, w)
                }
            },
        );
        // Standard heuristic-MO constraint handling: penalize weight
        // overruns heavily so the recovered front is feasible.
        let penalty = 1000.0 * (w - self.capacity).max(0.0);
        Evaluation::new(vec![pa - penalty, pb - penalty])
    }
}

/// Each bit 50/50 independently.
#[derive(Clone, Copy)]
struct RandomBinary { n: usize }
impl Initializer<Vec<bool>> for RandomBinary {
    fn initialize(&mut self, size: usize, rng: &mut Rng) -> Vec<Vec<bool>> {
        (0..size).map(|_| (0..self.n).map(|_| rng.random_bool(0.5)).collect()).collect()
    }
}

/// One-point crossover for binary chromosomes.
#[derive(Default)]
struct OnePointCrossoverBool;
impl Variation<Vec<bool>> for OnePointCrossoverBool {
    fn vary(&mut self, parents: &[Vec<bool>], rng: &mut Rng) -> Vec<Vec<bool>> {
        let (p1, p2) = (&parents[0], &parents[1]);
        let n = p1.len();
        let cut = rng.random_range(1..n);
        let mut c1 = Vec::with_capacity(n);
        let mut c2 = Vec::with_capacity(n);
        c1.extend_from_slice(&p1[..cut]); c1.extend_from_slice(&p2[cut..]);
        c2.extend_from_slice(&p2[..cut]); c2.extend_from_slice(&p1[cut..]);
        vec![c1, c2]
    }
}

fn main() {
  let profits_a = vec![0.0; N_ITEMS];
  let profits_b = vec![0.0; N_ITEMS];
  let weights   = vec![0.0; N_ITEMS];
    let problem = BiKnapsack {
        profits_a, profits_b, weights,
        capacity: 750.0, // ~half the total weight
    };

    let mut optimizer = Nsga2::new(
        Nsga2Config {
            population_size: 120,
            generations: 400,
            seed: 19,
        },
        RandomBinary { n: N_ITEMS },
        CompositeVariation {
            crossover: OnePointCrossoverBool,
            mutation: BitFlipMutation { probability: 1.0 / N_ITEMS as f64 },
        },
    );
    let result = optimizer.run(&problem);
    println!("Pareto-front size: {}", result.pareto_front.len());
}

Two things worth noting:

• OnePointCrossoverBool and RandomBinary are defined locally. They're tiny and common — a future PR could lift them into the library, but for now you write them inline.
• Constraint handling is a penalty. The factor 1000.0 is chosen so that even a 1-unit overrun beats any feasible solution by more than the entire profit range; the recovered front is entirely feasible. This is the standard heuristic-MO pattern (Deb 2001) and cheaper than a hard repair operator.

Three-objective JSS with NSGA-III

NSGA-III is designed for ≥ 3 objectives. NSGA-II's crowding distance degrades when most of the population is mutually non-dominated, which is the rule rather than the exception in higher dimensions; NSGA-III uses reference-point niching instead.

The example below adds tardiness to the standard (makespan, flow time) JSS pair. Tardiness needs due dates; the common heuristic is dⱼ = 1.3 × sum_of_processing_times(j).

use heuropt::prelude::*;
use rand::Rng as _;

const N_JOBS: usize = 10;
const N_MACHINES: usize = 5;

struct La01ThreeObjective {
    routing: [[usize; N_MACHINES]; N_JOBS],
    times:   [[f64;   N_MACHINES]; N_JOBS],
    due:     [f64; N_JOBS],
}

impl Problem for La01ThreeObjective {
    type Decision = Vec<usize>;

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![
            Objective::minimize("makespan"),
            Objective::minimize("total_flow_time"),
            Objective::minimize("total_tardiness"),
        ])
    }

    fn evaluate(&self, schedule: &Vec<usize>) -> Evaluation {
        let mut job_next = [0_usize; N_JOBS];
        let mut job_clock = [0.0_f64; N_JOBS];
        let mut machine_clock = [0.0_f64; N_MACHINES];
        for &job in schedule {
            let k = job_next[job];
            let m = self.routing[job][k];
            let t = self.times[job][k];
            let start = job_clock[job].max(machine_clock[m]);
            let end = start + t;
            job_clock[job] = end;
            machine_clock[m] = end;
            job_next[job] = k + 1;
        }
        let makespan = machine_clock.iter().cloned().fold(0.0_f64, f64::max);
        let flow_time: f64 = job_clock.iter().sum();
        let tardiness: f64 = job_clock.iter().zip(self.due.iter())
            .map(|(&c, &d)| (c - d).max(0.0))
            .sum();
        Evaluation::new(vec![makespan, flow_time, tardiness])
    }
}

/// Mix Insertion and Scramble per call — both preserve the multiset,
/// giving the search access to two complementary neighborhood moves.
#[derive(Default)]
struct InsertionOrScramble;
impl Variation<Vec<usize>> for InsertionOrScramble {
    fn vary(&mut self, parents: &[Vec<usize>], rng: &mut Rng) -> Vec<Vec<usize>> {
        if rng.random_bool(0.5) {
            InsertionMutation.vary(parents, rng)
        } else {
            ScrambleMutation.vary(parents, rng)
        }
    }
}

fn main() {
  let routing = [[0; N_MACHINES]; N_JOBS];
  let times   = [[0.0; N_MACHINES]; N_JOBS];
    let due = std::array::from_fn::<f64, N_JOBS, _>(
        |j| 1.3 * times[j].iter().sum::<f64>(),
    );
    let problem = La01ThreeObjective { routing, times, due };

    let mut optimizer = Nsga3::new(
        Nsga3Config {
            population_size: 120,
            generations: 600,
            reference_divisions: 12, // 91 Das-Dennis points in 3-D
            seed: 9,
        },
        ShuffledMultisetPermutation::new(vec![N_MACHINES; N_JOBS]),
        // Drop in a local PrecedenceOrderCrossover (POX) here for the
        // crossover slot if you want stronger mixing; see the
        // permutation recipe for the implementation.
        InsertionOrScramble,
    );
    let result = optimizer.run(&problem);

    println!("Pareto-front size: {}", result.pareto_front.len());
}

A few NSGA-III tips:

• reference_divisions controls how many reference points the algorithm spreads across the front. For M objectives, the Das–Dennis construction produces C(divisions + M - 1, M - 1) reference points. For M = 3 and divisions = 12 that's 91 points; pick a population size ≥ that.
• PrecedenceOrderCrossover (POX) belongs in the crossover slot for JSS. The strict-permutation crossovers (OX, PMX, CX, ERX) break the operation-string multiset. See Optimize a permutation for the local POX definition.

Comparing operators by hypervolume

For Pareto-front problems, single-objective fitness is the wrong comparison metric. Use hypervolume instead — the dominated area under the front, against a fixed reference point.

use heuropt::metrics::hypervolume_2d;

let ref_point = [40_000.0, 40_000.0]; // worse than anything you expect

for (name, crossover) in &[
    ("OX",  Box::new(OrderCrossover)             as Box<dyn Variation<Vec<usize>>>),
    ("PMX", Box::new(PartiallyMappedCrossover)   as _),
    ("CX",  Box::new(CycleCrossover)             as _),
    ("ERX", Box::new(EdgeRecombinationCrossover) as _),
] {
    let result = run_nsga2_with_crossover(crossover);
    let hv = hypervolume_2d(&result.pareto_front, &problem.objectives(), ref_point);
    println!("{name:>3}: hv = {hv:.0}");
}

This is the pattern in examples/tsp_operators_compare.rs. On the KroAB-25 instance it ranks ERX > OX > PMX > CX by hypervolume.

For ≥ 3 objectives, hypervolume in N dimensions is exponentially expensive; use hypervolume_2d when you can collapse to two objectives for the metric, or sample-based hypervolume from [HypE] otherwise.

Pareto-front tips

Constrain your search with Repair

heuropt models constraints with a single constraint_violation scalar on each Evaluation. That works for soft penalties. When constraints are hard and the search keeps generating infeasible decisions, the better pattern is repair: project each candidate back into the feasible region every time it leaves.

The Repair<D> trait is the abstraction. Two impls ship in the box; you can write your own for arbitrary geometry.

Built-in: ClampToBounds

For per-axis box constraints (lo ≤ xᵢ ≤ hi), pair ClampToBounds with any Variation to get a bounds-aware variant for free.

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

let bounds = vec![(-5.0, 5.0); 3];

// Without bounds, GaussianMutation can step outside the search box.
// ClampToBounds projects each variable back in.
let mut sigma = GaussianMutation { sigma: 0.5 };
let mut clamp = ClampToBounds::new(bounds.clone());

let mut rng = rng_from_seed(42);
let parent = vec![4.9, -4.9, 0.0];
let mut child = sigma.vary(std::slice::from_ref(&parent), &mut rng).pop().unwrap();
clamp.repair(&mut child);
// every entry of `child` is now within [-5, 5].
}

ClampToBounds is idempotent: applying it twice is the same as applying it once.

For most real problems you'd just use BoundedGaussianMutation which combines both in one operator.

Built-in: ProjectToSimplex

For budget constraints — "the components must sum to a fixed total and be non-negative" — ProjectToSimplex projects onto the probability simplex (or any scaled simplex).

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

let mut proj = ProjectToSimplex::new(1.0); // probability simplex
let mut x = vec![0.6, 0.5, -0.1, 0.3];     // sum 1.3, one negative
proj.repair(&mut x);
// x now sums to 1.0 and every entry is ≥ 0.
let s: f64 = x.iter().sum();
debug_assert!((s - 1.0).abs() < 1e-12);
debug_assert!(x.iter().all(|&v| v >= 0.0));
}

Use this for portfolio / resource-allocation problems where the decision is a vector of weights that must sum to a budget.

Custom repair

Anything that takes a &mut Vec<f64> (or any &mut D for your custom decision type) and returns a feasible version is a valid Repair. Implement the trait directly:

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

/// Force the largest variable to be at least `min_largest`.
struct AtLeastOneActive { min_largest: f64 }

impl Repair<Vec<f64>> for AtLeastOneActive {
    fn repair(&mut self, x: &mut Vec<f64>) {
        let max_idx = x.iter()
            .enumerate()
            .fold(0, |best, (i, &v)| {
                if v > x[best] { i } else { best }
            });
        if x[max_idx] < self.min_largest {
            x[max_idx] = self.min_largest;
        }
    }
}
}

Stochastic-ranking selection

When the feasible region is narrow — most of the search space is infeasible — the strict "feasibles always beat infeasibles" rule traps the search outside it. Runarsson & Yao's stochastic ranking breaks the trap by, on each pairwise comparison, using a probabilistic "compare by objective" instead of "compare by feasibility" with a small probability pf:

use heuropt::selection::tournament::stochastic_ranking_select;

let picks = stochastic_ranking_select(
    &population,
    &objectives,
    0.45,             // pf — Runarsson & Yao's canonical value
    count,
    &mut rng,
);

This is a drop-in replacement for tournament_select_single_objective in your custom optimizer or in a forked algorithm.

When to use which

Pick one answer off a Pareto front

A multi-objective optimizer hands you a front — a Pareto-optimal trade-off curve — not a single answer. Eventually you have to pick one point off it. There are several principled ways to do that; this recipe covers the most common: the a-posteriori weighted decision rule.

The pattern: optimize without baking your preferences into the search, then apply your preferences as a scoring function over the front.

This is exactly the pattern from examples/jiggly_tuning.rs (the USB-jiggler firmware tuning example).

The shape

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

struct Cost;
impl Problem for Cost {
    type Decision = Vec<f64>;
    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![Objective::minimize("a"), Objective::minimize("b"), Objective::minimize("c")])
    }
    fn evaluate(&self, _x: &Vec<f64>) -> Evaluation { Evaluation::new(vec![0.0,0.0,0.0]) }
}

let problem = Cost;
let mut opt = Nsga2::new(
    Nsga2Config { population_size: 100, generations: 200, seed: 42 },
    RealBounds::new(vec![(-1.0, 1.0); 4]),
    CompositeVariation {
        crossover: SimulatedBinaryCrossover::new(vec![(-1.0, 1.0); 4], 15.0, 0.5),
        mutation:  PolynomialMutation::new(vec![(-1.0, 1.0); 4], 20.0, 1.0),
    },
);
let result = opt.run(&problem);

// 1. Get the Pareto front.
let front = &result.pareto_front;

// 2. Define your preferences as a scoring function over (oriented)
//    objective values. Lower score = preferred.
let space = problem.objectives();
let weights = [1.0, 2.0, 0.5];

let scored: Vec<(f64, &Candidate<Vec<f64>>)> = front.iter()
    .map(|c| {
        let oriented = space.as_minimization(&c.evaluation.objectives);
        let score: f64 = oriented.iter().zip(&weights)
            .map(|(v, w)| v * w)
            .sum();
        (score, c)
    })
    .collect();

// 3. Pick the lowest-scoring point.
let best = scored.iter()
    .min_by(|a, b| a.0.partial_cmp(&b.0).unwrap())
    .unwrap();

println!("picked: {:?} with weighted score {:.3}",
    best.1.evaluation.objectives, best.0);
}

as_minimization returns the objective vector with maximized axes flipped to negative — so a single set of positive weights does the right thing whether each axis is min or max.

Why a-posteriori vs a-priori weighting

If you know your weights up front, you could just optimize the weighted sum directly with a single-objective algorithm. Why bother with the multi-objective dance?

Two reasons:

1. Weighted sum can't reach concave parts of the Pareto front. Any single-objective optimization with a linear scalarization converges to a point at the boundary of the convex hull. Concave front segments are unreachable. The multi-objective optimizer finds them.
2. Weights are usually wrong on the first try. Optimizing the front first lets you see what's actually possible before deciding how much each axis is worth. Run once, look at the trade-offs, adjust weights.

Penalty terms beyond linear weights

The jiggly example also adds a hinge penalty — a term that's zero inside an acceptable region and grows quadratically once you exceed some hard cap. Useful when one axis is "soft up to X, hard cap at Y":

#![allow(unused)]
fn main() {
fn hinge(x: f64, soft_cap: f64, hard_cap: f64) -> f64 {
    if x <= soft_cap { 0.0 }
    else if x >= hard_cap { f64::INFINITY }
    else {
        let t = (x - soft_cap) / (hard_cap - soft_cap);
        100.0 * t * t
    }
}
}

Compose linear weights + hinge penalties and you have a flexible scoring function over the front without re-running the optimizer.

Other strategies

• Knee point. Pick the point where small gains in one axis cost large losses in another — the "elbow" of the trade-off curve. Knea explicitly biases the search toward knees during the run.
• Reference-direction. Pick the point closest to a desired trade-off direction (a unit vector in objective space). Moead / Nsga3 use this internally during search; you can apply it post-hoc the same way.
• Random / interactive selection. Show the front to a user (perhaps via a plotting library), let them pick.

The right pick depends on the problem; the front itself doesn't prescribe one.

Explore your results in a webapp

Real Pareto fronts have 50–200+ candidates spanning 2–7+ objectives. Reading them as a wall of numbers in a terminal scales badly. Drop the result into heuropt-explorer to filter, brush, pin, and rank candidates interactively in the browser — parallel coordinates, scatter plots, sortable table, range filters, weighted ranking, knee-point detection.

This recipe shows the export side. The webapp is a static page; no install needed beyond a browser.

Enable the serde feature

[dependencies]
heuropt = { version = "0.10", features = ["serde"] }

The export uses serde_json under the hood, so the explorer module is gated on the existing serde feature.

Enrich your Problem (optional but worth it)

Two places to add display metadata that flows through to the explorer's axis labels and tooltips:

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

struct PickACar;

impl Problem for PickACar {
    type Decision = Vec<f64>;

    fn objectives(&self) -> ObjectiveSpace {
        ObjectiveSpace::new(vec![
            // `name` is the canonical short ID; `label` and `unit`
            // are display-only. The explorer renders axes as
            // `Price ($k)` instead of just `price`.
            Objective::minimize("price").with_label("Price").with_unit("$k"),
            Objective::minimize("zero_to_sixty").with_label("0-60 mph").with_unit("s"),
            Objective::minimize("fuel").with_label("Fuel").with_unit("gal/100mi"),
            Objective::minimize("noise").with_label("Idle noise").with_unit("dB"),
        ])
    }

    fn decision_schema(&self) -> Vec<DecisionVariable> {
        // Optional: provide name/label/unit/bounds per decision-variable
        // slot. If you skip this, the exporter falls back to `x[0]`,
        // `x[1]`, … with no units or bounds.
        vec![
            DecisionVariable::new("displacement")
                .with_label("Engine size").with_unit("L").with_bounds(1.0, 6.0),
            DecisionVariable::new("weight")
                .with_label("Curb weight").with_unit("kg").with_bounds(1100.0, 2200.0),
            DecisionVariable::new("drag")
                .with_label("Drag coefficient").with_unit("Cd").with_bounds(0.20, 0.40),
        ]
    }

    fn evaluate(&self, x: &Vec<f64>) -> Evaluation {
        // ... compute objectives ...
      Evaluation::new(vec![0.0, 0.0, 0.0, 0.0])
    }
}
}

Both Objective::with_label / with_unit and Problem::decision_schema are entirely optional — the rest of heuropt doesn't read them. They exist so the exported JSON describes itself well enough for a display tool to render readable axes.

Run the optimizer and write the JSON

The simplest call (no algorithm metadata in the export):

use heuropt::prelude::*;

let result = optimizer.run(&problem);
heuropt::explorer::ExplorerExport::from_result(&problem, &result)
    .to_file("results.json")
    .unwrap();

The richer call — pulls algorithm name + seed automatically from the AlgorithmInfo trait that every built-in algorithm implements:

use heuropt::prelude::*;

let started = std::time::Instant::now();
let result = optimizer.run(&problem);

let export = heuropt::explorer::ExplorerExport::from_result(&problem, &result)
    .with_algorithm_info(&optimizer)
    .with_problem_name("Pick a car")
    .with_wall_clock(started.elapsed().as_secs_f64());
export.to_file("results.json").unwrap();

There's also a one-liner if you don't need to set extra metadata:

heuropt::explorer::to_file("results.json", &problem, &optimizer, &result).unwrap();

Open it in the explorer

Visit https://swaits.github.io/heuropt-explorer/ and drag the JSON file onto the page. The explorer reads the units and labels you attached and renders parallel-coordinates / scatter / table views that respect them. Brushing on any axis filters the others; pinned candidates stay highlighted; the weight sliders let you rank the front by your priorities.

What's in the file

The full schema is documented in heuropt::explorer::ExplorerExport. The shape:

{
  "schema_version": 1,
  "run": {
    "problem_name": "Pick a car",
    "algorithm": "Nsga3",
    "seed": 42,
    "wall_clock_seconds": 0.097,
    "evaluations": 20100,
    "generations": 200
  },
  "objectives": [
    { "name": "price", "direction": "Minimize", "label": "Price", "unit": "$k" },
    ...
  ],
  "decision_variables": [
    { "name": "displacement", "label": "Engine size", "unit": "L", "min": 1.0, "max": 6.0 },
    ...
  ],
  "candidates": [
    {
      "decision": [1.0, 1505.0, 0.35],
      "objectives": [13.0, 7.0, 3.17, 63.0],
      "constraint_violation": 0.0,
      "feasible": true,
      "front_rank": 0,
      "in_pareto_front": true
    },
    ...
  ]
}

front_rank is computed by non_dominated_sort once at export time — 0 means on the Pareto front, higher numbers indicate deeper layers.

Custom decision types

Out of the box, Vec<f64>, Vec<bool>, Vec<usize>, and Vec<i64> work as decisions. For a custom decision type, implement heuropt::explorer::ToDecisionValues:

struct MyDecision { color: String, count: u32 }

impl heuropt::explorer::ToDecisionValues for MyDecision {
    fn to_decision_values(&self) -> Vec<serde_json::Value> {
        vec![
            serde_json::Value::String(self.color.clone()),
            serde_json::Value::Number(self.count.into()),
        ]
    }
}

The explorer renders strings as categorical axes and numbers as continuous.

Worked example

examples/pick_a_car.rs ships with the crate. It implements the problem above, runs NSGA-III for 200 generations, and writes pick_a_car.json ready to load:

cargo run --release --example pick_a_car --features serde

Write your own algorithm

Implement Optimizer<P> and you're done. There are no other traits to think about, no internal hooks to register. The example walks through a tiny hill-climber that reads almost identically to the canonical pseudocode.

The trait

pub trait Optimizer<P>
where
    P: Problem,
{
    fn run(&mut self, problem: &P) -> OptimizationResult<P::Decision>;
}

That's it. You own your config, your RNG, your main loop, and your OptimizationResult construction.

A minimal hill-climber

#![allow(unused)]
fn main() {
use heuropt::prelude::*;

pub struct MyHillClimber<I, V> {
    pub iterations: usize,
    pub seed: u64,
    pub initializer: I,
    pub variation: V,
}

impl<P, I, V> Optimizer<P> for MyHillClimber<I, V>
where
    P: Problem,
    P::Decision: Clone,
    I: Initializer<P::Decision>,
    V: Variation<P::Decision>,
{
    fn run(&mut self, problem: &P) -> OptimizationResult<P::Decision> {
        let mut rng = rng_from_seed(self.seed);
        let objectives = problem.objectives();
        assert!(objectives.is_single_objective(), "MyHillClimber is single-objective only");

        // Start with one initial decision.
        let init_decisions = self.initializer.initialize(1, &mut rng);
        let init = init_decisions.into_iter().next().unwrap();
        let mut current = Candidate::new(init.clone(), problem.evaluate(&init));
        let mut evaluations: usize = 1;

        for _ in 0..self.iterations {
            let children = self.variation.vary(std::slice::from_ref(&current.decision), &mut rng);
            for child_decision in children {
                let child_eval = problem.evaluate(&child_decision);
                evaluations += 1;
                let child = Candidate::new(child_decision, child_eval);
                if better(&child.evaluation, &current.evaluation, &objectives) {
                    current = child;
                }
            }
        }

        let pareto_front = vec![current.clone()];
        let best = Some(current.clone());
        OptimizationResult::new(
            Population::new(vec![current]),
            pareto_front,
            best,
            evaluations,
            self.iterations,
        )
    }
}

fn better(a: &Evaluation, b: &Evaluation, objectives: &ObjectiveSpace) -> bool {
    let am = objectives.as_minimization(&a.objectives);
    let bm = objectives.as_minimization(&b.objectives);
    am[0] < bm[0]
}
}

Things to notice

• Rng is one concrete type. No generics — call rng_from_seed and pass &mut rng everywhere it's needed.
• Initializer<D> sources the starting point(s).
• Variation<D> generates children from parents. For the hill-climber it's called with one parent.
• OptimizationResult carries the final population, the Pareto front (just the best for single-objective), the best candidate, the total evaluations, and the iteration count.
• as_minimization flips maximize-axis values so your comparison logic only ever needs to deal with "lower is better."

Adding parallel evaluation

If your algorithm batch-evaluates candidates per generation, use the crate's internal helper. From inside heuropt source you can call evaluate_batch(problem, decisions); from outside you'd use rayon directly behind a feature flag, the same way the built-in algorithms do.

#[cfg(feature = "parallel")]
fn batch_eval<P>(problem: &P, decisions: Vec<P::Decision>) -> Vec<Candidate<P::Decision>>
where P: Problem + Sync, P::Decision: Send,
{
    use rayon::prelude::*;
    decisions.into_par_iter()
        .map(|d| Candidate::new(d.clone(), problem.evaluate(&d)))
        .collect()
}

#[cfg(not(feature = "parallel"))]
fn batch_eval<P>(problem: &P, decisions: Vec<P::Decision>) -> Vec<Candidate<P::Decision>>
where P: Problem,
{
    decisions.into_iter()
        .map(|d| Candidate::new(d.clone(), problem.evaluate(&d)))
        .collect()
}

To stay bit-identical between serial and parallel modes, keep the RNG and selection on the main thread; only the evaluations run in parallel.

What's not in the trait

• No iteration / step API. The optimizer owns its loop.
• No callbacks. A future minor release may add an observer hook; for now you'd run the algorithm to completion and process the result.
• No error type. Invalid configuration panics with a clear message; this matches the style of the built-in algorithms.
• No async on the trait. Optimizer<P> is synchronous. For async evaluation, implement AsyncProblem on your problem and use the run_async(&problem, concurrency) method that comes with the async feature. See the Async evaluation cookbook recipe.

The smallness is the point: you should be able to read a built-in algorithm and write your own in an afternoon. See examples/custom_optimizer.rs for a slightly more polished version of the hill-climber above.

Comparison with other libraries

heuropt is one of many heuristic-optimization libraries. This chapter is an honest, opinionated comparison to help you choose.

The columns:

• Lang — primary implementation language.
• Algorithms — rough catalog count.
• Multi-obj — built-in support for Pareto-based multi-objective optimization.
• Surrogates — built-in Bayesian / TPE / multi-fidelity.
• Determinism — seeded reproducibility as a first-class property.
• Async / async-eval — first-class async runtime support.

When to pick heuropt

• You're working in Rust and want a single, dependency-light crate for evolutionary / metaheuristic optimization.
• You need multi-objective or many-objective algorithms (12+ Pareto-aware methods in the catalog) AND you don't want to glue Python into your Rust pipeline.
• You want bit-identical determinism: same seed produces same output, on every machine, across releases unless explicitly noted otherwise.
• You want a small, readable codebase — every algorithm is written for clarity, no trait-object plumbing, no GATs in user- facing APIs. Reading Random Search should be enough to write a new optimizer.
• You have IO-bound evaluations — calling an HTTP service, an RPC, or a subprocess — and want first-class async fn evaluate support. heuropt is the only mainstream optimization library that ships this (see Async evaluation).

When not to pick heuropt

• You need gradient-based optimization. Use argmin (Rust) or scipy.optimize (Python) — heuropt is gradient-free by design.
• You need GPU-accelerated evaluations. heuropt's evaluate function runs on CPU; use Python (jax/torch) or roll your own GPU pipeline.
• You need distributed multi-machine evaluation. heuropt parallelizes within one process via rayon. Distribution is up to you (split the seeds across machines, aggregate).
• You're comfortable in Python and pymoo / optuna already cover your problem. heuropt's value-add over pymoo is mostly that it's Rust — if that doesn't matter to you, the Python ecosystem has more battle-tested integrations.

Algorithm coverage at a glance

heuropt covers the same major Pareto MOEAs as pymoo and MOEA Framework: NSGA-II/III, SPEA2, IBEA, MOEA/D, MOPSO, SMS-EMOA, HypE, AGE-MOEA, GrEA, KnEA, RVEA, PESA-II, ε-MOEA, PAES.

The expensive-evaluation regime: Bayesian Optimization + TPE + Hyperband. This is comparable to optuna's coverage but in pure Rust.

The single-objective continuous catalog (CMA-ES, IPOP-CMA-ES, sNES, DE, PSO, GA, TLBO, (1+1)-ES, Nelder-Mead, Random Search, Hill Climber, Simulated Annealing) covers the canonical baselines and several modern variants.

What heuropt does not ship that some libraries do:

• Re-themed metaphor metaheuristics (Whale Optimization, Grey Wolf, Bat, Firefly, Harris Hawks, etc.). These are cut from the catalog deliberately — they are mostly DE/PSO with new names. If you specifically need one, please open an issue with citations.
• Non-evolutionary global optimizers like dual annealing or basin-hopping (use scipy.optimize for those).
• A web UI / dashboard like optuna's. heuropt is library-only.

Speed

heuropt's hot paths (Pareto utilities, hypervolume, key inner loops) are heavily optimized — see the perf entry in the v0.4.0 CHANGELOG. On the comparison harness in examples/compare.rs (10-seed mean, 30 000 evaluations on DTLZ2), the total wall-clock time across 12 algorithms is ~5 seconds. Per-algorithm timings are in examples/compare-results.md.

For comparison-shopping speed against Python libraries, the gap is typically 10×–100× in heuropt's favor for compute-bound evaluate functions, because Rust skips the Python-loop overhead. If your evaluate calls into NumPy/PyTorch and those are the bottleneck, the gap shrinks substantially.

Honest weakness: ecosystem

The biggest thing pymoo / optuna / DEAP have that heuropt doesn't: community + plug-ins + tutorials. They've been around longer and have rich third-party integrations (visualization, MLflow, Hyperband+BO hybrids, distributed runners). heuropt is younger; the core is solid but the ecosystem is small.

If you adopt heuropt and miss a thing, the project is small enough that contributions land fast. See CONTRIBUTING.md.

Stability and SemVer

heuropt is pre-1.0. The public API may change between minor versions. This page sets explicit expectations.

What "public API" means in heuropt

The crate's public surface is everything re-exported from heuropt::prelude plus the items reachable from heuropt::core, heuropt::traits, heuropt::operators, heuropt::algorithms, heuropt::pareto, heuropt::metrics, and heuropt::selection.

Items in heuropt::internal (e.g. the Cholesky / eigendecomposition helpers) are not public API. They may change between any two versions — use them at your own risk.

SemVer in heuropt 0.x

While we are pre-1.0:

• Minor bumps (0.10 → 0.11) may break the public API. The CHANGELOG calls out everything that changed, and a migration guide in this book documents the move.
• Patch bumps (0.10.0 → 0.10.1) only contain bug fixes, performance improvements, and additive non-breaking features. No deprecations, no removals.

What's actually likely to change before 1.0

In rough order of likelihood:

1. Algorithm config structs may gain fields. All current configs are public-field structs; adding a non-Default field is a breaking change. We may switch to builder patterns to avoid this class of break, or we may add #[non_exhaustive].
2. Some operators may move between operators and pareto as the boundary between "things that produce candidates" and "Pareto utilities" gets clearer.

What is not likely to change:

• The Problem trait shape.
• The AsyncProblem / AsyncPartialProblem trait shapes.
• The Variation / Initializer / Repair traits.
• The Optimizer<P> trait — single run method, no callbacks.
• The Evaluation / Candidate / Population / OptimizationResult data types.
• The seeded determinism property.

What "bit-identical" means for stability

heuropt promises that a given algorithm + seed + config produces the same numeric output on the same minor version of heuropt.

Across minor versions, output may change if an algorithm's implementation changes (e.g. a perf rewrite that reorders floating-point operations, or a new feature that changes the RNG-consumption pattern). The CHANGELOG calls this out explicitly when it happens. As of v0.8, the entire history of perf optimizations has been bit-identical against the v0.3.0 reference.

MSRV (minimum supported Rust version)

heuropt's MSRV is 1.85 as of v0.10. This is tested in CI against every PR.

MSRV bumps are treated as patch-bump-eligible (they don't break the public API). When the MSRV is bumped, the CHANGELOG entry for that release will note the new MSRV.

Feature-flag stability

The current optional features:

• serde — adds Serialize / Deserialize derives on the core data types.
• parallel — rayon-backed parallel population evaluation.
• async — AsyncProblem + AsyncPartialProblem traits, plus a run_async method on every algorithm in the catalog, for IO-bound evaluations.

Features added in 0.x can be renamed or removed in any minor bump that documents the change. Removing a feature is treated like a breaking API change.

How to track changes

• CHANGELOG.md — the canonical record of changes per release.
• Migration guides — per-release, in this book at migration.
• GitHub releases — each tag has release notes.
• Watch the repo — https://github.com/swaits/heuropt — to be notified of new releases.

Migration guides

Per-release notes for upgrading between heuropt versions. Skip the sections that don't apply to your starting version.

To 0.10

From 0.9.x

Almost additive. Bumping heuropt = "0.10" recompiles without touching most code. The one breaking change is the value returned by AlgorithmInfo::name():

If you pattern-matched on those strings (e.g. for branching display logic), update to the new canonical strings. They now match the literature and will be stable going forward.

What's new and additive:

• AlgorithmInfo::full_name(&self) -> &'static str — academic long form ("Non-dominated Sorting Genetic Algorithm II"). Defaults to name() for algorithms whose long and short forms coincide.
• ExplorerExport's RunMeta gained algorithm_full_name: Option<String>. Schema version stays at 1 (the new field is #[serde(default)]); display tools can use the long form as a hover tooltip on the short name.

To 0.9

From 0.8.x

Additive only. Bumping heuropt = "0.9" works for all 0.8.x code untouched. The new surfaces ship behind the existing serde feature.

What's new:

• heuropt::explorer module (gated on serde) — turns an OptimizationResult into a self-describing JSON file that the heuropt-explorer webapp can load. See the Explore your results recipe.
• Objective gained optional label and unit fields with fluent builders .with_label("…") / .with_unit("…"). Existing Objective::minimize("…") / Objective::maximize("…") are unchanged. The serde representation is forward- and backward- compatible (new fields are #[serde(default)]).
• Problem trait gained a default-empty fn decision_schema(&self) -> Vec<DecisionVariable> method. Existing impls compile untouched; override it to provide pretty names / labels / units / bounds for the explorer.
• heuropt::traits::AlgorithmInfo — every built-in algorithm exposes its short canonical name ("Nsga3", …) and its seed. Used by the explorer JSON export.

If you don't want any of this, no migration needed — just bump the version.

To 0.8

From 0.5.x

Additive feature only. Bumping heuropt = "0.8" is enough for any code that doesn't need async evaluation. To opt into async, enable the new feature flag:

heuropt = { version = "0.8", features = ["async"] }

What changed:

• New async feature flag, gated on the futures crate.
• New core::async_problem::AsyncProblem trait — mirrors Problem but with async fn evaluate_async.
• New core::async_problem::AsyncPartialProblem trait — mirrors PartialProblem for multi-fidelity (Hyperband) workloads.
• run_async(&problem, concurrency).await on every algorithm in the catalog (33 of them) for IO-bound evaluations.
• New cookbook recipe: Async evaluation.

From 0.7.x

0.7.0 introduced an experimental observability layer (Snapshot, Observer, run_with, MaxTime, TargetFitness, Stagnation, Periodic, AnyOf, AllOf, TracingObserver) and three additional Pareto metrics (igd, igd_plus, r2). All of those were rolled back in 0.8.0 — the design didn't bake long enough and they shipped half-wired (run_with was overridden on only 3 of 35 algorithms). The tracing feature flag is also gone.

If your code uses any of those APIs, the migration is:

• Remove all run_with(&problem, &mut observer) calls and replace with run(&problem).
• Remove all uses of Observer, Snapshot, ControlFlow, MaxTime, MaxIterations, TargetFitness, Stagnation, Periodic, AnyOf, AllOf, TracingObserver.
• Remove all uses of metrics::igd::igd, metrics::igd::igd_plus, metrics::r2::r2.
• Remove Population::as_slice() calls (the method is gone).
• Drop the tracing feature from your Cargo.toml if you had it.

Stop conditions can still be implemented by wrapping run in a loop with a custom RNG-driven termination, or by wrapping the algorithm yourself; observers may return as a public API in a future release once the design has settled.

The async work introduced in 0.7.0 (AsyncProblem + run_async) survived and is broadened in 0.8: every algorithm in the catalog now has a run_async (0.7.0 only had it on three of them), and multi-fidelity problems get a parallel AsyncPartialProblem trait that Hyperband's run_async consumes. Existing call sites continue to work unchanged.

To 0.5

From 0.4.x

No public-API changes. v0.5 is a documentation-and-polish release. Bumping heuropt = "0.5" in your Cargo.toml is enough.

What changed:

• Added a comprehensive mdbook user guide (this book).
• Added runnable rustdoc examples on every public algorithm, operator, metric, and Pareto utility.
• Added real-world examples/portfolio.rs, examples/hyperparam_tuning.rs, and examples/scheduling.rs.
• Added CONTRIBUTING.md, SECURITY.md, CODE_OF_CONDUCT.md (Builder's Code of Conduct), GitHub issue templates, and PR template.

The full list is in CHANGELOG.md.

From earlier than 0.4

If you're coming from 0.3.x or earlier, also read the older sections below.

To 0.4

From 0.3.x

No public-API changes. v0.4 was a testing-infrastructure expansion + perf pass. Same cargo update story.

The compare-harness wall-clock got 3.27× faster on v0.4 with bit-identical quality metrics, so any benchmark numbers you have from v0.3 are still numerically accurate but will run faster.

To 0.3

From 0.2.x

Additive only. New algorithms (Bayesian Optimization, TPE, (1+1)-ES, IPOP-CMA-ES, sNES, Nelder-Mead, Hyperband), new operators (LevyMutation, ClampToBounds, ProjectToSimplex), new traits (PartialProblem, Repair<D>).

CmaEsConfig gained an initial_mean: Option<Vec<f64>> field; existing call sites need a .. CmaEsConfig { initial_mean: None, .. } update.

To 0.2

From 0.1.x

Additive. New algorithms across the catalog (Hill Climber, SA, GA, PSO, CMA-ES, Tabu Search, Ant Colony, UMDA, TLBO, MOPSO, IBEA, SMS-EMOA, HypE, RVEA, PESA-II, ε-MOEA, AGE-MOEA, GrEA, KnEA), new operators (SimulatedBinaryCrossover, PolynomialMutation, CompositeVariation, BoundedGaussianMutation), and the hypervolume_nd metric.

Optimizer<P> impls now require P: Sync and P::Decision: Send (this enables the parallel feature without changing the public trait surface). Any normal Problem you've written satisfies these bounds automatically.