Wierstra et al. 2008/2014 NES with the diagonal-covariance "separable"
variant (sNES). Different theoretical foundation from CMA-ES: rather
than tracking a full covariance matrix and adapting it through
evolution paths, sNES updates the sampling distribution's parameters
by following the natural gradient of expected fitness.
Each generation:
- Sample λ offspring from N(μ, diag(σ²))
- Rank-shape the fitnesses (utility weights from the standard NES table)
- Update μ along the natural gradient: μ ← μ + η_μ · σ · sum(u_i · z_i)
- Update σ multiplicatively: σ_j ← σ_j · exp(η_σ/2 · sum(u_i · (z_i,j² - 1)))
Vec<f64> decisions only, single-objective only. The diagonal covariance
makes per-step cost O(λ·n) instead of CMA-ES's O(λ·n²) — much faster on
high-dimensional problems where full-covariance tracking is expensive
or numerically fragile, at the cost of being unable to handle strongly
rotated landscapes.
The first sample-efficient algorithm in heuropt. Bayesian optimization
maintains a Gaussian-process surrogate of the objective and at each
step picks the next decision by maximizing an acquisition function on
that surrogate, so the evaluation budget is used surgically.
Implementation:
- **Kernel**: anisotropic RBF (squared-exponential) with per-axis
length scales, signal variance, and a small noise/jitter floor.
Hyperparameters are exposed in the config; a future version can add
marginal-likelihood maximization.
- **Posterior**: standard formulation. Cholesky factorizes K (using
the new internal helper); mean and variance predictions follow.
- **Acquisition**: Expected Improvement against the best observed
feasible point. Optimized by best-of-N random sampling — simple,
predictable cost, no inner-optimizer footgun.
- **Initial design**: `initial_samples` uniform-random points in
bounds before the BO loop starts.
- **Constraints**: feasibility-aware EI — best observed value uses
only feasible points; infeasible candidates are penalized.
Vec<f64> decisions, single-objective only. Targets the regime no
existing heuropt algorithm covers: 50–500 evaluations on an
expensive black-box function (CFD sim, ML training run, real-world
measurement).
Tests cover convergence on the 1-D sphere within a tight evaluation
budget (~30 evals get to f < 1e-6 — vs population-based methods
needing thousands), deterministic reruns, and panic on
multi-objective + dim mismatches.
Hand-rolled `A = L · L^T` factorization plus forward/backward triangular
solves, used by the upcoming Bayesian Optimization implementation for
the GP posterior. Same f64 row-major Vec<Vec<f64>> interface as the
existing Jacobi eigen helper so we don't pull in nalgebra for one
algorithm.
Returns Err on non-positive-definite input (a small jitter is the
typical caller-side fix). Tested against the standard 2x2 case, the
3x3 known-result case, A·x = b round-trip, and the SPD-failure case.
Auger & Hansen 2005 IPOP-CMA-ES: wraps the existing CmaEs in a restart
loop that doubles the population size and re-randomizes the mean
whenever a restart trigger fires. Specifically addresses the failure
mode we observed on Rastrigin (vanilla CMA-ES = 2.3 vs DE = 0).
Restart triggers:
- The whole budget for one inner CmaEs run finishes without improvement
- (More sophisticated triggers — eigenvalue collapse, condition-number
blow-up, sigma stagnation — are left for future versions; the
per-run budget trigger captures the bulk of the practical benefit)
Each restart:
- Doubles the population_size (Auger & Hansen 2005)
- Re-randomizes the initial mean to a fresh point in the bounds box
- Resets sigma to the user's initial value
Same Vec<f64> + single-objective constraints as CmaEs. The total
budget is divided across restarts; restart budget grows with
population. Tests verify it beats vanilla CMA-ES on Rastrigin.
Nelder & Mead 1965: gradient-free local optimizer that maintains a
simplex of n+1 points in n-D and at each iteration replaces the worst
vertex by one of {reflect, expand, outside-contract, inside-contract,
shrink} relative to the centroid of the rest. The five standard
coefficients (reflection α=1, expansion γ=2, contraction ρ=0.5,
shrinkage σ=0.5) are exposed in the config but default to canonical
values so users can leave them alone.
Single-objective only, Vec<f64> only, bounds enforced by clamping
each new vertex. Termination is purely iteration-count for v0.2;
"vertices have collapsed" stopping is a future enhancement.
Filling a real gap: heuropt had population-based local search
(SimulatedAnnealing, HillClimber) but no classical direct-search
algorithm. Excellent for low-dim smooth-ish problems where a
population is overkill.
Rechenberg 1973's elemental evolution strategy: one parent, one child
each generation, accept the child if it is no worse, and adapt the
mutation step size by tracking the success rate. If more than 1/5 of
recent moves were accepted the search is too cautious — multiply σ by
`step_increase` (typical 1.22). Below 1/5 — divide by the same factor.
At 1/5 — leave it alone. The success window has length `adaptation_period`.
Single-objective only. Vec<f64> only. Generic Gaussian step bounded by
the embedded `RealBounds`.
Why ship it: it's the smallest possible self-adapting evolution strategy
and a useful pedagogical / baseline endpoint. Pairs well as the budget
floor ("give me anything cheaper than CMA-ES").
Zhang, Tian & Jin 2015 KnEA: many-objective MOEA that biases survival
selection toward 'knee points' on the Pareto front — points where a
small improvement in one objective costs a large degradation in
another.
Each generation:
- NSGA-II-like loop with offspring + non_dominated_sort
- For the splitting front, identify knee points by perpendicular
distance from the hyperplane connecting the front's extreme points.
Members further from the hyperplane (= more 'kneeness') are preferred.
- Survival keeps every knee-tagged member; if room remains, fill from
remaining members by largest perpendicular distance.
Knee points are intuitively the most attractive points on a Pareto
front when no preference information is available. KnEA pushes the
search toward them at the cost of less uniform front coverage.
Yang, Li, Liu & Zheng 2013 GrEA: many-objective MOEA whose secondary
ranking is a grid-based diversity score instead of crowding distance
or reference vectors.
Each generation:
- NSGA-II-like loop with offspring + non_dominated_sort
- For the splitting front:
- Translate by ideal/nadir; partition objective space into a
(`grid_divisions` per axis) grid
- For every member compute three grid scores:
- GR (grid rank) = sum of grid coordinates (closer to ideal = lower)
- GCD (grid crowding distance) = #neighbors within 1 grid unit (in any axis)
- GCPD (grid coordinate point distance) = max coord - min coord
- Sort F_l ascending by GR, then by GCD, then by GCPD
- Take the top `n - already_selected` survivors
GrEA's grid-based niching is a different lens from NSGA-III's reference
points and RVEA's reference vectors — particularly effective on
non-convex fronts where reference-vector approaches struggle.
Panichella 2019 AGE-MOEA: a many-objective MOEA that *infers* the
front's geometry (its L_p shape, where p = 1 is linear, p = 2 is
spherical, p < 1 is convex etc.) from the current non-dominated set
and uses that estimate to drive both proximity and diversity in
survival selection.
Each generation:
- NSGA-II-like loop: random parent selection + variation + evaluation
- Combine + non_dominated_sort
- Fill front-by-front; for the splitting front:
- Translate by ideal point z*
- Find extreme points by ASF (same as NSGA-III) and intercepts
- Estimate the geometry parameter p by minimizing
\|f − ideal\|_p constancy on the extreme points
- Score every member by survival_score = (proximity_to_ideal) +
(1 / nearest-neighbor distance in the same L_p frame)
- Keep the top scorers
The geometry estimation is the novel contribution; with 3+ objectives
it produces fronts whose spread better matches the true shape than
NSGA-III's reference points (which assume a known geometry).
Rao 2011 TLBO: parameter-free single-objective optimizer for Vec<f64>.
The selling point — uniquely among the metaheuristics we ship — is that
it has NO algorithm-specific hyperparameters: no F, CR, w, c1, c2, σ,
mutation rate, etc. Just population_size and generations.
Each generation has two phases:
- **Teacher phase**: identify the best individual (the 'teacher'). For
every learner, compute a 'mean' learner and try replacing it with a
candidate moved toward the teacher by a random fraction, scaled by
the gap between teacher and (TF · mean), where TF ∈ {1, 2}.
- **Learner phase**: each learner picks a random partner and tries
moving toward the better one of the pair. Only successful moves are
kept.
Single-objective only, Vec<f64> only, bounds enforced via clamping.
Tests cover Sphere1D convergence, deterministic reruns, and panic on
multi-objective.
Lévy-flight perturbation: each variable receives a step drawn from a
heavy-tailed Lévy(α) distribution rather than a Normal. The result is
"mostly small steps with rare big jumps," which gives a more
exploratory mutation than Gaussian without abandoning local search.
Decision type: Vec<f64>, with optional bounds (clamped per-axis if
`bounds` is non-empty). The step is sampled via Mantegna's algorithm
which generates Lévy(α) by combining two Normal samples and taking
the right power, controlled by the tail exponent `alpha` (typical
1.5; 1 is heavy, 2 collapses to Normal).
This is the only genuinely-different mutation kernel from Cuckoo
Search and other Lévy-flight metaheuristics; ship it as a Variation
operator usable from any algorithm rather than as a separate
algorithm.
Replaces strict Pareto dominance with ε-dominance: A ε-dominates B when
`floor(A_i / ε) ≤ floor(B_i / ε)` for every objective and strictly
less in at least one (minimization frame). The result is a regular
discretization of objective space — at most one archive member per
ε-box — so the front spreads out automatically and the archive size
self-limits without truncation tricks.
Steady-state design: each generation samples one parent from the main
population and one from the ε-archive, applies variation, evaluates
the child, and offers it to both archives. Every member's
ε-coordinates and the box-tie rules are precomputed each insertion.
Tests: produces a front on Schaffer N.1 with reasonable spread,
deterministic reruns, panic on `epsilon[i] <= 0.0` and on
`epsilon.len() != objectives.len()`.
Corne, Jerram, Knowles & Oates 2001: divides objective space into a
hyperbox grid and uses per-box population counts to drive selection
toward sparsely-populated regions.
Each generation:
- Maintain an external archive of non-dominated members
- Build a hyperbox grid (`grid_divisions` per axis on the archive's
current axis ranges); count members per box
- Selection picks two parents by region-based tournament: choose two
random non-empty boxes and take a uniform-random member from the
one with fewer occupants
- Variation produces an offspring; insert into archive, dropping
dominated members and (if archive overflows) the most-crowded
occupant of the most-occupied box
Tests cover non-empty front on Schaffer N.1, deterministic reruns,
and panic on `archive_size == 0`.
Cheng, Jin, Olhofer & Sendhoff 2016 RVEA: many-objective MOEA built
around a fixed set of Das–Dennis reference vectors. Each generation:
- Generate offspring via random parent selection + variation +
evaluation
- Combine population + offspring; translate by ideal point z*
- Associate every member with the reference vector whose angle to
the translated objective vector is smallest
- For each occupied vector, keep the member with the smallest
Angle-Penalized Distance (APD) score; the rest are dropped
- APD = (1 + α(t)·θ_max·γ) · |f − z*| where γ is the angle to the
associated reference and α(t) = (t / t_max)^2 anneals the angle
penalty over the run
This produces well-spread fronts at high objective counts where
Pareto-rank methods (NSGA-II, SPEA2) lose discrimination.
Bader & Zitzler 2011: HypE estimates hypervolume contributions via
Monte Carlo sampling instead of computing them exactly. The point of
the trick is that exact hypervolume becomes prohibitively expensive
beyond ~5 objectives, while MC sampling stays cheap and accurate
enough at any dimension.
Each generation:
- Generate offspring via parent selection + variation + evaluation
- Combine, run non_dominated_sort, fill front-by-front
- For the splitting front, estimate each member's HV contribution
by drawing `n_samples` uniform points in the box [ideal, reference]
and counting how many points are dominated by *exactly* one front
member — that count, divided by n_samples and multiplied by the
box volume, is the member's expected unique HV contribution.
- Drop members one at a time from the splitting front by smallest
estimated contribution.
Public API matches the rest of the MO algorithms (Config + Optimizer).
The reference point is supplied in the config so the user controls
the integration domain. Tests cover non-empty front, deterministic
reruns, and panic on dim-mismatched reference.
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
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.
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.
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.
- 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.
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.
- 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).
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