From 2b453464e78fa454bbe2b23b949bae9d09ef562e Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Thu, 14 May 2026 07:13:05 -0600 Subject: [PATCH] refactor(compare): realign + sort tables, add combinatorial problems MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The terminal output was misaligned: headers and data were right-aligned with hardcoded column widths and separator lengths, and the mean ± std cells contained the non-ASCII `±` (plus `ε`, `↑`, `↓`) -- on any terminal that renders those at a non-1 column width the columns drift, and the hardcoded `-`.repeat(n) separators didn't match the real table width anyway. Changes: - New `print_table` helper: column widths derived from the actual cell contents (header + every row), separator length computed to match. - All table cells are now ASCII: `+/-` instead of `±`, `eps-MOEA` instead of `ε-MOEA`, arrows dropped from headers. - Every table is sorted best-first by its primary quality metric. - Added three combinatorial / sequencing problems with their own (permutation- / bitstring-native) algorithm rosters: a convex-position ring TSP (known optimum), FT06 job-shop makespan (known optimum 55), and a bi-objective 0/1 knapsack scored by hypervolume. - Expanded every problem's preamble: what it is, why it's hard, and the best-known / optimal result. - Regenerated examples/compare-results.md to match. Continuous-problem quality metrics are unchanged (bit-identical to prior snapshots); only ms columns and row order move. Co-Authored-By: Claude Opus 4.7 (1M context) --- examples/compare-results.md | 276 ++++++---- examples/compare.rs | 1032 ++++++++++++++++++++++++++++++++--- 2 files changed, 1130 insertions(+), 178 deletions(-) diff --git a/examples/compare-results.md b/examples/compare-results.md index 68c8720..9d2a5e9 100644 --- a/examples/compare-results.md +++ b/examples/compare-results.md @@ -1,142 +1,234 @@ # `compare` example — reference output -Snapshot from `cargo run --release --example compare` after the v0.4.0 -perf pass landed (2026-05-05). 10 seeds per algorithm per problem. +Snapshot from `cargo run --release --example compare`, refreshed 2026-05-14 +for heuropt v0.10.0. 10 seeds per algorithm per problem. -The **quality metrics** (hypervolume / spacing / mean L2 / mean dist / -front size) are bit-identical to the v0.3.0 snapshot — the v0.4.0 -optimization work was strictly CPU-time, never algorithmic. The **ms -columns** reflect the v0.4.0 numbers; total compare-harness wall-clock -dropped from ~18.6 s to ~5.7 s (3.27× faster). +Each table is **sorted best-first** by its primary quality metric. The +live terminal output uses ASCII `+/-` for the mean ± std cells (so column +alignment can't be broken by a terminal that renders `±` at an odd +width); this doc uses `±` since markdown renders it fine. -Wall-clock numbers are from the development machine and will vary; -the *relative* numbers across algorithms are the interesting part. +The **continuous-problem quality metrics** are bit-identical to the +v0.3.0–v0.4.0 snapshots — every optimization pass so far (including the +Phase B CPU work) has been verified bit-identical by the `run()` snapshot +tests. The **ms columns** are the post-Phase-B numbers; SMS-EMOA on DTLZ2 +in particular fell ~2.7× from the `hypervolume_nd` rework. + +This refresh also adds three **combinatorial / sequencing** problems — +TSP, job-shop scheduling, and a bi-objective knapsack — which exercise the +permutation and bitstring operators and a different algorithm roster (the +real-vector methods can't run them). + +Wall-clock numbers are from the development machine and will vary; the +*relative* numbers across algorithms are the interesting part. --- ## ZDT1 (dim=30, 25000 evals/run × 10 seeds) -Two-objective benchmark with a smooth Pareto front along -`f₂ = 1 − √f₁`. Hypervolume reference point: `[11, 11]`. +Zitzler-Deb-Thiele 2-objective benchmark: 30 real variables, one smooth +convex Pareto front `f₂ = 1 − √f₁`. Hard because 29 of 30 variables must +collapse to 0 before the front is even reachable, and only then can the +population spread along it. Optimum: mean L2 → 0 (the front is known +exactly). Sorted by hypervolume (reference `[11, 11]`). -| algorithm | hypervolume ↑ | spacing ↓ | mean L2 ↓ | front | ms | +| algorithm | hypervolume ↑ | spacing ↓ | mean L2 ↓ | front | ms | |---|---|---|---|---|---| -| RandomSearch | 99.5691 ± 0.94 | 0.0937 ± 0.03 | 2.3621 ± 0.14 | 28 | 94 | -| PAES | 104.1887 ± 0.90 | 0.0351 ± 0.01 | 1.3195 ± 0.06 | 33 | 30 | -| MOPSO | **120.6149 ± 0.05** | 0.0125 ± 0.00 | **0.0005 ± 0.00** | 100 | 89 | -| SPEA2 | 118.0823 ± 0.60 | 0.0111 ± 0.00 | 0.2408 ± 0.05 | 97 | 234 | -| PESA-II | 119.3670 ± 0.33 | **0.0095 ± 0.00** | 0.0802 ± 0.04 | 100 | 73 | -| ε-MOEA | 118.8742 ± 0.68 | 0.0167 ± 0.01 | 0.0493 ± 0.02 | 45 | 50 | -| IBEA | 120.0167 ± 0.31 | 0.0130 ± 0.00 | 0.0448 ± 0.02 | 73 | 138 | -| HypE | 105.6489 ± 0.98 | 0.0266 ± 0.01 | 1.4820 ± 0.10 | 72 | 38 | -| SMS-EMOA | 102.8871 ± 1.05 | 0.0263 ± 0.00 | 1.4937 ± 0.12 | 40 | 67 | -| RVEA | 111.7151 ± 1.82 | 0.0308 ± 0.01 | 0.8399 ± 0.16 | 47 | 65 | -| NSGA-II | 118.3336 ± 0.78 | 0.0112 ± 0.00 | 0.1891 ± 0.06 | 96 | 67 | -| NSGA-III | 115.1612 ± 0.47 | 0.0139 ± 0.00 | 0.4314 ± 0.06 | 86 | 70 | -| MOEA/D | 119.9450 ± 0.50 | 0.0118 ± 0.00 | 0.0065 ± 0.00 | 96 | 28 | +| MOPSO | **120.6149 ± 0.0529** | 0.0125 ± 0.0025 | **0.0005 ± 0.0001** | 100 | 80 | +| IBEA | 120.0167 ± 0.3112 | 0.0130 ± 0.0027 | 0.0448 ± 0.0168 | 73 | 130 | +| MOEA/D | 119.9450 ± 0.4953 | 0.0118 ± 0.0013 | 0.0065 ± 0.0020 | 96 | 27 | +| PESA-II | 119.3670 ± 0.3261 | **0.0095 ± 0.0011** | 0.0802 ± 0.0354 | 100 | 67 | +| eps-MOEA | 118.8742 ± 0.6835 | 0.0167 ± 0.0058 | 0.0493 ± 0.0227 | 45 | 46 | +| NSGA-II | 118.3336 ± 0.7750 | 0.0112 ± 0.0022 | 0.1891 ± 0.0599 | 96 | 40 | +| SPEA2 | 118.0823 ± 0.5973 | 0.0111 ± 0.0023 | 0.2408 ± 0.0509 | 97 | 226 | +| NSGA-III | 115.1612 ± 0.4745 | 0.0139 ± 0.0029 | 0.4314 ± 0.0582 | 86 | 47 | +| RVEA | 111.7151 ± 1.8195 | 0.0308 ± 0.0099 | 0.8399 ± 0.1569 | 47 | 62 | +| HypE | 105.6489 ± 0.9789 | 0.0266 ± 0.0053 | 1.4820 ± 0.1003 | 72 | 30 | +| PAES | 104.1887 ± 0.8953 | 0.0351 ± 0.0067 | 1.3195 ± 0.0558 | 33 | 27 | +| SMS-EMOA | 102.8871 ± 1.0543 | 0.0263 ± 0.0039 | 1.4937 ± 0.1192 | 40 | 54 | +| RandomSearch | 99.5691 ± 0.9383 | 0.0937 ± 0.0347 | 2.3621 ± 0.1428 | 28 | 88 | **MOPSO and MOEA/D dominate** convergence (mean L2 to true front ≤ 0.01). PESA-II edges spacing. ## ZDT3 (dim=30, 25000 evals × 10 seeds) -Disconnected Pareto front; tests an algorithm's ability to maintain -spread across gaps. +Zitzler-Deb-Thiele 2-objective with a **disconnected** front: five +separate arcs rather than one curve. Hard because an algorithm has to +discover and populate every arc while not stranding solutions in the +dominated gaps between them. -| algorithm | hypervolume ↑ | spacing ↓ | front | ms | +| algorithm | hypervolume ↑ | spacing ↓ | front | ms | |---|---|---|---|---| -| NSGA-II | 123.1826 ± 1.58 | 0.0092 ± 0.00 | 98 | 68 | -| MOEA/D | 125.2413 ± 2.16 | 0.0198 ± 0.00 | 92 | 28 | -| **IBEA** | **126.2072 ± 1.23** | 0.0164 ± 0.00 | 48 | 135 | -| AGE-MOEA | 119.5132 ± 1.27 | 0.0136 ± 0.00 | 90 | 199 | +| **IBEA** | **126.2072 ± 1.2280** | 0.0164 ± 0.0036 | 48 | 131 | +| MOEA/D | 125.2413 ± 2.1647 | 0.0198 ± 0.0043 | 92 | 27 | +| NSGA-II | 123.1826 ± 1.5829 | **0.0092 ± 0.0020** | 98 | 40 | +| AGE-MOEA | 119.5132 ± 1.2732 | 0.0136 ± 0.0023 | 90 | 171 | -## DTLZ2 (3-obj, dim=12, 30000 evals × 10 seeds) +## DTLZ2 (3-obj, dim=12, 30000 evals/run × 10 seeds) -Spherical Pareto front. Mean dist = `|‖f‖ − 1|`. +Deb-Thiele-Laumanns-Zitzler 3-objective; the Pareto front is the +unit-sphere octant (`Σf² = 1, all f ≥ 0`) — a curved 2-D surface embedded +in 3-D objective space. `mean dist = |‖f‖ − 1|`, so 0 means perfectly on +the sphere (the known optimum). -| algorithm | mean dist ↓ | spacing ↓ | front | ms | +| algorithm | mean dist ↓ | spacing ↓ | front | ms | |---|---|---|---|---| -| RandomSearch | 0.3949 ± 0.02 | 0.0797 ± 0.01 | 239 | 520 | -| MOPSO | 0.0566 ± 0.00 | 0.0687 ± 0.01 | 100 | 71 | -| NSGA-II | 0.0332 ± 0.01 | 0.0577 ± 0.01 | 92 | 104 | -| SPEA2 | 0.0368 ± 0.00 | **0.0288 ± 0.00** | 92 | 534 | -| PESA-II | 0.0395 ± 0.00 | 0.0616 ± 0.01 | 100 | 396 | -| ε-MOEA | 0.0325 ± 0.01 | 0.0572 ± 0.02 | 136 | 89 | -| **IBEA** | **0.0014 ± 0.00** | 0.0607 ± 0.00 | 87 | 156 | -| HypE | 0.0113 ± 0.00 | 0.0269 ± 0.02 | 80 | 53 | -| SMS-EMOA | 0.0484 ± 0.01 | 0.0764 ± 0.01 | 40 | 1218 | -| RVEA | 0.0510 ± 0.00 | 0.0631 ± 0.00 | 68 | 73 | -| NSGA-III | 0.0197 ± 0.00 | 0.0735 ± 0.01 | 92 | 137 | -| MOEA/D | 0.0037 ± 0.00 | 0.0886 ± 0.00 | 78 | 24 | +| **IBEA** | **0.0014 ± 0.0002** | 0.0607 ± 0.0047 | 87 | 148 | +| MOEA/D | 0.0037 ± 0.0003 | 0.0886 ± 0.0024 | 78 | 23 | +| HypE | 0.0113 ± 0.0033 | **0.0269 ± 0.0172** | 80 | 41 | +| NSGA-III | 0.0197 ± 0.0015 | 0.0735 ± 0.0052 | 92 | 91 | +| eps-MOEA | 0.0325 ± 0.0104 | 0.0572 ± 0.0170 | 136 | 88 | +| NSGA-II | 0.0332 ± 0.0068 | 0.0577 ± 0.0109 | 92 | 60 | +| SPEA2 | 0.0368 ± 0.0021 | 0.0288 ± 0.0038 | 92 | 530 | +| PESA-II | 0.0395 ± 0.0033 | 0.0616 ± 0.0051 | 100 | 372 | +| SMS-EMOA | 0.0484 ± 0.0134 | 0.0764 ± 0.0081 | 40 | 483 | +| RVEA | 0.0510 ± 0.0044 | 0.0631 ± 0.0024 | 68 | 66 | +| MOPSO | 0.0566 ± 0.0048 | 0.0687 ± 0.0084 | 100 | 66 | +| RandomSearch | 0.3949 ± 0.0152 | 0.0797 ± 0.0083 | 239 | 530 | -**IBEA wins decisively** (15× closer to the true front than NSGA-III). +**IBEA wins decisively** (14× closer to the true front than NSGA-III). +SMS-EMOA's wall-clock fell ~2.7× from the v0.4.0 snapshot — the +`hypervolume_nd` rework. ## DTLZ1 (3-obj, dim=7, 30000 evals × 10 seeds) -Linear simplex Pareto front (`Σf = 0.5`). +Deb-Thiele-Laumanns-Zitzler 3-objective; the Pareto front is the linear +simplex `Σf = 0.5` in the positive octant. Hard because a deceptive +multimodal `g` term riddles the approach with a huge number of local +fronts — only fully-converged runs land on the simplex. -| algorithm | mean dist ↓ | spacing ↓ | front | ms | +| algorithm | mean dist ↓ | spacing ↓ | front | ms | |---|---|---|---|---| -| NSGA-III | 5.9130 ± 2.82 | 0.4375 ± 0.22 | 92 | 133 | -| MOEA/D | 2.8022 ± 1.78 | 0.2279 ± 0.22 | 78 | 21 | -| AGE-MOEA | 4.5395 ± 2.21 | 0.3930 ± 0.29 | 90 | 247 | -| **GrEA** | **1.7725 ± 0.99** | **0.0719 ± 0.04** | 72 | 104 | +| **GrEA** | **1.7725 ± 0.9897** | **0.0719 ± 0.0438** | 72 | 62 | +| MOEA/D | 2.8022 ± 1.7807 | 0.2279 ± 0.2247 | 78 | 22 | +| AGE-MOEA | 4.5395 ± 2.2114 | 0.3930 ± 0.2864 | 90 | 193 | +| NSGA-III | 5.9130 ± 2.8212 | 0.4375 ± 0.2212 | 92 | 81 | **GrEA shines on linear fronts** — the grid-based niching matches the geometry better than reference points. ## Rastrigin (dim=5, 50000 evals/run × 10 seeds) -Multimodal trap. Global minimum f = 0 at the origin. +Highly multimodal trap: `f = 10n + Σ(xᵢ² − 10·cos(2π·xᵢ))`. Hard because a +near-quadratic global bowl is overlaid with ~10⁵ regularly spaced local +minima — any greedy step lands in the nearest dimple. Global optimum +`f = 0` at the origin. -| algorithm | best f | ms | +| algorithm | best f | ms | |---|---|---| -| RandomSearch | 1.1064e1 ± 2.54 | 14 | -| HillClimber | 1.5966e1 ± 6.25 | 6 | -| **(1+1)-ES** | **0.0000e0 ± 0.00** | 4 | -| SimulatedAnneal | 3.8540e0 ± 1.48 | 7 | -| PAES | 1.5966e1 ± 6.25 | 10 | -| GA | 7.0913e-8 ± 5.50e-8 | 16 | -| PSO | 7.9598e-1 ± 8.67e-1 | 5 | -| NSGA-II | 4.9270e-5 ± 5.04e-5 | 83 | -| **DE** | **0.0000e0 ± 0.00** | 6 | -| CMA-ES | 2.3453e0 ± 1.49 | 11 | -| **IPOP-CMA-ES** | 1.3423e-1 ± 2.71e-1 | 66 | +| **(1+1)-ES** | **0.0000e0 ± 0.00e0** | 4 | +| **DE** | **0.0000e0 ± 0.00e0** | 6 | +| GA | 7.0913e-8 ± 5.50e-8 | 15 | +| NSGA-II | 4.9270e-5 ± 5.04e-5 | 60 | +| IPOP-CMA-ES | 1.3423e-1 ± 2.71e-1 | 61 | +| PSO | 7.9598e-1 ± 8.67e-1 | 5 | +| CMA-ES | 2.3453e0 ± 1.49e0 | 10 | +| SimulatedAnneal | 3.8540e0 ± 1.48e0 | 7 | +| RandomSearch | 1.1064e1 ± 2.54e0 | 14 | +| HillClimber | 1.5966e1 ± 6.25e0 | 6 | +| PAES | 1.5966e1 ± 6.25e0 | 10 | -(1+1)-ES and DE tie for f = 0. **IPOP-CMA-ES drops vanilla CMA-ES from +(1+1)-ES and DE tie for `f = 0`. **IPOP-CMA-ES drops vanilla CMA-ES from 2.35 → 0.13** — the restart logic does what it should. ## Rosenbrock (dim=5, 30000 evals × 10 seeds) -Smooth non-convex valley. - -| algorithm | best f | ms | -|---|---|---| -| DE | 3.3345e-1 ± 3.01e-1 | 2 | -| PSO | 8.2124e-1 ± 1.58e0 | 2 | -| **CMA-ES** | **3.6207e-29 ± 2.35e-29** | 5 | -| TLBO | 1.8458e-3 ± 1.91e-3 | 1 | -| (1+1)-ES | 2.2115e0 ± 2.70e0 | 1 | -| **Nelder-Mead** | **0.0000e0 ± 0.00** | 1 | -| BO (60 evals) | 3.1725e3 ± 2.92e3 | 40 | - -Nelder-Mead **= 0 exactly**, CMA-ES at machine epsilon. BO at only 60 -evaluations is honestly bad on 5-D Rosenbrock (no kernel -hyperparameter tuning) — included as a reminder that BO needs more -evaluations than a smooth problem actually requires for these other -methods. - -## Ackley (dim=5, 30000 evals × 10 seeds) - -Smoother multimodal landscape than Rastrigin. +Rosenbrock's banana valley: `f = Σ(100·(xᵢ₊₁ − xᵢ²)² + (1 − xᵢ)²)`. Hard +because the minimum sits in a long, bent, near-flat valley — easy to +enter, very slow to crawl along to the tip. Global optimum `f = 0` at the +all-ones point. | algorithm | best f | ms | |---|---|---| -| DE | 4.4409e-16 ± 0.00 | 4 | +| **Nelder-Mead** | **0.0000e0 ± 0.00e0** | 1 | +| CMA-ES | 3.6207e-29 ± 2.35e-29 | 5 | +| TLBO | 1.8458e-3 ± 1.91e-3 | 1 | +| DE | 3.3345e-1 ± 3.01e-1 | 2 | +| PSO | 8.2124e-1 ± 1.58e0 | 2 | +| (1+1)-ES | 2.2115e0 ± 2.70e0 | 1 | +| BO (60 evals) | 3.1725e3 ± 2.92e3 | 39 | + +Nelder-Mead **= 0 exactly**, CMA-ES at machine epsilon. BO at only 60 +evaluations is honestly bad on 5-D Rosenbrock (no kernel hyperparameter +tuning) — included as a reminder that BO needs more evaluations than a +smooth problem actually requires for these other methods. + +## Ackley (dim=5, 30000 evals × 10 seeds) + +Ackley's function: a near-flat outer plateau with shallow ripples +surrounding a single deep, narrow global basin. Hard because the gradient +is almost zero far from the optimum, giving local search little to +follow. Global optimum `f = 0` at the origin. + +| algorithm | best f | ms | +|---|---|---| +| **DE** | **4.4409e-16 ± 0.00e0** | 3 | | PSO | 1.5099e-15 ± 1.63e-15 | 3 | -| CMA-ES | 1.5099e-15 ± 1.63e-15 | 6 | +| CMA-ES | 1.5099e-15 ± 1.63e-15 | 5 | | TLBO | 2.2204e-15 ± 1.78e-15 | 2 | -| BO (60 evals) | 1.9622e1 ± 1.23 | 40 | +| BO (60 evals) | 1.9622e1 ± 1.23e0 | 38 | All conventional methods reach machine precision. BO at 60 evals struggles — same caveat as Rosenbrock. + +--- + +## TSP ring-15 (8000 evals/run × 10 seeds) + +15 equally-spaced cities on the unit circle; minimize the closed tour +length. The space is `(15−1)!/2` distinct tours, but cities in convex +position have no 2-opt local optima — so this instance cleanly separates +methods with good neighbourhood moves (inversion = 2-opt) from blind +recombination / sampling. Known optimum (the polygon perimeter): +**6.2374**. + +| algorithm | tour length ↓ | ms | +|---|---|---| +| **HillClimber** | **6.2374 ± 0.0000** | 0 | +| **SimulatedAnneal** | **6.2374 ± 0.0000** | 0 | +| **TabuSearch** | **6.2374 ± 0.0000** | 0 | +| **AntColony** | **6.2374 ± 0.0000** | 8 | +| GA | 7.0133 ± 0.6725 | 2 | +| RandomSearch | 12.1474 ± 0.6797 | 1 | + +Every local-search method (and Ant Colony) hits the exact optimum — as +theory predicts for convex-position TSP under 2-opt. The GA's order +crossover drifts off the optimum, and random sampling is hopeless. + +## JSS FT06 (8000 evals/run × 10 seeds) + +Fisher & Thompson 1963 6-job × 6-machine job-shop; minimize makespan. +Hard because every job has a fixed machine order, so swapping two +operations can ripple delays across the whole schedule. Known optimum: +**55**. + +| algorithm | makespan ↓ | ms | +|---|---|---| +| **SimulatedAnneal** | **55.2000 ± 0.6000** | 1 | +| TabuSearch | 55.9000 ± 1.4457 | 1 | +| GA | 56.0000 ± 1.5492 | 4 | +| RandomSearch | 58.5000 ± 1.2042 | 4 | +| HillClimber | 62.5000 ± 4.3186 | 0 | + +Simulated annealing gets within 0.4% of the known optimum on average; +greedy hill-climbing stalls in operation-order local optima. + +## Knapsack (30 items, bi-objective, 20000 evals/run × 10 seeds) + +Zitzler-Thiele style 0/1 knapsack: two profit vectors, one capacity (half +the total weight). Hard because the two profit objectives conflict and +the capacity constraint carves feasible regions out of the `2³⁰` +bitstrings. No closed-form optimum; scored by hypervolume vs reference +`[0, 0]` (higher is better). + +| algorithm | hypervolume ↑ | front | ms | +|---|---|---|---| +| **NSGA-II** | **1360468.3 ± 11619.6** | 100 | 39 | +| SPEA2 | 1355615.5 ± 9266.8 | 100 | 213 | +| IBEA | 1352595.5 ± 10183.6 | 99 | 101 | +| NSGA-III | 1346446.0 ± 6922.1 | 100 | 39 | +| RandomSearch | 1118233.1 ± 34150.3 | 9 | 17 | + +The three Pareto EAs land within ~1% of each other; random search finds a +front of only ~9 points and trails badly. diff --git a/examples/compare.rs b/examples/compare.rs index 117414c..a26fa29 100644 --- a/examples/compare.rs +++ b/examples/compare.rs @@ -12,6 +12,9 @@ use std::f64::consts::PI; use std::time::Instant; +use rand::Rng as _; + +use heuropt::core::rng::Rng; use heuropt::metrics::{hypervolume::hypervolume_2d, spacing::spacing}; use heuropt::prelude::*; @@ -283,6 +286,47 @@ fn mean_std(values: &[f64]) -> (f64, f64) { (mean, var.sqrt()) } +/// Print an aligned text table: column 0 left-justified, the rest +/// right-justified. Column widths are derived from the actual cell +/// contents (header *and* every row), so the separator and all rows line +/// up no matter how the value magnitudes vary. +/// +/// All cells must be ASCII — width is measured with `str::len`, so a +/// multi-byte character (e.g. `±`, `ε`) would silently break alignment on +/// terminals that render it at a different column width. Callers format +/// `mean +/- std` rather than `mean ± std` for exactly this reason. +fn print_table(header: &[&str], rows: &[Vec]) { + let ncols = header.len(); + let mut widths: Vec = header.iter().map(|h| h.len()).collect(); + for row in rows { + for c in 0..ncols { + widths[c] = widths[c].max(row[c].len()); + } + } + let fmt_row = |cells: &[String]| -> String { + let mut out = String::new(); + for c in 0..ncols { + if c > 0 { + out.push_str(" "); + } + let w = widths[c]; + if c == 0 { + out.push_str(&format!("{:w$}", cells[c])); + } + } + out + }; + let header_owned: Vec = header.iter().map(|s| s.to_string()).collect(); + let header_line = fmt_row(&header_owned); + println!("{header_line}"); + println!("{}", "-".repeat(header_line.len())); + for row in rows { + println!("{}", fmt_row(row)); + } +} + // ----------------------------------------------------------------------------- // ZDT1 algorithm runners // ----------------------------------------------------------------------------- @@ -1659,13 +1703,12 @@ fn mean_distance_to_dtlz1_front(front: &[Candidate>]) -> f64 { fn run_zdt1_comparison() { println!("== ZDT1 (dim={ZDT1_DIM}, {ZDT1_BUDGET} evals/run × {SEEDS} seeds) =="); - println!("metric arrows: hypervolume↑ (higher better), others↓ (lower better)"); + println!("Zitzler-Deb-Thiele 2-objective benchmark: {ZDT1_DIM} real variables, one smooth"); + println!("convex Pareto front f2 = 1 - sqrt(f1). Hard because 29 of 30 variables must"); + println!("collapse to 0 before the front is even reachable, and only then can the"); + println!("population spread along it. Optimum: mean L2 -> 0 (the front is known exactly)."); + println!("sorted best-first by hypervolume (higher better; spacing / mean L2 lower better)"); println!(); - println!( - "{:<14} {:>16} {:>14} {:>14} {:>10} {:>10}", - "algorithm", "hypervolume", "spacing", "mean L2", "front", "ms", - ); - println!("{}", "-".repeat(82)); let zdt1 = Zdt1 { dim: ZDT1_DIM }; let zdt1_objs = zdt1.objectives(); @@ -1677,7 +1720,7 @@ fn run_zdt1_comparison() { ("MOPSO", zdt1_mopso), ("SPEA2", zdt1_spea2), ("PESA-II", zdt1_pesa2), - ("ε-MOEA", zdt1_epsilon_moea), + ("eps-MOEA", zdt1_epsilon_moea), ("IBEA", zdt1_ibea), ("HypE", zdt1_hype), ("SMS-EMOA", zdt1_sms_emoa), @@ -1687,6 +1730,7 @@ fn run_zdt1_comparison() { ("MOEA/D", zdt1_moead), ]; + let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let hv: Vec = runs @@ -1707,28 +1751,43 @@ fn run_zdt1_comparison() { let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); - println!( - "{:<14} {:>16} {:>14} {:>14} {:>10} {:>10}", - name, - format!("{hv_m:.4}±{hv_s:.4}"), - format!("{sp_m:.4}±{sp_s:.4}"), - format!("{l2_m:.4}±{l2_s:.4}"), - format!("{fs_m:.0}"), - format!("{ms_m:.0}"), - ); + rows.push(( + hv_m, + vec![ + name.to_string(), + format!("{hv_m:.4}+/-{hv_s:.4}"), + format!("{sp_m:.4}+/-{sp_s:.4}"), + format!("{l2_m:.4}+/-{l2_s:.4}"), + format!("{fs_m:.0}"), + format!("{ms_m:.0}"), + ], + )); } + // Higher hypervolume is better. + rows.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table( + &[ + "algorithm", + "hypervolume", + "spacing", + "mean L2", + "front", + "ms", + ], + &table, + ); } fn run_dtlz2_comparison() { println!(); println!("== DTLZ2 (3-obj, dim={DTLZ2_DIM}, {DTLZ2_BUDGET} evals/run × {SEEDS} seeds) =="); - println!("Pareto front: unit sphere octant (Σf²=1, all f≥0); 'mean dist' is |‖f‖−1|"); + println!("Deb-Thiele-Laumanns-Zitzler 3-objective; the Pareto front is the unit-sphere"); + println!("octant (Σf² = 1, all f >= 0) -- a curved 2-D surface embedded in 3-D objective"); + println!("space. Hard because crowding/spacing must work in a higher dimension."); + println!("'mean dist' = |‖f‖ - 1|, so 0 means perfectly on the sphere (the known optimum)."); + println!("sorted best-first by mean dist (lower is better)"); println!(); - println!( - "{:<14} {:>16} {:>14} {:>10} {:>10}", - "algorithm", "mean dist↓", "spacing↓", "front", "ms", - ); - println!("{}", "-".repeat(70)); let dtlz2 = dtlz2_problem(); let dtlz2_objs = dtlz2.objectives(); @@ -1740,7 +1799,7 @@ fn run_dtlz2_comparison() { ("NSGA-II", dtlz2_nsga2), ("SPEA2", dtlz2_spea2), ("PESA-II", dtlz2_pesa2), - ("ε-MOEA", dtlz2_epsilon_moea), + ("eps-MOEA", dtlz2_epsilon_moea), ("IBEA", dtlz2_ibea), ("HypE", dtlz2_hype), ("SMS-EMOA", dtlz2_sms_emoa), @@ -1749,6 +1808,7 @@ fn run_dtlz2_comparison() { ("MOEA/D", dtlz2_moead), ]; + let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let dist: Vec = runs @@ -1767,24 +1827,39 @@ fn run_dtlz2_comparison() { let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); - println!( - "{:<14} {:>16} {:>14} {:>10} {:>10}", - name, - format!("{d_m:.4}±{d_s:.4}"), - format!("{sp_m:.4}±{sp_s:.4}"), - format!("{fs_m:.0}"), - format!("{ms_m:.0}"), - ); + rows.push(( + d_m, + vec![ + name.to_string(), + format!("{d_m:.4}+/-{d_s:.4}"), + format!("{sp_m:.4}+/-{sp_s:.4}"), + format!("{fs_m:.0}"), + format!("{ms_m:.0}"), + ], + )); } + // Lower mean distance to the true front is better. + rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table( + &["algorithm", "mean dist", "spacing", "front", "ms"], + &table, + ); } fn run_rastrigin_comparison() { println!(); println!("== Rastrigin (dim={RASTRIGIN_DIM}, {RASTRIGIN_BUDGET} evals/run × {SEEDS} seeds) =="); - println!("global minimum: f = 0 (lower is better)"); + println!( + "Highly multimodal trap: f = 10n + Σ(x_i² - 10·cos(2π·x_i)), {RASTRIGIN_DIM} dimensions." + ); + println!( + "Hard because a near-quadratic global bowl is overlaid with ~10^{RASTRIGIN_DIM} regularly" + ); + println!("spaced local minima -- any greedy step lands in the nearest dimple."); + println!("Global optimum: f = 0 at the origin."); + println!("sorted best-first (lower is better)"); println!(); - println!("{:<14} {:>20} {:>10}", "algorithm", "best f", "ms"); - println!("{}", "-".repeat(48)); type Runner = fn(u64) -> SoRun; let runners: &[(&str, Runner)] = &[ @@ -1801,6 +1876,7 @@ fn run_rastrigin_comparison() { ("IPOP-CMA-ES", rastrigin_ipop_cma_es), ]; + let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let best: Vec = runs.iter().map(|r| r.best_value).collect(); @@ -1809,22 +1885,31 @@ fn run_rastrigin_comparison() { let (b_m, b_s) = mean_std(&best); let (ms_m, _) = mean_std(&ms); - println!( - "{:<14} {:>20} {:>10}", - name, - format!("{b_m:.4e} ± {b_s:.2e}"), - format!("{ms_m:.0}"), - ); + rows.push(( + b_m, + vec![ + name.to_string(), + format!("{b_m:.4e} +/- {b_s:.2e}"), + format!("{ms_m:.0}"), + ], + )); } + // Lower best f is better. + rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table(&["algorithm", "best f", "ms"], &table); } fn run_rosenbrock_comparison() { println!(); println!("== Rosenbrock (dim={ROSENBROCK_DIM}, {ROSENBROCK_BUDGET} evals × {SEEDS} seeds) =="); - println!("smooth non-convex valley; global minimum f = 0 at all-ones"); + println!( + "Rosenbrock's banana valley: f = Σ(100·(x_{{i+1}} - x_i²)² + (1 - x_i)²), {ROSENBROCK_DIM} dims." + ); + println!("Hard because the minimum sits in a long, bent, near-flat valley -- easy to enter,"); + println!("very slow to crawl along to the tip. Global optimum: f = 0 at the all-ones point."); + println!("sorted best-first (lower is better)"); println!(); - println!("{:<14} {:>20} {:>10}", "algorithm", "best f", "ms"); - println!("{}", "-".repeat(48)); type Runner = fn(u64) -> SoRun; let runners: &[(&str, Runner)] = &[ ("DE", rosenbrock_de), @@ -1835,28 +1920,38 @@ fn run_rosenbrock_comparison() { ("Nelder-Mead", rosenbrock_nelder_mead), ("BO (60 evals)", rosenbrock_bo), ]; + let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let best: Vec = runs.iter().map(|r| r.best_value).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (b_m, b_s) = mean_std(&best); let (ms_m, _) = mean_std(&ms); - println!( - "{:<14} {:>20} {:>10}", - name, - format!("{b_m:.4e} ± {b_s:.2e}"), - format!("{ms_m:.0}"), - ); + rows.push(( + b_m, + vec![ + name.to_string(), + format!("{b_m:.4e} +/- {b_s:.2e}"), + format!("{ms_m:.0}"), + ], + )); } + rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table(&["algorithm", "best f", "ms"], &table); } fn run_ackley_comparison() { println!(); println!("== Ackley (dim={ACKLEY_DIM}, {ACKLEY_BUDGET} evals × {SEEDS} seeds) =="); - println!("smoother multimodal landscape than Rastrigin; global minimum f = 0 at origin"); + println!("Ackley's function: a near-flat outer plateau with shallow ripples surrounding a"); + println!( + "single deep, narrow global basin, {ACKLEY_DIM} dimensions. Hard because the gradient is" + ); + println!("almost zero far from the optimum, giving local search little to follow."); + println!("Global optimum: f = 0 at the origin."); + println!("sorted best-first (lower is better)"); println!(); - println!("{:<14} {:>20} {:>10}", "algorithm", "best f", "ms"); - println!("{}", "-".repeat(48)); type Runner = fn(u64) -> SoRun; let runners: &[(&str, Runner)] = &[ ("DE", ackley_de), @@ -1865,31 +1960,36 @@ fn run_ackley_comparison() { ("TLBO", ackley_tlbo), ("BO (60 evals)", ackley_bo), ]; + let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let best: Vec = runs.iter().map(|r| r.best_value).collect(); let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); let (b_m, b_s) = mean_std(&best); let (ms_m, _) = mean_std(&ms); - println!( - "{:<14} {:>20} {:>10}", - name, - format!("{b_m:.4e} ± {b_s:.2e}"), - format!("{ms_m:.0}"), - ); + rows.push(( + b_m, + vec![ + name.to_string(), + format!("{b_m:.4e} +/- {b_s:.2e}"), + format!("{ms_m:.0}"), + ], + )); } + rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table(&["algorithm", "best f", "ms"], &table); } fn run_zdt3_comparison() { println!(); println!("== ZDT3 (dim={ZDT3_DIM}, {ZDT3_BUDGET} evals × {SEEDS} seeds) =="); - println!("disconnected Pareto front (not contiguous); spread across gaps matters"); + println!("Zitzler-Deb-Thiele 2-objective with a DISCONNECTED front: five separate arcs"); + println!("rather than one curve. Hard because an algorithm has to discover and populate"); + println!("every arc while not stranding solutions in the dominated gaps between them."); + println!("Optimum: cover all five arcs; scored by hypervolume vs [11, 11]."); + println!("sorted best-first by hypervolume (higher better; spacing lower better)"); println!(); - println!( - "{:<14} {:>16} {:>14} {:>10} {:>10}", - "algorithm", "hypervolume↑", "spacing↓", "front", "ms", - ); - println!("{}", "-".repeat(70)); let problem = zdt3_problem(); let objs = problem.objectives(); type Runner = fn(u64) -> MoRun; @@ -1899,6 +1999,7 @@ fn run_zdt3_comparison() { ("IBEA", zdt3_ibea), ("AGE-MOEA", zdt3_age_moea), ]; + let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let hv: Vec = runs @@ -1912,27 +2013,34 @@ fn run_zdt3_comparison() { let (sp_m, sp_s) = mean_std(&sp); let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); - println!( - "{:<14} {:>16} {:>14} {:>10} {:>10}", - name, - format!("{hv_m:.4}±{hv_s:.4}"), - format!("{sp_m:.4}±{sp_s:.4}"), - format!("{fs_m:.0}"), - format!("{ms_m:.0}"), - ); + rows.push(( + hv_m, + vec![ + name.to_string(), + format!("{hv_m:.4}+/-{hv_s:.4}"), + format!("{sp_m:.4}+/-{sp_s:.4}"), + format!("{fs_m:.0}"), + format!("{ms_m:.0}"), + ], + )); } + rows.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table( + &["algorithm", "hypervolume", "spacing", "front", "ms"], + &table, + ); } fn run_dtlz1_comparison() { println!(); println!("== DTLZ1 (3-obj, dim={DTLZ1_DIM}, {DTLZ1_BUDGET} evals × {SEEDS} seeds) =="); - println!("Pareto front: linear simplex Σf=0.5 in the positive octant"); + println!("Deb-Thiele-Laumanns-Zitzler 3-objective; the Pareto front is the linear simplex"); + println!("Σf = 0.5 in the positive octant. Hard because a deceptive multimodal 'g' term"); + println!("riddles the approach with a huge number of local fronts -- only fully-converged"); + println!("runs land on the simplex. Optimum: mean dist -> 0."); + println!("sorted best-first by mean dist (lower is better)"); println!(); - println!( - "{:<14} {:>16} {:>14} {:>10} {:>10}", - "algorithm", "mean dist↓", "spacing↓", "front", "ms", - ); - println!("{}", "-".repeat(70)); let problem = dtlz1_problem(); let objs = problem.objectives(); type Runner = fn(u64) -> MoRun; @@ -1942,6 +2050,7 @@ fn run_dtlz1_comparison() { ("AGE-MOEA", dtlz1_age_moea), ("GrEA", dtlz1_grea), ]; + let mut rows: Vec<(f64, Vec)> = Vec::new(); for (name, runner) in runners { let runs: Vec = (0..SEEDS).map(runner).collect(); let dist: Vec = runs @@ -1955,15 +2064,763 @@ fn run_dtlz1_comparison() { let (sp_m, sp_s) = mean_std(&sp); let (fs_m, _) = mean_std(&fs); let (ms_m, _) = mean_std(&ms); - println!( - "{:<14} {:>16} {:>14} {:>10} {:>10}", - name, - format!("{d_m:.4}±{d_s:.4}"), - format!("{sp_m:.4}±{sp_s:.4}"), - format!("{fs_m:.0}"), - format!("{ms_m:.0}"), - ); + rows.push(( + d_m, + vec![ + name.to_string(), + format!("{d_m:.4}+/-{d_s:.4}"), + format!("{sp_m:.4}+/-{sp_s:.4}"), + format!("{fs_m:.0}"), + format!("{ms_m:.0}"), + ], + )); } + rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table( + &["algorithm", "mean dist", "spacing", "front", "ms"], + &table, + ); +} + +// ============================================================================= +// Combinatorial / sequencing problems +// +// These use a different decision encoding — `Vec` permutations and +// `Vec` bitstrings — than the continuous problems above, so the +// applicable-algorithm roster differs: the real-vector methods (CMA-ES, DE, +// PSO, …) cannot run here, while permutation-native methods (Ant Colony, +// Tabu Search) and the EDA-style ones can. +// ============================================================================= + +const TSP_CITIES: usize = 15; +const TSP_BUDGET: usize = 8_000; +const JSS_BUDGET: usize = 8_000; +const KNAPSACK_BUDGET: usize = 20_000; + +// ---- Single-objective TSP: a regular polygon on the unit circle ------------- + +/// `TSP_CITIES` points equally spaced on the unit circle. The optimal tour +/// just visits them in angular order; its length is the regular-polygon +/// perimeter `2·n·sin(π/n)`, which gives an exact known-best baseline. +struct RingTsp { + distances: Vec>, +} +impl RingTsp { + fn new() -> Self { + let pts: Vec<(f64, f64)> = (0..TSP_CITIES) + .map(|i| { + let a = 2.0 * PI * (i as f64) / (TSP_CITIES as f64); + (a.cos(), a.sin()) + }) + .collect(); + let distances = (0..TSP_CITIES) + .map(|i| { + (0..TSP_CITIES) + .map(|j| { + let (xi, yi) = pts[i]; + let (xj, yj) = pts[j]; + ((xi - xj).powi(2) + (yi - yj).powi(2)).sqrt() + }) + .collect() + }) + .collect(); + Self { distances } + } + /// Known optimal tour length — the regular-polygon perimeter. + fn optimal_length() -> f64 { + 2.0 * TSP_CITIES as f64 * (PI / TSP_CITIES as f64).sin() + } +} +impl Problem for RingTsp { + type Decision = Vec; + fn objectives(&self) -> ObjectiveSpace { + ObjectiveSpace::new(vec![Objective::minimize("tour_length")]) + } + fn evaluate(&self, tour: &Vec) -> Evaluation { + let n = tour.len(); + let mut total = 0.0; + for i in 0..n { + total += self.distances[tour[i]][tour[(i + 1) % n]]; + } + Evaluation::new(vec![total]) + } +} + +fn tsp_random(seed: u64) -> SoRun { + let problem = RingTsp::new(); + let mut opt = RandomSearch::new( + RandomSearchConfig { + iterations: TSP_BUDGET, + batch_size: 1, + seed, + }, + ShuffledPermutation { n: TSP_CITIES }, + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn tsp_hill_climber(seed: u64) -> SoRun { + let problem = RingTsp::new(); + let mut opt = HillClimber::new( + HillClimberConfig { + iterations: TSP_BUDGET, + seed, + }, + ShuffledPermutation { n: TSP_CITIES }, + InversionMutation, + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn tsp_simulated_annealing(seed: u64) -> SoRun { + let problem = RingTsp::new(); + let mut opt = SimulatedAnnealing::new( + SimulatedAnnealingConfig { + iterations: TSP_BUDGET, + initial_temperature: 1.0, + final_temperature: 1e-3, + seed, + }, + ShuffledPermutation { n: TSP_CITIES }, + InversionMutation, + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn tsp_tabu_search(seed: u64) -> SoRun { + let problem = RingTsp::new(); + // Each step considers 16 inversion (2-opt-style) neighbours. + let neighbors = |tour: &Vec, rng: &mut Rng| { + let mut m = InversionMutation; + (0..16) + .map(|_| m.vary(std::slice::from_ref(tour), rng).pop().unwrap()) + .collect() + }; + let mut opt = TabuSearch::new( + TabuSearchConfig { + iterations: TSP_BUDGET / 16, + tabu_tenure: 20, + seed, + }, + ShuffledPermutation { n: TSP_CITIES }, + neighbors, + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn tsp_genetic_algorithm(seed: u64) -> SoRun { + let problem = RingTsp::new(); + let mut opt = GeneticAlgorithm::new( + GeneticAlgorithmConfig { + population_size: 80, + generations: TSP_BUDGET / 80, + tournament_size: 3, + elitism: 2, + seed, + }, + ShuffledPermutation { n: TSP_CITIES }, + CompositeVariation { + crossover: OrderCrossover, + mutation: InversionMutation, + }, + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn tsp_ant_colony(seed: u64) -> SoRun { + let problem = RingTsp::new(); + let distances = problem.distances.clone(); + let mut opt = AntColonyTsp::new( + AntColonyTspConfig { + ants: 20, + generations: TSP_BUDGET / 20, + alpha: 1.0, + beta: 3.0, + evaporation: 0.5, + deposit: 1.0, + initial_pheromone: 1.0, + seed, + }, + distances, + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +// ---- Single-objective job-shop scheduling: Fisher & Thompson FT06 ----------- + +const JSS_JOBS: usize = 6; +const JSS_MACHINES: usize = 6; + +/// FT06 routing — machine id of the k-th operation of job j. +const FT06_MACHINE: [[usize; JSS_MACHINES]; JSS_JOBS] = [ + [2, 0, 1, 3, 5, 4], + [1, 2, 4, 5, 0, 3], + [2, 3, 5, 0, 1, 4], + [1, 0, 2, 3, 4, 5], + [2, 1, 4, 5, 0, 3], + [1, 3, 5, 0, 4, 2], +]; +/// FT06 processing times — duration of the k-th operation of job j. +const FT06_TIME: [[f64; JSS_MACHINES]; JSS_JOBS] = [ + [1.0, 3.0, 6.0, 7.0, 3.0, 6.0], + [8.0, 5.0, 10.0, 10.0, 10.0, 4.0], + [5.0, 4.0, 8.0, 9.0, 1.0, 7.0], + [5.0, 5.0, 5.0, 3.0, 8.0, 9.0], + [9.0, 3.0, 5.0, 4.0, 3.0, 1.0], + [3.0, 3.0, 9.0, 10.0, 4.0, 1.0], +]; +/// FT06's optimal makespan is a long-settled benchmark value. +const FT06_OPTIMAL_MAKESPAN: f64 = 55.0; + +/// FT06 job-shop, makespan objective. Operation-string encoding: a length-36 +/// multiset where job id `j` appears `JSS_MACHINES` times; the k-th +/// occurrence of `j` is its k-th operation. +struct Ft06Makespan; +impl Problem for Ft06Makespan { + type Decision = Vec; + fn objectives(&self) -> ObjectiveSpace { + ObjectiveSpace::new(vec![Objective::minimize("makespan")]) + } + fn evaluate(&self, schedule: &Vec) -> Evaluation { + let mut job_next = [0_usize; JSS_JOBS]; + let mut job_clock = [0.0_f64; JSS_JOBS]; + let mut machine_clock = [0.0_f64; JSS_MACHINES]; + for &job in schedule { + let k = job_next[job]; + let m = FT06_MACHINE[job][k]; + let t = FT06_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]) + } +} + +/// Precedence-Order Crossover — multiset-preserving recombination for the +/// operation-string encoding (the strict-permutation crossovers would +/// corrupt the multiset). +#[derive(Debug, Clone, Copy, Default)] +struct PrecedenceOrderCrossover; +impl Variation> for PrecedenceOrderCrossover { + fn vary(&mut self, parents: &[Vec], rng: &mut Rng) -> Vec> { + assert!(parents.len() >= 2, "POX requires 2 parents"); + let (p1, p2) = (&parents[0], &parents[1]); + let mut in_j1 = [false; JSS_JOBS]; + loop { + for slot in &mut in_j1 { + *slot = rng.random_bool(0.5); + } + let c = in_j1.iter().filter(|&&b| b).count(); + if c > 0 && c < JSS_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 { + 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 fill_idx = 0; + for &v in filler { + if !in_donor_set[v] { + while fill_idx < n && child[fill_idx] != usize::MAX { + fill_idx += 1; + } + child[fill_idx] = v; + fill_idx += 1; + } + } + child +} + +fn jss_initializer() -> ShuffledMultisetPermutation { + ShuffledMultisetPermutation::new(vec![JSS_MACHINES; JSS_JOBS]) +} + +fn jss_random(seed: u64) -> SoRun { + let mut opt = RandomSearch::new( + RandomSearchConfig { + iterations: JSS_BUDGET, + batch_size: 1, + seed, + }, + jss_initializer(), + ); + let t0 = Instant::now(); + let result = opt.run(&Ft06Makespan); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn jss_hill_climber(seed: u64) -> SoRun { + let mut opt = HillClimber::new( + HillClimberConfig { + iterations: JSS_BUDGET, + seed, + }, + jss_initializer(), + InsertionMutation, + ); + let t0 = Instant::now(); + let result = opt.run(&Ft06Makespan); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn jss_simulated_annealing(seed: u64) -> SoRun { + let mut opt = SimulatedAnnealing::new( + SimulatedAnnealingConfig { + iterations: JSS_BUDGET, + initial_temperature: 5.0, + final_temperature: 1e-2, + seed, + }, + jss_initializer(), + InsertionMutation, + ); + let t0 = Instant::now(); + let result = opt.run(&Ft06Makespan); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn jss_tabu_search(seed: u64) -> SoRun { + let neighbors = |schedule: &Vec, rng: &mut Rng| { + let mut m = InsertionMutation; + (0..16) + .map(|_| m.vary(std::slice::from_ref(schedule), rng).pop().unwrap()) + .collect() + }; + let mut opt = TabuSearch::new( + TabuSearchConfig { + iterations: JSS_BUDGET / 16, + tabu_tenure: 20, + seed, + }, + jss_initializer(), + neighbors, + ); + let t0 = Instant::now(); + let result = opt.run(&Ft06Makespan); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +fn jss_genetic_algorithm(seed: u64) -> SoRun { + let mut opt = GeneticAlgorithm::new( + GeneticAlgorithmConfig { + population_size: 80, + generations: JSS_BUDGET / 80, + tournament_size: 3, + elitism: 2, + seed, + }, + jss_initializer(), + CompositeVariation { + crossover: PrecedenceOrderCrossover, + mutation: InsertionMutation, + }, + ); + let t0 = Instant::now(); + let result = opt.run(&Ft06Makespan); + SoRun { + best_value: result.best.unwrap().evaluation.objectives[0], + wall_ms: t0.elapsed().as_millis(), + } +} + +// ---- Bi-objective 0/1 knapsack: Zitzler & Thiele style ---------------------- + +const KNAPSACK_N: usize = 30; +const KP_PROFIT_A: [f64; KNAPSACK_N] = [ + 61.0, 17.0, 92.0, 49.0, 73.0, 28.0, 84.0, 36.0, 55.0, 78.0, 23.0, 91.0, 12.0, 67.0, 45.0, 58.0, + 33.0, 71.0, 14.0, 26.0, 87.0, 42.0, 19.0, 65.0, 30.0, 51.0, 79.0, 22.0, 47.0, 88.0, +]; +const KP_PROFIT_B: [f64; KNAPSACK_N] = [ + 24.0, 81.0, 16.0, 67.0, 29.0, 73.0, 41.0, 60.0, 52.0, 19.0, 77.0, 34.0, 95.0, 22.0, 71.0, 88.0, + 56.0, 27.0, 64.0, 90.0, 18.0, 43.0, 79.0, 31.0, 85.0, 25.0, 38.0, 92.0, 70.0, 13.0, +]; +const KP_WEIGHT: [f64; KNAPSACK_N] = [ + 35.0, 58.0, 22.0, 71.0, 14.0, 86.0, 31.0, 53.0, 78.0, 19.0, 44.0, 16.0, 67.0, 88.0, 25.0, 51.0, + 33.0, 74.0, 12.0, 47.0, 63.0, 28.0, 91.0, 36.0, 55.0, 17.0, 82.0, 41.0, 24.0, 68.0, +]; +/// Hypervolume reference point for the knapsack front. Both objectives are +/// maximized, so in the minimization-oriented frame any positive profit is +/// below 0 — `[0, 0]` is dominated by every feasible solution. +const KNAPSACK_REFERENCE: [f64; 2] = [0.0, 0.0]; + +/// Bi-objective 0/1 knapsack: two profit vectors, one shared capacity. +/// Weight overruns are penalized in both objectives so the Pareto front is +/// composed of feasible solutions. +struct BiKnapsack { + capacity: f64, +} +impl BiKnapsack { + fn new() -> Self { + Self { + capacity: 0.5 * KP_WEIGHT.iter().sum::(), + } + } +} +impl Problem for BiKnapsack { + type Decision = Vec; + fn objectives(&self) -> ObjectiveSpace { + ObjectiveSpace::new(vec![ + Objective::maximize("profit_a"), + Objective::maximize("profit_b"), + ]) + } + fn evaluate(&self, take: &Vec) -> Evaluation { + let (mut pa, mut pb, mut w) = (0.0, 0.0, 0.0); + for (i, &t) in take.iter().enumerate() { + if t { + pa += KP_PROFIT_A[i]; + pb += KP_PROFIT_B[i]; + w += KP_WEIGHT[i]; + } + } + let penalty = 1000.0 * (w - self.capacity).max(0.0); + Evaluation::new(vec![pa - penalty, pb - penalty]) + } +} + +/// Random binary initializer — each bit 50/50 independently. +#[derive(Debug, Clone, Copy)] +struct RandomBinary { + n: usize, +} +impl Initializer> for RandomBinary { + fn initialize(&mut self, size: usize, rng: &mut Rng) -> Vec> { + (0..size) + .map(|_| (0..self.n).map(|_| rng.random_bool(0.5)).collect()) + .collect() + } +} + +/// One-point crossover for binary chromosomes. +#[derive(Debug, Clone, Copy, Default)] +struct OnePointCrossoverBool; +impl Variation> for OnePointCrossoverBool { + fn vary(&mut self, parents: &[Vec], rng: &mut Rng) -> Vec> { + assert!( + parents.len() >= 2, + "OnePointCrossoverBool requires 2 parents" + ); + let (p1, p2) = (&parents[0], &parents[1]); + let n = p1.len(); + if n < 2 { + return vec![p1.clone(), p2.clone()]; + } + 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 knap_binary_variation() -> CompositeVariation { + CompositeVariation { + crossover: OnePointCrossoverBool, + mutation: BitFlipMutation { + probability: 1.0 / KNAPSACK_N as f64, + }, + } +} + +/// One run's knapsack metrics. The bi-objective front lives over `Vec` +/// decisions, so `MoRun` (which fixes `Vec`) does not fit — we keep +/// only the aggregate quality numbers. +#[derive(Clone)] +struct KnapRun { + hypervolume: f64, + front_size: usize, + wall_ms: u128, +} + +fn knap_run(front: &[Candidate>], problem: &BiKnapsack, wall_ms: u128) -> KnapRun { + KnapRun { + hypervolume: hypervolume_2d(front, &problem.objectives(), KNAPSACK_REFERENCE), + front_size: front.len(), + wall_ms, + } +} + +fn knapsack_random(seed: u64) -> KnapRun { + let problem = BiKnapsack::new(); + let mut opt = RandomSearch::new( + RandomSearchConfig { + iterations: KNAPSACK_BUDGET, + batch_size: 1, + seed, + }, + RandomBinary { n: KNAPSACK_N }, + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + let ms = t0.elapsed().as_millis(); + knap_run(&result.pareto_front, &problem, ms) +} + +fn knapsack_nsga2(seed: u64) -> KnapRun { + let problem = BiKnapsack::new(); + let mut opt = Nsga2::new( + Nsga2Config { + population_size: 100, + generations: KNAPSACK_BUDGET / 100, + seed, + }, + RandomBinary { n: KNAPSACK_N }, + knap_binary_variation(), + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + let ms = t0.elapsed().as_millis(); + knap_run(&result.pareto_front, &problem, ms) +} + +fn knapsack_spea2(seed: u64) -> KnapRun { + let problem = BiKnapsack::new(); + let mut opt = Spea2::new( + Spea2Config { + population_size: 100, + archive_size: 100, + generations: KNAPSACK_BUDGET / 100, + seed, + }, + RandomBinary { n: KNAPSACK_N }, + knap_binary_variation(), + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + let ms = t0.elapsed().as_millis(); + knap_run(&result.pareto_front, &problem, ms) +} + +fn knapsack_nsga3(seed: u64) -> KnapRun { + let problem = BiKnapsack::new(); + let mut opt = Nsga3::new( + Nsga3Config { + population_size: 100, + generations: KNAPSACK_BUDGET / 100, + reference_divisions: 20, + seed, + }, + RandomBinary { n: KNAPSACK_N }, + knap_binary_variation(), + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + let ms = t0.elapsed().as_millis(); + knap_run(&result.pareto_front, &problem, ms) +} + +fn knapsack_ibea(seed: u64) -> KnapRun { + let problem = BiKnapsack::new(); + let mut opt = Ibea::new( + IbeaConfig { + population_size: 100, + generations: KNAPSACK_BUDGET / 100, + kappa: 0.05, + seed, + }, + RandomBinary { n: KNAPSACK_N }, + knap_binary_variation(), + ); + let t0 = Instant::now(); + let result = opt.run(&problem); + let ms = t0.elapsed().as_millis(); + knap_run(&result.pareto_front, &problem, ms) +} + +// ---- Combinatorial comparison runners -------------------------------------- + +fn run_tsp_comparison() { + println!(); + println!("== TSP ring-{TSP_CITIES} ({TSP_BUDGET} evals/run × {SEEDS} seeds) =="); + println!( + "{TSP_CITIES} equally-spaced cities on the unit circle; minimize the closed tour length." + ); + println!( + "The space is ({TSP_CITIES}-1)!/2 distinct tours, but cities in convex position have no" + ); + println!("2-opt local optima -- so this instance cleanly separates methods with good"); + println!("neighbourhood moves (inversion = 2-opt) from blind recombination / sampling."); + println!( + "Known optimum (the polygon perimeter): {:.4}", + RingTsp::optimal_length() + ); + println!("sorted best-first by tour length (lower is better)"); + println!(); + + type Runner = fn(u64) -> SoRun; + let runners: &[(&str, Runner)] = &[ + ("RandomSearch", tsp_random), + ("HillClimber", tsp_hill_climber), + ("SimulatedAnneal", tsp_simulated_annealing), + ("TabuSearch", tsp_tabu_search), + ("GA", tsp_genetic_algorithm), + ("AntColony", tsp_ant_colony), + ]; + + let mut rows: Vec<(f64, Vec)> = Vec::new(); + for (name, runner) in runners { + let runs: Vec = (0..SEEDS).map(runner).collect(); + let best: Vec = runs.iter().map(|r| r.best_value).collect(); + let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); + let (b_m, b_s) = mean_std(&best); + let (ms_m, _) = mean_std(&ms); + rows.push(( + b_m, + vec![ + name.to_string(), + format!("{b_m:.4}+/-{b_s:.4}"), + format!("{ms_m:.0}"), + ], + )); + } + rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table(&["algorithm", "tour length", "ms"], &table); +} + +fn run_jss_comparison() { + println!(); + println!("== JSS FT06 ({JSS_BUDGET} evals/run × {SEEDS} seeds) =="); + println!("Fisher & Thompson 1963 6-job × 6-machine job-shop; minimize makespan."); + println!("Hard because every job has a fixed machine order, so swapping two"); + println!( + "operations can ripple delays across the whole schedule. Known optimum: {FT06_OPTIMAL_MAKESPAN:.0}" + ); + println!("sorted best-first by makespan (lower is better)"); + println!(); + + type Runner = fn(u64) -> SoRun; + let runners: &[(&str, Runner)] = &[ + ("RandomSearch", jss_random), + ("HillClimber", jss_hill_climber), + ("SimulatedAnneal", jss_simulated_annealing), + ("TabuSearch", jss_tabu_search), + ("GA", jss_genetic_algorithm), + ]; + + let mut rows: Vec<(f64, Vec)> = Vec::new(); + for (name, runner) in runners { + let runs: Vec = (0..SEEDS).map(runner).collect(); + let best: Vec = runs.iter().map(|r| r.best_value).collect(); + let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); + let (b_m, b_s) = mean_std(&best); + let (ms_m, _) = mean_std(&ms); + rows.push(( + b_m, + vec![ + name.to_string(), + format!("{b_m:.4}+/-{b_s:.4}"), + format!("{ms_m:.0}"), + ], + )); + } + rows.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table(&["algorithm", "makespan", "ms"], &table); +} + +fn run_knapsack_comparison() { + println!(); + println!( + "== Knapsack ({KNAPSACK_N} items, bi-objective, {KNAPSACK_BUDGET} evals/run × {SEEDS} seeds) ==" + ); + println!("Zitzler-Thiele style 0/1 knapsack: two profit vectors, one capacity"); + println!("(half the total weight). Hard because the two profit objectives"); + println!("conflict and the capacity constraint carves feasible regions out of"); + println!("the 2^{KNAPSACK_N} bitstrings. No closed-form optimum; scored by hypervolume"); + println!("vs reference {KNAPSACK_REFERENCE:?} (higher is better)."); + println!("sorted best-first by hypervolume"); + println!(); + + type Runner = fn(u64) -> KnapRun; + let runners: &[(&str, Runner)] = &[ + ("RandomSearch", knapsack_random), + ("NSGA-II", knapsack_nsga2), + ("SPEA2", knapsack_spea2), + ("NSGA-III", knapsack_nsga3), + ("IBEA", knapsack_ibea), + ]; + + let mut rows: Vec<(f64, Vec)> = Vec::new(); + for (name, runner) in runners { + let runs: Vec = (0..SEEDS).map(runner).collect(); + let hv: Vec = runs.iter().map(|r| r.hypervolume).collect(); + let fs: Vec = runs.iter().map(|r| r.front_size as f64).collect(); + let ms: Vec = runs.iter().map(|r| r.wall_ms as f64).collect(); + let (hv_m, hv_s) = mean_std(&hv); + let (fs_m, _) = mean_std(&fs); + let (ms_m, _) = mean_std(&ms); + rows.push(( + hv_m, + vec![ + name.to_string(), + format!("{hv_m:.1}+/-{hv_s:.1}"), + format!("{fs_m:.0}"), + format!("{ms_m:.0}"), + ], + )); + } + rows.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal)); + let table: Vec> = rows.into_iter().map(|(_, r)| r).collect(); + print_table(&["algorithm", "hypervolume", "front", "ms"], &table); } fn main() { @@ -1974,4 +2831,7 @@ fn main() { run_rastrigin_comparison(); run_rosenbrock_comparison(); run_ackley_comparison(); + run_tsp_comparison(); + run_jss_comparison(); + run_knapsack_comparison(); }