Beume, Naujoks & Emmerich 2007: a steady-state MOEA that uses
hypervolume contribution as the secondary survival selection criterion.
Each generation:
- Generate ONE child via parent selection + variation + evaluation.
- Combine population + child, run non_dominated_sort.
- The discarded individual is the worst-front member with the
smallest hypervolume contribution (computed via the new
hypervolume_nd_from_evaluations helper).
Selection-quality is excellent at moderate objective counts (2–4) at
the cost of higher per-step compute (each survival selection requires
N+1 hypervolume evaluations of size ≤ N each). Best paired with a
tightly-bounded objective space — the user supplies a fixed reference
point in the config.
Tests: produces a non-empty front on Schaffer N.1, deterministic
reruns, panic on `population_size == 0`, panic on
`reference_point.len() != objectives.len()`.
Generalizes the existing 2-D hypervolume to arbitrary M ≥ 1 dimensions
using the standard recursive Hypervolume-by-Slicing-Objectives (HSO)
algorithm from While et al. 2006:
- For M = 1: return reference[0] - min(points[0])
- For M = 2: sort by axis 0, sweep accumulating rectangles (matches
hypervolume_2d's existing exact behavior)
- For M ≥ 3: sort by the last axis, peel off slices of increasing
thickness and recursively compute the (M−1)-dimensional HV of each
slice's projected non-dominated subset
Direction-aware: minimization-oriented input is the entry point, so
maximize objectives are negated by the caller via
`ObjectiveSpace::as_minimization` before the recursion runs.
Tested against:
- the existing 2-D analytical case (3 points → area 6)
- a known 3-D unit-cube case (1 point at origin, ref [1,1,1] → 1)
- empty front → 0
- agreement with hypervolume_2d on random 2-D fronts
A substantial README section walking newcomers through choosing an
optimizer. Defines the terminology as it comes up — single- vs multi-
vs many-objective, Pareto front, dominance, multimodality, evaluation
cost — so a reader who has never touched heuristic optimization can
still pick a sensible starting algorithm.
Five-step decision flow:
1. What is the decision?
2. How many objectives?
3. What's the landscape like? (multimodal, smooth, discrete)
4. How expensive is each evaluation?
5. Are there constraints?
Each branch ends with 1–3 algorithm recommendations and a one-line
rationale, plus a compact "quick reference" table at the bottom for
returning users.
Adds runners for HillClimber, SimulatedAnnealing, GeneticAlgorithm,
ParticleSwarm, CmaEs, and Umda to `examples/compare.rs`. Rastrigin
section now compares 8 single-objective optimizers against each other
on a fixed evaluation budget.
The MO sections (ZDT1, DTLZ2) are unchanged for now — MOPSO and IBEA
get added in a follow-up commit so each algorithm's debut shows up
clearly in the harness.
Mühlenbein 1997 UMDA: simplest Estimation-of-Distribution Algorithm for
`Vec<bool>` problems. Each generation:
- Evaluate the current population
- Select the top μ members by fitness
- Estimate per-bit marginal probability p_i = (count of 1s at bit i in
the μ-best) / μ
- Sample population_size new individuals from the resulting product-of-
Bernoullis distribution
Single-objective only. Bit-wise probabilities are clamped to
`[1 / (2 · μ), 1 - 1 / (2 · μ)]` to keep the population from collapsing
to a deterministic single string before convergence is meaningful
(standard Laplace-style smoothing for UMDA).
Tests: solves OneMax (maximize Σ bits) on a 20-bit instance,
deterministic reruns, panic on multi-objective.
Dorigo-style Ant System for permutation problems on a complete graph:
each generation, every ant constructs a tour by probabilistically
picking the next node from those it has not yet visited, weighted by
`τ_ij^α · η_ij^β` where τ is the pheromone level on edge (i, j) and
η is the heuristic desirability (1 / distance, here). After all ants
finish, pheromone evaporates by a factor `(1 - ρ)` and is reinforced
on each ant's tour proportional to that tour's quality.
Decision type is `Vec<usize>` (a permutation of 0..n_cities). The user
supplies a distance matrix and the n_cities is inferred. Single-objective
only (the cost is total tour length, which the Problem evaluates).
Tests build a 5-city ring and verify ACO finds a near-optimal tour,
plus deterministic reruns and panic on multi-objective.
Zitzler & Künzli 2004 IBEA: replaces Pareto-rank + crowding fitness
with a single scalar fitness derived from a binary quality indicator
(here, the additive ε-indicator). Loses no information at three or
more objectives the way crowding distance does.
Algorithm:
- For every (i, j) pair compute I(i, j) = max_k (f_k(i) - f_k(j)) on
minimization-oriented objectives.
- Fitness F(i) = -Σ_{j≠i} exp(-I(j, i) / κ).
- Each generation: combine parents + offspring, iteratively remove the
lowest-F member (cleanly recomputing the contribution of the dropped
member from each surviving member's fitness) until population_size
remain.
- Parent selection: binary tournament on F (higher wins).
Bounds-aware operators recommended (SBX + PolyMut).
Tests: produces a non-empty front on Schaffer N.1, deterministic
reruns, panic on `population_size == 0`.
Coello, Pulido & Lechuga 2004 MOPSO: PSO adapted for multi-objective
optimization via an external Pareto archive used as the source of
swarm leaders.
Each generation:
- Evaluate every particle's current position
- Insert non-dominated members into the archive (using ParetoArchive)
- For each particle, pick a leader from the archive (uniform random
among archive members)
- Update velocity using inertia + cognitive (toward pbest) + social
(toward leader)
- Update positions, clamp to bounds
- Refresh personal bests using Pareto comparison: pbest is replaced
only when the new position dominates it; on non-dominated, keep
with 50/50 random tiebreak
Vec<f64> decisions only. Truncates the archive to `archive_size` via
the existing simple-tail truncation. Tests: produces a non-empty
front on Schaffer N.1, deterministic reruns, panic on
single-objective.
Hansen & Ostermeier 2001 CMA-ES, the canonical real-valued
single-objective stochastic optimizer. Implements the full (μ/μ_w, λ)
update with rank-μ + rank-1 covariance updates and cumulative step-size
adaptation:
- Sample λ offspring from N(mean, σ² · C)
- Select the μ best, weight them, recompute mean
- Update evolution paths p_σ (step size) and p_c (covariance)
- Rank-1 update of C from p_c, plus rank-μ update from selected offspring
- Adapt σ via |p_σ| / E‖N(0,I)‖
Eigendecomposition (used to convert C into its B·D form for sampling
N(0, σ²·C)) goes through the new internal Jacobi helper, recomputed
every `eigen_decomposition_period` generations to amortize cost.
Vec<f64> decisions only. Bounds taken from a `RealBounds` field; mean
and offspring are clamped per dimension. Single-objective only.
Hyperparameters use the standard CMA-ES defaults (μ=λ/2, weights from
Hansen's tutorial, c_σ, c_c, c_1, c_μ, d_σ all formulae from §7.1).
Tests cover: convergence on Sphere1D and 5-D Rosenbrock, deterministic
reruns, panic on multi-objective, panic on `population_size < 4`.
Hand-rolled symmetric-matrix eigendecomposition via the cyclic Jacobi
rotation method. Returns sorted (eigenvalue, eigenvector) pairs in
descending order. Pure f64 row-major `Vec<Vec<f64>>` interface so we
don't pull in nalgebra for one algorithm.
Lives in `src/internal/eigen.rs` (new module). Used by the upcoming
CMA-ES implementation to maintain the covariance matrix's
eigendecomposition each generation. Tested against the standard
2x2 case, the diagonal case, and a known 3x3 result.
Glover 1986 tabu search for single-objective problems. Generic over
decision type — the user supplies a neighbor-generator closure that
produces a finite list of candidate moves from the current incumbent
(e.g. all 2-swaps for a permutation, or N Gaussian-perturbed copies of
a real vector). Each iteration picks the best non-tabu neighbor (with
an aspiration override that lets a tabu move through if it beats the
best-seen-ever incumbent) and adds the chosen move's decision to a
fixed-size FIFO tabu list.
Single-objective only. Tracks the best-seen-ever incumbent across the
run, returned as the result. Generic over the decision `D: Hash + Eq`
so the tabu list can match by full decision (simple and correct;
move-based tabu is left for users to implement themselves via a
custom decision wrapper).
Eberhart & Kennedy 1995 PSO with the standard inertia-weight update:
v[i,t+1] = w·v[i,t] + c1·r1·(pbest[i] - x[i,t]) + c2·r2·(gbest - x[i,t])
x[i,t+1] = clamp(x[i,t] + v[i,t+1], bounds)
Single-objective only, `Vec<f64>` decisions only (PSO's velocity vector
needs a Euclidean structure that doesn't generalize cleanly to bool/perm).
Velocities are clamped to ±(hi - lo) per dim to keep particles from
exploding off into space.
Config exposes the four standard knobs — swarm size, generations,
inertia w, cognitive c1, social c2 — plus a seed. Tests cover
convergence on Sphere1D, deterministic reruns, and panic on
multi-objective.
Canonical generational GA with elitism: each generation runs binary
tournament selection (using `tournament_select_single_objective`) on
the current population, applies the variation operator pair-wise to
produce offspring, evaluates them, then replaces the population while
preserving the top `elitism` members from the previous generation
(elitism prevents fitness regression on a single seed).
Single-objective only. Generic over decision type — pair with
`SimulatedBinaryCrossover + PolynomialMutation` for real-valued,
single-point crossover + bit-flip for binary, etc.
Tests: convergence on Sphere1D, deterministic reruns, panic on
multi-objective, panic on `population_size < 2`, panic on
`elitism > population_size`.
Classic Kirkpatrick et al. 1983 SA: hill climber that also accepts
worse moves with probability `exp(-Δ/T)` where T anneals geometrically
from `initial_temperature` to `final_temperature` over the iteration
count.
Single-objective only. Generic over decision type — works on real
vectors, bool vectors, permutations, anything. Tracks the best-seen
incumbent across the run (not just the last accepted move) so the
result reflects the actual best ever visited, not where the random
walk happened to end.
Tests cover: convergence on Sphere1D under reasonable hyperparameters,
deterministic reruns, panic on multi-objective, panic on
non-positive temperatures.
The simplest possible local search: start from one initializer-sampled
decision, repeatedly mutate it via the variation operator, and keep the
child only when it is strictly better than the current incumbent (with
the standard feasible-beats-infeasible / lower-violation tiebreaks
when relevant).
Single-objective only — panics with a clear message if the problem
exposes more than one objective. Deterministic under a seed. Returns
a population/front of size one (the current incumbent) so it slots
into the comparison harness like any other optimizer.
The Python tune_runtime.py treats the morning boot press and the 13:00
post-lunch re-tap as 'free' and only counts extra warning-phase taps.
That undercounts what the user actually presses each day and breaks
any comparison against a stated 'presses/day' comfort cap.
Updated `simulate_one` to count every press the user makes:
- boot press at workday start (always +1)
- 13:00 re-login press when the workday continues past lunch (+1)
- per-minute Bernoulli warning-phase presses (already counted)
- death-restart press: each time the device transitions running→dead
during workday and the user is at-desk (not at lunch), the user
presses to restart the cycle (warning press and death-restart for
the same cycle are mutually exclusive — extending via warning press
prevents that cycle's death)
With baseline now ~2 presses/day already mandatory, the hinge/cap
shift up too: PRESS_HINGE_LOW = 2.5/d (full reward up to baseline +
half a warning press) and PRESS_COMFORT_CAP = 3.5/d (rejected above).
Output 'Why' bullet now reports the total directly and notes the
component breakdown so the number is interpretable against the
new thresholds.
Tweak the a-posteriori scoring to match the user's stated preferences:
- Reweight: lunch_sleep 30%, after_hours 25%, work_fail 20%,
presses 15% (with hinge below), balance 10%.
- Press term is now a hinge instead of a normalized minimize:
* <= 2 presses/day → score 1.0 (no penalty)
* 2 → 3 presses/day → linear ramp from 1.0 to 0.0
* > 3 presses/day → -inf (excluded; comfort cap)
- New balance term: bonus for longer warning phases. Computed as
min(YA - RA, RA - FRA), saturated at 10 minutes. So a 5/5/X split
scores 0.5, an 8/8/X split scores 0.8, and 10/10/X or wider saturates
at 1.0.
Constants moved to module scope so the printout in main and the
scoring function stay in sync.
Adds an a-posteriori decision step to the jiggly example. After NSGA-III
produces the Pareto front, we apply a weighted-sum score over each
objective normalized to [0, 1] across the front (best→1, worst→0,
direction-aware), and report the top three plus a clear recommendation.
The weights are stated explicitly with rationale, not buried in code:
work_fail 45% — screen sleeping mid-meeting is the worst failure
lunch_sleep 30% — the actual design goal
presses 15% — UX friction the user feels
after_hours 10% — minor, mostly screen burn
This is the standard structure for picking a single answer out of a
Pareto set without losing the front itself: someone with different
weights can read the front and pick differently, but we surface a
specific recommendation with reasoning rather than leaving the user
to stare at 84 incomparable rows. Identical normalization could be
swapped for TOPSIS or knee-point detection later if useful.
Port of `scripts/tune_runtime.py` from ~/Code/jiggly: optimize the four
lifecycle constants of a USB mouse-jiggler firmware so the screen sleeps
during the user's lunch hour rather than failing during work.
The Python script grid-searches against a single composite score that
linearly combines several genuinely conflicting goals — a workaround
for the fact that grid search needs one number to rank by. heuropt has
the actual right tool, so this example is structured as a 4-objective
NSGA-III run that surfaces the Pareto front of legitimate tradeoffs:
1. minimize work-time failures (mean_work_sleep)
2. maximize lunch sleep (mean_lunch)
3. minimize human button presses (mean_presses)
4. minimize after-hours waste (mean_after)
Decision: 4-element `Vec<f64>` for (RT, YA, RA, FRA), continuous-relaxed
and rounded to integer minutes inside `evaluate`. The firmware ordering
constraint YA > RA > FRA > 0 is encoded as `constraint_violation` so
heuropt's feasible-beats-infeasible logic handles it for free.
Solver: NSGA-III with M=4, H=6 → 84 reference points, matching the
population size. Each evaluate runs a 1,000-workday Monte Carlo, so
the example is also a deliberately meaty evaluator that benefits from
`--features parallel`.
Output is in jiggly's native units — RT as Xh00m, thresholds in plain
minutes, sleep durations as Xh00m / Mm, probabilities as percentages —
and contrasts the Pareto front against:
- the four extreme single-axis winners (most lunch / fewest work fails /
fewest presses / least after-hours)
- the firmware's currently-shipping defaults (which sit inside the
front as a balanced compromise)
Two new runners — `zdt1_moead` and `dtlz2_moead` — using the same
SBX + PolyMut variation as the other Pareto-based methods. Reference
divisions chosen so the implied population size is comparable to the
other algorithms in each section (99 → 100 weights for ZDT1; 12 → 91
weights for DTLZ2).
Implementation of Zhang & Li 2007 MOEA/D — the canonical
decomposition-based MOEA. Different paradigm from Pareto-dominance
algorithms: each subproblem is a scalarized single-objective problem
defined by a Das–Dennis weight vector, and subproblems with similar
weight vectors form neighborhoods that share genetic material.
Each generation iterates over every weight vector `i`:
1. Pick two parents uniformly from the T-nearest neighbors of weight i
(T = neighborhood_size).
2. Apply variation, evaluate the child.
3. Update the ideal point z* with the child's objectives.
4. Walk the entire neighborhood: for each j, if the child's
Tchebycheff value g(child | w_j, z*) <= g(current[j] | w_j, z*),
replace current[j] with the child.
Tchebycheff scalarization:
g(f | w, z*) = max_k w_k · |f_k - z*_k|
(With the standard `w_k = 1e-6` floor when a weight is zero, so the
max well-defined.)
Public API:
MoeadConfig {
generations,
reference_divisions, // Das-Dennis H, also fixes population size
neighborhood_size, // T
seed,
}
Moead { config, initializer, variation }
impl<P, I, V> Optimizer<P> for Moead<I, V>
Population size equals the number of weight vectors generated by
das_dennis(num_objectives, reference_divisions). Re-exported from the
prelude. Tests cover non-empty Pareto front, deterministic reruns,
and panic on `reference_divisions` that would yield zero weights.
- nsga3: drop redundant `.into_iter()` in extend call; use
`#[allow(clippy::needless_range_loop)]` on the back-substitution
loop where `j` indexes into the matrix; remove an unneeded
`return` keyword in a closure.
- spea2: switch `pool.extend(x.drain(..))` to `pool.append(&mut x)`.
- examples/compare.rs DTLZ2 evaluator: same `needless_range_loop`
silencer on the inner cosine product loop.
NSGA-III's value over NSGA-II shows up at 3+ objectives, where
crowding distance loses its diversity signal. Adds a third comparison
section to `examples/compare.rs`:
DTLZ2 (3-objective, 12-D, the textbook benchmark for many-objective
algorithms): unit-sphere-octant Pareto front. Compares RandomSearch,
NSGA-II, SPEA2, and NSGA-III on:
- mean distance from front points to the unit sphere
(closed-form: |1 - sqrt(f1² + f2² + f3²)|),
- spacing,
- front size,
- wall-clock ms.
NSGA-III config: H=12 reference divisions (91 reference points,
matching the canonical setup from Deb & Jain 2014).
Also wires NSGA-III into the existing ZDT1 (2-objective) section even
though it's not its sweet spot — useful as a regression check that the
algorithm at least keeps up with NSGA-II on bi-objective problems.
Implementation of Deb & Jain 2014 NSGA-III — the canonical
many-objective MOEA. Replaces NSGA-II's crowding-distance niching
with a structured reference-point niching procedure that scales to
3+ objectives where crowding distance loses its diversity signal.
Each generation:
1. Random parent selection + variation + offspring evaluation, same as
NSGA-II.
2. Combine + non_dominated_sort, fill the next population front-by-
front until the next front would overflow (the splitting front F_l).
3. Survival on F_l uses reference-point niching:
- Translate by the ideal point z* (per-axis min in oriented space).
- Compute extreme points by ASF and intercepts; normalize by
intercepts (with a robust fallback to per-axis range if extreme
points are degenerate).
- Associate every member of the working pool with the closest
reference direction by perpendicular distance.
- Iteratively pick from F_l: prefer the niche with the smallest
count among references that have F_l candidates; if the niche is
empty in the already-selected set, take the closest associated
member by perpendicular distance, otherwise pick uniformly from
the niche.
Public API:
Nsga3Config { population_size, generations, reference_divisions, seed }
Nsga3 { config, initializer, variation }
impl<P, I, V> Optimizer<P> for Nsga3<I, V>
Re-exported from the prelude. Tests cover non-empty Pareto front,
exact final population size, deterministic reruns, and panic on
`population_size == 0`. Uses the existing tests_support problems.
The standard structured weight/reference vector generator for
many-objective MOEAs (NSGA-III, MOEA/D). Generates (H+M-1 choose M-1)
points uniformly distributed on the unit simplex by enumerating all
integer compositions of `divisions` into `num_objectives` parts and
dividing each by `divisions`.
Lives in src/pareto/reference_points.rs. Re-exported from the prelude
as `das_dennis`.
Tests cover: M=2/H=4 → 5 points along the diagonal; M=3/H=12 → 91
points (the canonical NSGA-III 3-objective ref set); each generated
point has exactly M components summing to 1 within float tolerance.
Implementation of Zitzler, Laumanns, Thiele 2001 SPEA2 — the classic
Pareto MOEA built around an explicit external archive of fixed size.
Each generation:
1. Combine the current population and the archive into one pool.
2. For every member, compute strength S(i) = number of others that
member dominates, then raw fitness R(i) = sum of S(j) over members
j that dominate i.
3. Add a density estimator D(i) = 1/(σ_k + 2) where σ_k is the distance
to the k-th nearest neighbor (k = floor(sqrt(|pool|))) in
minimization-oriented objective space.
4. Final fitness F(i) = R(i) + D(i); lower is better.
5. Build the next archive by taking every non-dominated member
(R(i) == 0). If too many, prune by repeatedly removing the member
with the smallest k-th-nearest-neighbor distance. If too few, fill
from the rest sorted by F ascending.
6. Generate the next population by binary tournament on F (lower wins),
then variation, then evaluation.
Public API mirrors the other algorithms:
Spea2Config { population_size, archive_size, generations, seed }
Spea2 { config, initializer, variation }
impl<P, I, V> Optimizer<P> for Spea2<I, V>
Re-exported from the prelude. Tests cover archive size invariants,
non-empty Pareto front on Schaffer N.1, deterministic reruns under
the same seed, and panic on population_size == 0.
A comparison example that runs every applicable optimizer on ZDT1 and
Rastrigin across N seeds and reports mean ± stddev for each quality
metric. Designed so a new algorithm slots in by adding a single runner
function — no harness changes needed.
ZDT1 (multi-objective, dim=30):
Reports hypervolume_2d (against ref point [1.1, 1.1]), spacing, mean
L2 distance to the analytical Pareto front, front size, and wall-clock
ms. RandomSearch, PAES, and NSGA-II all use bounds-aware operators
(RealBounds, BoundedGaussianMutation, SBX+PolyMut) so the Problem
itself stays unclamped — apples-to-apples.
Rastrigin (single-objective, dim=5):
Reports mean ± stddev best objective and ms. RandomSearch, PAES,
NSGA-II (degenerate single-obj case), and DE.
Default budget: 10 seeds × 25,000 evaluations on ZDT1, × 50,000 on
Rastrigin. Run with:
cargo run --release --example compare
- PolynomialMutation::vary: `#[allow(clippy::needless_range_loop)]`
on the per-dimension loop — body indexes both `self.bounds[j]` and
`child[j]` so a range index is the cleanest option.
- Operator tests: replace `x >= lo && x <= hi` with
`(lo..=hi).contains(&x)` per clippy's manual_range_contains lint.
Replace the v0.1 `GaussianMutation` + clamp-inside-`evaluate` setup
with the canonical NSGA-II operator pair: SBX (η_c=15, per-var prob 0.5)
followed by PolynomialMutation (η_m=20, per-var prob 1/dim), composed
via `CompositeVariation`. Both are bounds-aware on their own, so the
in-evaluate clamping is dropped.
Result on ZDT1 (dim=30, pop=100, gens=1000, seed=42): mean L2 distance
to the analytical Pareto front is 0.00152 — comfortably within the
published NSGA-II range for this benchmark.
Note on the previous number: the v0.1 setup reported 0.00072 at 40k
evals, but that was an artifact of clamping inside `evaluate`. Out-of-
bounds Gaussian mutations on `x[0]` were snapping to 0, which
coincides with the ZDT1 Pareto-front extreme (f1=0). The new operator
pair has no such free lunch — it runs the actual NSGA-II algorithm —
and the new measurement is what honest convergence on ZDT1 actually
looks like.
Generations bumped from 400 to 1000 (40k → 100k evaluations) to give
the operators headroom; matches the budget DE uses for Rastrigin so
the example feels balanced.
Generic two-stage Variation operator: runs an inner crossover-style
operator on the parents, then applies an inner mutation-style operator
to each resulting child. Lets users build the canonical NSGA-II
operator stack — `SimulatedBinaryCrossover` followed by
`PolynomialMutation` — by composing the existing primitives instead
of bundling a one-off SbxPolyMut struct.
Lives in src/operators/composite.rs to keep type-specific operator
files unchanged. Generic over decision type and over both inner
operators.
Deb's standard real-valued mutation pair to SBX, used together by
canonical NSGA-II. For each variable, with probability
`per_variable_probability` (typical: 1/n where n is dim), perturb the
parent value by a polynomial-distributed delta scaled by the bound
range, then clamp.
Per-dim formula:
- `u ~ U[0, 1)`
- `δ = (2u)^(1/(η+1)) − 1` if `u < 0.5` else `1 − (2(1−u))^(1/(η+1))`
- `child[j] = parent[j] + δ · (hi − lo)`, clamped to bounds
`eta` is the distribution index (typical 20; smaller → more spread).
This is the simple bound-rescale form; the bound-aware δ_q variant from
the full paper is left as a future refinement.
Always returns one child. Tests cover: child stays in bounds with high
sigma-equivalent eta, per_variable_probability=0 returns the parent
unchanged, and standard panics.
Deb & Agrawal's standard real-valued crossover for NSGA-II. Takes two
parents, returns two children; per dimension, with
`per_variable_probability`, mixes the parents using a polynomial
spread parameter \\(\\beta\\) drawn from a distribution controlled by
`eta` (the distribution index — typical values 10–30, default 15).
Children are clamped to per-variable bounds.
Per-dim formula (Deb & Agrawal 1995):
- `u ~ U[0, 1)`
- `β = (2u)^(1/(η+1))` if `u ≤ 0.5` else `(1 / (2(1-u)))^(1/(η+1))`
- `c1 = 0.5·((1+β)·p1 + (1-β)·p2)`, `c2 = 0.5·((1-β)·p1 + (1+β)·p2)`
This is the simple compute-then-clamp form; the bounds-aware
β formulation from the full paper is left as a future refinement.
Tests cover: two children for two parents, output lengths preserved,
all variables clamped to bounds, and per_variable_probability=0
returns the parents unchanged.
A bounded variant of GaussianMutation: same Gaussian noise applied to
the first parent, but every variable is clamped to its per-dimension
inclusive bound. Useful as a drop-in for problems that need feasibility
maintained across generations rather than relying on
clamp-inside-evaluate.
Panics on `sigma <= 0.0`, on no parents, and on construction if any
`(lo, hi)` has `lo > hi`. Decision length must match the bounds
length when called.
Adds a `parallel` Cargo feature that pulls in rayon and parallelizes
the only step that's actually expensive in practice — calls to
`Problem::evaluate` — across the population. RNG-driven steps (parent
and donor selection, variation, replacement decisions) stay serial, so
seeded runs remain deterministic regardless of feature state, and the
default and `--features parallel` builds produce bit-identical
results.
Wiring:
- New `algorithms::parallel_eval::evaluate_batch` helper with two
cfg-gated implementations (rayon's `into_par_iter` when the feature
is on, plain `into_iter` otherwise). Both preserve input order, so
pareto_front and crowding-distance decisions remain reproducible.
- `RandomSearch`, `Nsga2`, and `DifferentialEvolution` now route
population/offspring evaluation through the helper. NSGA-II's main
loop is restructured into a serial selection-and-variation phase
followed by a parallel-friendly batch evaluation phase.
- DE's per-target loop is restructured into three phases (serial trial
construction → batch evaluation → serial replacement). Side effect
of the restructuring: DE is now the canonical synchronous DE/rand/1/bin
rather than the asynchronous variant where target `i+1` sees `i`'s
in-flight update. Synchronous is the textbook formulation, so this
is a small correctness improvement on top of the parallelism enable.
- PAES stays serial — its main loop has a sequential dependency on the
current candidate and would gain nothing from rayon.
Cost: algorithm impls now require `P: Sync` and `P::Decision: Send`
unconditionally so a single impl serves both feature modes. This is a
small bound tightening that any plain-data Problem already satisfies; in
return the public `Problem` trait itself stays unchanged and the
default build picks up no new dependencies.
Verified:
- `cargo test` and `cargo test --features parallel` both pass; the
Nsga2 `deterministic_with_same_seed` test confirms reproducibility.
- `cargo run --release --example benchmarks` and the same with
`--features parallel` produce bit-identical ZDT1 / Rastrigin
results.
Two canonical optimization benchmarks in a single runnable example:
- ZDT1 (Zitzler-Deb-Thiele 1): 30-D, two minimization objectives,
closed-form Pareto front \\(f_2 = 1 - \\sqrt{f_1}\\) for
\\(f_1 \\in [0, 1]\\). Solved with NSGA-II.
- Rastrigin: highly multimodal single-objective, global minimum
\\(f = 0\\) at the origin. Solved with DE.
Both are public-domain mathematical formulas. Implemented as Problem
impls in examples/benchmarks.rs; main() runs each, prints front /
best, and (for ZDT1) reports the mean L2 distance from the known
analytical Pareto front so the example doubles as a sanity check on
solution quality.
- pareto/crowding.rs: rewrite the inner loop to iterate per-objective
via index_axis-style indexing on `oriented` rather than naming an
unused loop variable `k`.
- operators/{binary,permutation}.rs tests: pass parents via
`std::slice::from_ref` instead of `&[parent.clone()]` to avoid the
cloned_ref_to_slice_refs lint.
Pure cleanup — no behavior change, all 83 unit tests + 2 doctests still
pass.
Adds:
- README.md following spec §19.1 (what / install / define problem /
run NSGA-II / custom optimizer / current algorithms / design
philosophy).
- A short-but-runnable crate-level //! example in lib.rs for
`cargo doc` (spec §19.2).
The three runnable examples called out in spec §18.5 / §19. All open
with `use heuropt::prelude::*;` so they double as a check that the
prelude is sufficient on its own:
- toy_nsga2.rs: Schaffer N.1 solved with NSGA-II.
- random_search.rs: 2D sphere solved with RandomSearch.
- custom_optimizer.rs: a minimal hill-climber implementing
`Optimizer<P>` directly, demonstrating spec §2.3.
Exact 2D dominated hypervolume against a fixed reference point. Sorts
points by the first minimization-oriented objective ascending, then
sweeps and accumulates the dominated rectangle area against the
reference. Points that don't strictly dominate the reference are
ignored. Panics with a clear message if the objective space does not
have exactly two objectives (spec §14.2).
Tests cover a known-area front, the no-coverage case, and the panic on
non-2D problems.
Standard Schott spacing: for each front point compute the Manhattan
distance to its nearest neighbor on minimization-oriented objective
values; the spacing metric is the population standard deviation of
those nearest-neighbor distances.
Returns 0.0 for empty or single-point fronts (spec §14.1).
Optional v1 algorithm requested by the user (spec §12.4):
- Vec<f64> decisions only.
- Single-objective only — panics with a clear message otherwise.
- Standard DE/rand/1/bin: for each target i, sample distinct r1, r2, r3;
mutant = x[r1] + F * (x[r2] - x[r3]); apply binomial crossover with at
least one forced index; greedy replacement on direction-correct
comparison.
- Bounds taken from the embedded RealBounds (mutants are clamped to the
per-variable range so the trial vector stays feasible).
- Seed-deterministic; tests verify reproducibility, that DE improves on
the initial random population for a sphere problem, and that
multi-objective use panics.
Standard (μ+λ) NSGA-II with binary tournament parent selection on
(rank, crowding distance) and elitist survival selection on the combined
parent + offspring population (spec §12.3):
1. Initialize population_size random decisions.
2. Each generation: select parents by binary tournament (rank ↑ then
crowding ↓ then random), apply variation, evaluate offspring,
combine, non_dominated_sort, fill the next population front-by-front
trimming the partial last front by crowding distance descending.
3. Return final population, Pareto front, best (None for >1 objective),
evaluation count, and generation count.
Internal Nsga2Entry { candidate, rank, crowding_distance } stays
private. Panics with clear messages on `population_size == 0` or
empty `vary` output. Tests cover population length, evaluation count,
non-empty front, and full determinism with the same seed (spec §18.4).
A readable v1 PAES (spec §12.2):
- Single starting decision from the initializer.
- Each iteration mutates the current decision via the Variation operator,
evaluates the child, and pareto_compares to the current.
- Dominating children become current; for non-dominated comparisons we
move to the child (acceptable v1 behavior per spec).
- Both current and child are inserted into a ParetoArchive truncated
to `archive_size` (simple tail-truncation in v1).
The final result returns the archive as both `population` and
`pareto_front`. Tests verify the archive never exceeds
`archive_size`.
The reference baseline and the spec's recommended starting example. Per
iteration it asks the initializer for `batch_size` decisions, evaluates
each, and accumulates them. At the end it returns the full population
plus the Pareto front and (if single-objective) the best feasible
candidate. `generations` equals `iterations`; `evaluations` equals
`iterations * batch_size` (spec §12.1).
Includes a tiny single-objective sphere test problem under
`tests_support` that later algorithm tests will reuse.
`select_random` samples `count` decisions with replacement and clones
them out of the population (spec §10.1).
`tournament_select_single_objective` runs binary-or-larger tournaments
with the spec's tiebreak order: feasible beats infeasible, lower
violation among infeasibles, and direction-correct objective comparison
among feasibles. Panics if not exactly one objective (spec §10.2).
Selection helpers stay under `heuropt::selection` and are not part of
the prelude (spec §15).
Variation that clones the first parent (a Vec<usize> permutation) and
swaps two distinct random indices when len >= 2 (spec §11.4). Tests
confirm the multiset of contents is preserved.
Variation that clones the first parent and flips each bit independently
with probability `probability`. Panics if probability is outside [0, 1]
(spec §11.3).
Tests verify that probability=0 produces an unchanged child and
probability=1 flips every bit (spec §18.3).
swaits