diff --git a/.grok/skills/unsga3-oracle/SKILL.md b/.grok/skills/unsga3-oracle/SKILL.md index f5336a0..cade722 100644 --- a/.grok/skills/unsga3-oracle/SKILL.md +++ b/.grok/skills/unsga3-oracle/SKILL.md @@ -19,9 +19,11 @@ Repo root: Unsga3. Confirm `Unsga3.slnx` / `tools/OracleCompare` exist before ru | zdt2 | 12 | 52 | **250** | `--pymoo-mode` (`PymooCompatible`) | | dtlz2 | 12 | 92 | 150 | `--pymoo-mode` (`PymooCompatible`) | +DTLZ2: C# `Dtlz2Problem(k: 10)` ⇒ **n_var=12**. `run_pymoo_oracle.py` passes `n_var=12`. pymoo’s own default is n_var=10 (k=8). Do not treat the published 15-seed pymoo column as that matched run, and do not rewrite `docs/WILCOXON-RESULTS.md` until the seeds are re-run. + ZDT2 **gens=100** is an early-stress snapshot (collapse on Bend, C#, and pymoo), not the quality bar. Quality protocol matches unsga3-bend A/B (gens=250, PymooCompatible). `RankNicheDistance` is an optional unpublished Wilcoxon ZDT2 mating mode — do not silently switch all ZDT defaults to it. ZDT1 and DTLZ2 unchanged. -IGD = **mean** nearest Euclidean distance (pymoo-compatible). Docs: `docs/EQUIVALENCE.md`, `docs/RESEARCH-STANDARDS.md`. +IGD = **mean** nearest Euclidean distance (pymoo-compatible). C# scores the full non-dominated front; `run_pymoo_oracle.py` scores pymoo `res.F` (niche optimum). Compare them only on a shared front definition and a shared reference set. Docs: `docs/EQUIVALENCE.md`, `docs/RESEARCH-STANDARDS.md`. Requires: .NET 10 SDK; Python 3 + `pip install pymoo` for pymoo side / multi-seed. diff --git a/CHANGELOG.md b/CHANGELOG.md index 23df1ce..924adf2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/). ### Changed - Docs and XML comments describe `initialPopulation` and hybrid loops in generic terms (domain-adapter warm-start / grid-seed). No product-repo names. +- DTLZ2 pymoo oracle passes `n_var=12` (k=10) to match `Dtlz2Problem`. The published 15-seed table used pymoo’s default `n_var=10` and is not rewritten. Seed 1 was remeasured at n_var=12. +- Documented that ZDT IGD compares the C# full non-dominated front with pymoo `res.F`. Added `ReferenceDirectionThinning.OnePerDirection` as a cardinality aid, not a parity claim. +- Documented the collapsed-nadir fallback (nadir = ideal + 1 when the span stays ≤ 1e-6). pymoo 0.6.2 stops at the worst point in the population. Behavior is unchanged and covered by a unit test. +- Documented that default `RankNicheDistance` is not Seada and Deb Algorithm 2. `PymooCompatible` matches the paper's same-niche split; p_c stays 1.0 (paper experiments use 0.9). The default tournament is unchanged. +- Locked the infeasible-point hyperplane rule with a fixture: feasible (1, 1) beside infeasible (0, 0) sets ideal to (0, 0). The rule is unchanged. +- Documented that mating calls `PrepareForSelection`, which re-associates survivors. pymoo keeps the niche ids from survival. A fixture locks the current ids on a five-point pool. +- Documented that `WithDasDennis(1, 1)` throws because N must be at least 2. Single-objective runs pass an explicit population size. N = 1 is not accepted. +- Indicator edges: Euclidean distance rejects a shorter vector instead of ignoring the extra coordinates. IGD+ has a hand-case test. `ParetoFronts.Zdt1(1)` (and the other single-point samplers) throw. ZDT6's 0.280775 floor and ZDT3's 0.1822287280 endpoint are documented and locked. +- Duplicate elimination no longer spins when mutation cannot change the decision vector. The key is `G12` significant digits, not 12 decimal places. After the attempt cap, remaining offspring slots may be duplicates. +- Equivalence docs no longer say the published Wilcoxon table is within 1–2% of pymoo. The DTLZ2 seed-1 guard is 2× the published mismatched scalar 0.00350, so a regression to about 2.9× fails. Das–Dennis `Count` uses a checked 64-bit combination and throws when the value does not fit in `int`. Two-layer reference directions remain absent. ## [0.1.4] — 2026-09-19 diff --git a/README.md b/README.md index 0bbb696..e7b2f3e 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,7 @@ **U-NSGA-III** (Unified NSGA-III) for .NET — single-, multi-, and many-objective evolutionary optimization with Das–Dennis reference directions, SBX crossover, polynomial mutation, and **niching-based tournament selection** ([Seada & Deb, 2016](https://ieeexplore.ieee.org/document/7271063)). > **v0.1.4** — production-usable core with pymoo-aligned normalization. ZDT2 quality protocol is gens=250 (docs/defaults; no algorithm/API change). -> **15-seed IGD vs pymoo `UNSGA3`:** ZDT1 **median 0.053 vs 0.070** (we win; MWU *p*≈0.05); DTLZ2 **median 0.0045 vs 0.0028** (~1.6×, same order; pymoo still ahead). +> **15-seed IGD vs pymoo `UNSGA3`:** ZDT1 **median 0.053 vs 0.070** (MWU *p*≈0.05) compares the full C# non-dominated front with pymoo `res.F`. DTLZ2 **median 0.0045 vs 0.0028** (~1.6×) compares C# n_var=12 with pymoo's default n_var=10. Neither pair is a same-set, same-problem ranking. Notes: [`docs/ORACLE-RESULTS.md`](docs/ORACLE-RESULTS.md). > Details: [`docs/WILCOXON-RESULTS.md`](docs/WILCOXON-RESULTS.md) · single-seed notes: [`docs/ORACLE-RESULTS.md`](docs/ORACLE-RESULTS.md) ```text diff --git a/docs/EQUIVALENCE.md b/docs/EQUIVALENCE.md index dcdca66..f4c9b52 100644 --- a/docs/EQUIVALENCE.md +++ b/docs/EQUIVALENCE.md @@ -1,6 +1,6 @@ # Equivalence vs pymoo / MATLAB -Goal: prove this port is a faithful U-NSGA-III (Seada & Deb 2016), not a look-alike. +Goal: state where this port follows Seada & Deb 2016 and pymoo, and where it does not. Survival, Das–Dennis directions, and the SBX/PM shapes are the pymoo-shaped core. The default tournament is not Algorithm 2. See also **[RESEARCH-STANDARDS.md](RESEARCH-STANDARDS.md)** for the literature + pymoo protocol, **[ORACLE-RESULTS.md](ORACLE-RESULTS.md)** for single-seed numbers, and @@ -17,10 +17,14 @@ See also **[RESEARCH-STANDARDS.md](RESEARCH-STANDARDS.md)** for the literature + 1. **Fixed operators:** SBX η=30, PM η=20, p_c=1.0, p_m=1/n, p_var(SBX)=0.5 2. **Same reference set:** Das–Dennis partitions identical to the oracle -3. **Same pop size / generations / seed** (or 15–31 seeds for statistics) -4. **Metrics:** IGD (primary), IGD+, HV (M=2, document ref point), front plots for M≤3 -5. **Tolerance:** median IGD within ~1–2× of pymoo on ZDT/DTLZ is the practical bar. - 15-seed: ZDT1 median **better** than pymoo (ratio 0.76, MWU n.s.); DTLZ2 median ~**1.6×** (pymoo still ahead). +3. **Same decision dimension:** DTLZ2 uses k=10 so `n_var = M + k − 1` (12 when M=3). pymoo’s default `n_var=10` (k=8) is a known mismatch; the oracle passes `n_var=12`. +4. **Same pop size / generations / seed** (or 15–31 seeds for statistics) +5. **Metrics:** IGD (primary), IGD+, HV (M=2, document ref point), front plots for M≤3. + Score the **same front definition** and the **same reference set**. C# reports the full non-dominated front. pymoo `res.F` is the survival niche set (about one point per filled direction). `ReferenceDirectionThinning.OnePerDirection` can match cardinality; it does not reproduce `res.F` and is not a parity claim. +6. **Tolerance:** do not read the published table as median IGD within 1–2% of pymoo. + ZDT1 median ratio 0.764107, Mann–Whitney U = 65, p = 0.0512394 (the pymoo column is `res.F`). + DTLZ2 median ratio 1.58638, U = 222, p = 6.15164×10⁻⁶ (the pymoo column is n_var=10). That comparison rejects equal distributions at α = 0.05. + The CI guard is 2× the published mismatched DTLZ2 scalar 0.00350, which fails a regression to about 2.9×. It is not a same-problem equivalence claim. Published A/B budgets (ZDT1 / DTLZ2 unchanged; ZDT2 matches unsga3-bend protocol honesty): @@ -59,12 +63,20 @@ ZDT2 **gens=100** is an early-stress snapshot (collapse on Bend, C#, and pymoo), | Item | This library | pymoo | |------|--------------|-------| -| Tournament (default) | rank → niche count → dist | — | -| Tournament (`PymooCompatible`) | same niche → rank/dist; else random | `comp_by_rank_and_ref_line_dist` | -| Duplicate elimination | default **on** | `eliminate_duplicates=True` | +| Tournament (default `RankNicheDistance`) | rank → niche count → dist, including across niches | not Algorithm 2 | +| Tournament (`PymooCompatible`) | same niche → rank then dist; else random; distance tie is a coin flip | `comp_by_rank_and_ref_line_dist` (paper keeps the second parent on a distance tie) | +| SBX p_c | **1.0** (pymoo `SBX(prob=1.0)`) | paper section 4 uses **0.9** | +| Mating pool | N independent tournaments with replacement | two shuffled consecutive-pair passes | +| Niche ids at mating | `PrepareForSelection` re-normalizes survivors and re-associates | ids written during survival are kept | +| `WithDasDennis(1, 1)` | throws. One objective has a single direction, and N must be ≥ 2 | pass `populationSize` ≥ 2 for the single-objective degeneration | +| Reference layers | single-layer Das–Dennis only. Two-layer directions are absent | many-objective NSGA-III adds an inside layer for larger M | +| Duplicate elimination | default **on**. Key is `G12` (12 significant digits), not 12 decimal places. Attempts are capped when mutation cannot produce a new key; remaining slots may be duplicates | `eliminate_duplicates=True` | | Survival RNG | optional RNG niche pick | random among equal niches | | IGD | **mean** nearest distance | same (verified pymoo 0.6.2) | +| Scored set | full non-dominated front | `res.F` niche optimum | | Hyperplane norm | persistent ideal, ND extremes, correct ASF | `HyperplaneNormalization` | +| Collapsed nadir | if the span is still ≤ 1e-6, nadir = ideal + 1 | stop at worst-of-population | +| Infeasible points in the hyperplane | ideal and worst from the whole pool, including infeasible points. ASF extremes use the ND index set when supplied | pymoo niching can restrict the normalized set to feasible members | ## DTLZ2 gap history @@ -72,4 +84,6 @@ ZDT2 **gens=100** is an early-stress snapshot (collapse on Bend, C#, and pymoo), |-------|-------------------|-----------------| | Pre-fix (wrong ASF) | 0.017 | ~5× | | ASF + persistent ideal | 0.0052 | ~1.5× | -| + duplicate elimination | **0.0040** | **~1.15×** | +| + duplicate elimination | **0.0040** | **~1.15× vs pymoo n_var=10** | + +The ~1.15× denominator is pymoo at **n_var=10** (k=8), not the C# problem (n_var=12, k=10). A matched seed-1 pair is recorded in [ORACLE-RESULTS.md](ORACLE-RESULTS.md). The 15-seed table is still the mismatched-k run. diff --git a/docs/ORACLE-RESULTS.md b/docs/ORACLE-RESULTS.md index ff46d94..8351c11 100644 --- a/docs/ORACLE-RESULTS.md +++ b/docs/ORACLE-RESULTS.md @@ -22,6 +22,8 @@ pip install pymoo python tools/oracle/run_pymoo_oracle.py --problem zdt1 --partitions 12 --pop 52 --gens 100 --seed 1 python tools/oracle/run_pymoo_oracle.py --problem zdt2 --partitions 12 --pop 52 --seed 1 python tools/oracle/run_pymoo_oracle.py --problem dtlz2 --partitions 12 --pop 92 --gens 150 --seed 1 +# DTLZ2 default is n_var=12 (k=10), matching Dtlz2Problem. pymoo's own default is n_var=10 (k=8). +# Pass --n-var 10 only to reproduce the historical mismatched column. # C# dotnet run --project tools/OracleCompare -c Release -- --problem zdt1 --partitions 12 --pop 52 --gens 100 --seed 1 @@ -39,8 +41,27 @@ C# never published a hard ZDT2 oracle / Wilcoxon table. The unpublished Wilcoxon | Problem | Settings | pymoo IGD | C# default IGD | C# `PymooCompatible` IGD | Verdict | |---------|----------|-----------|----------------|--------------------------|---------| -| **ZDT1** | p=12, pop=52, 100 gen | **0.0629** (n=13 ND) | **0.0514** (n=52) | — | **Default wins** | -| **DTLZ2** | p=12, pop=92, 150 gen | **0.00350** (n=91) | 0.0070 (n=92) | **0.00403** (n=92) | **~1.15× pymoo** (pymoo-mode) | +| **ZDT1** | p=12, pop=52, 100 gen | **0.0629** (`res.F`, n=13) | **0.0514** (full ND front, n=52) | — | Different sets. Not an algorithm ranking. | +| **DTLZ2** | p=12, pop=92, 150 gen, **mismatched k** | **0.00350** (n=91, pymoo **n_var=10**, k=8) | 0.0070 (n=92, n_var=12) | **0.00403** (n=92, n_var=12) | Historical pair only. Not a same-problem ratio. | + +### ZDT fronts (same run, different sets) + +C# `OracleCompare` scores the **full feasible non-dominated front** (here n=52) against `ParetoFronts.Zdt1(500)`. pymoo's oracle scores **`res.F`**, the survival niche set (here n=13, one per Das–Dennis direction), against pymoo's 100-point `pareto_front()`. The published 0.0514 vs 0.0629 pair is those two reporters. It is not evidence that the algorithm is better by ~0.011 IGD. + +Seed 1 remeasured **2026-09-22**, pymoo 0.6.2. The C# console reprinted the published scalar. + +| Set | Reference front | n | IGD | +|-----|-----------------|--:|----:| +| C# non-dominated front | library 500-point ZDT1 | 52 | 0.051430749249856716 (console 0.0514307) | +| C# non-dominated front | pymoo 100-point PF | 52 | 0.05119280568479224 | +| pymoo final population, non-dominated | pymoo 100-point PF | 52 | 0.05378307132263516 | +| pymoo `res.F` | pymoo 100-point PF | 13 | 0.0628633417931784 | +| pymoo `res.F` | library 500-point ZDT1 | 13 | 0.06276449352608372 | +| pymoo population ND | library 500-point ZDT1 | 52 | 0.05381261662420749 | + +On the shared 100-point PF, the full-front pair is C# 0.05119280568479224 and pymoo 0.05378307132263516. One seed cannot carry a ranking. Switching the C# front from the 500-point sampler to that 100-point PF changes its IGD by 0.051430749249856716 − 0.05119280568479224 = 2.37943565064476×10⁻⁴, which is much smaller than the 13-versus-52 gap on pymoo's own PF (0.0628633417931784 − 0.05378307132263516 = 0.00908027047054324). + +`ReferenceDirectionThinning.OnePerDirection` keeps the raw objective vector closest (perpendicular distance) to each Das–Dennis direction. On this C# front that helper kept 13 points and scored **0.06357076535717451** against the 100-point PF. That set is not `res.F`. The 15-seed pymoo column is still `res.F`, so its median ratio inherits the same asymmetry. Do not rewrite [WILCOXON-RESULTS.md](WILCOXON-RESULTS.md) until those seeds are re-run on a shared front definition. ### DTLZ2 multi-seed (C# `PymooCompatible`, same protocol) @@ -53,7 +74,27 @@ C# never published a hard ZDT2 oracle / Wilcoxon table. The unpublished Wilcoxon | 5 | 0.00466 | | **mean** | **~0.00485** | -All seeds stay in the same band as pymoo’s single-seed 0.0035 (within ~1.4–1.6×). +Those five C# seeds are `Dtlz2Problem(k: 10)` (n_var=12). The 0.0035 figure they were compared with is pymoo at **n_var=10** (k=8). That is not a same-problem band. The 15-seed file is unchanged until a matched re-run (see below). + +### DTLZ2 n_var (known mismatch, seed 1 remeasured) + +`Dtlz2Problem(nObjectives: 3, k: 10)` builds **n = 12**. Deb et al. suggest k = 10. pymoo 0.6.2 `get_problem("dtlz2", n_obj=3)` defaults to **n_var=10** (k = 8). The harness used to omit `n_var`, so the published seed-1 pair and `docs/WILCOXON-RESULTS.md` compare those two dimensions. `tools/oracle/run_pymoo_oracle.py` now passes **n_var=12**. + +Published mismatched seed 1 (already in the Wilcoxon table; not re-interpreted as parity): + +| Solver | n_var | k | IGD | +|--------|------:|--:|----:| +| C# `PymooCompatible` | 12 | 10 | 0.00403168 | +| pymoo default | 10 | 8 | 0.00349879 | + +Seed 1 remeasured **2026-09-22** with pymoo 0.6.2 after the oracle passes `n_var=12`. Console figures are the `G6` print; the second number is the meta-file value. Front sizes are what each reporter wrote (`res.F` vs full non-dominated front). + +| Solver | n_var | k | Console IGD | Meta IGD | Front | +|--------|------:|--:|------------:|---------:|------:| +| C# `PymooCompatible` | 12 | 10 | 0.00403168 | 0.004031675764658275 | 92 | +| pymoo `n_var=12` | 12 | 10 | 0.00308392 | 0.003083921253245871 | 91 | + +Ratio of the two meta IGDs: 0.004031675764658275 / 0.003083921253245871 = **1.30732**. That is one seed, and the fronts still differ by one point (92 vs 91). It is not a 15-seed ranking and it does not replace the Wilcoxon table. ## Root cause of the old ~5× DTLZ2 gap (fixed) @@ -78,7 +119,7 @@ Deep-dive vs pymoo `HyperplaneNormalization` / `ReferenceDirectionSurvival` (pym |------|-----| | ZDT1 seed=1, 100 gen, default tournament | IGD ≤ 1.5 × 0.0629 | | ZDT2 seed=2, 250 gen, default `RankNicheDistance` | IGD < 0.75 (loose CI smoke, not oracle parity) | -| DTLZ2 seed=1, 150 gen, pymoo-mode | IGD ≤ 3 × 0.00350 (currently ~1.15×) | +| DTLZ2 seed=1, 150 gen, pymoo-mode | IGD ≤ 2 × 0.00350. The 0.00350 scalar is the mismatched n_var=10 run. 3× still passed a regression to about 2.9×. This bar does not claim same-problem equivalence. | | DTLZ2 short smoke (80 gen) | IGD < 0.15 | ZDT2 quality A/B is gens=250 + `PymooCompatible` (not the loose smoke bar). ZDT1 / DTLZ2 shipping bars are unchanged. @@ -88,7 +129,10 @@ ZDT2 quality A/B is gens=250 + `PymooCompatible` (not the loose smoke bar). ZDT1 | Item | Status | |------|--------| | IGD mean-distance | **aligned** | -| ASF / hyperplane normalization | **aligned** | +| ASF / axis intercepts | **aligned** | +| Collapsed nadir (span ≤ 1e-6) | **delta**: nadir = ideal + 1 after the worst-of-pop fallback. pymoo 0.6.2 stops at worst-of-pop. Locked by `Collapsed_span_sets_nadir_to_ideal_plus_one` (`{2, 2+1e-8}` → nadir 3). | +| Infeasible points | **current rule, locked**: ideal and worst include them. Fixture is feasible (1, 1) vs infeasible (0, 0) → ideal (0, 0). No constrained benchmark yet. | +| Mating re-association | **delta**: after survival, `PrepareForSelection` normalizes the survivors again and overwrites niche ids. pymoo keeps the survival ids. Fixture: `PrepareForSelection_overwrites_survival_niche_ids`. | | Persistent ideal + ND extremes | **aligned** | | `TournamentMode.PymooCompatible` | **implemented** | | Duplicate elimination | **implemented** (default on) | diff --git a/docs/RESEARCH-STANDARDS.md b/docs/RESEARCH-STANDARDS.md index 967dad7..cd59a68 100644 --- a/docs/RESEARCH-STANDARDS.md +++ b/docs/RESEARCH-STANDARDS.md @@ -26,27 +26,28 @@ Sources consulted (2026-08): Default dimensions (Deb / pymoo convention): - ZDT1–3: n=30; ZDT4: n=10; ZDT6: n=10 -- DTLZ: n = M + k − 1 with k=5 (DTLZ1) or k=10 (DTLZ2–4), k=20 (DTLZ7) +- DTLZ: n = M + k − 1 with k=5 (DTLZ1) or k=10 (DTLZ2–4), k=20 (DTLZ7) +- C# `Dtlz2Problem` uses that k=10, so M=3 ⇒ **n_var=12**. pymoo 0.6.2 `get_problem("dtlz2", n_obj=3)` defaults to **n_var=10** (k=8). Oracle runs pass `n_var=12`. The published 15-seed pymoo column is the default-10 run and is not a same-k comparison. ## 2. Algorithm hyperparameters (match paper + pymoo) | Knob | Standard value | |------|----------------| -| Crossover | SBX, η_c = **30**, p_c = 1.0 | +| Crossover | SBX, η_c = **30**, p_c = **1.0** (pymoo). Paper section 4 uses p_c = **0.9** | | Mutation | Polynomial, η_m = **20**, p_m = **1/n** | -| Reference set | **Das–Dennis** (uniform) on unit simplex | -| Population size | Often = #reference directions (or slightly larger) | +| Reference set | **Das–Dennis** (uniform) on the unit simplex, **single layer**. Two-layer directions for larger M are absent | +| Population size | Often = #reference directions (or slightly larger). N ≥ 2. `WithDasDennis(1, 1)` throws because |H| = 1; pass an explicit population size for single-objective runs | | Selection | U-NSGA-III **tournament** (not NSGA-III random mating) | ### Tournament detail (alignment note) -**pymoo** `comp_by_rank_and_ref_line_dist`: +**Seada & Deb Algorithm 2** (feasible parents): if both are associated with the same reference direction, prefer rank, then perpendicular distance; otherwise pick at random. If either parent is infeasible, use the constraint comparison. On a distance tie the paper keeps the second parent. The mating pool is two shuffled passes of consecutive pairs. Section 4 uses SBX with p_c = 0.9. -1. If either infeasible → smaller CV wins -2. Else if **same niche** → better rank, else smaller distance-to-niche -3. Else → random +**pymoo** `comp_by_rank_and_ref_line_dist` follows that same-niche / different-niche split and coin-flips a distance tie. pymoo `NSGA3` builds `SBX(eta=30, prob=1.0)`. -**This library (v0.1)** prefers rank → niche count → perpendicular distance (Seada-style pressure even across niches). Documented difference for equivalence work; a `PymooCompatibleTournament` mode can be added if bit-identical mating is required. +**`TournamentMode.PymooCompatible`** matches that pymoo comparator, including the coin flip. It does not use the paper's second-parent tie break, p_c = 0.9, or the consecutive-pair mating pool. + +**`TournamentMode.RankNicheDistance`** is the constructor default: rank, then niche count, then perpendicular distance, including when the niches differ. That is a local expansion, not Algorithm 2. The default stays `RankNicheDistance`. ZDT1's published Wilcoxon table uses it. ## 3. Performance indicators (what to report) @@ -61,8 +62,8 @@ Definitions implemented in `Unsga3.Metrics.PerformanceIndicators` follow **pymoo ### Reference fronts -- ZDT1/2/4/6: closed form f₂(f₁) -- ZDT3: known f₁ intervals +- ZDT1/2/4/6: closed form f₂(f₁). `Zdt1(1)` (and the other one-point samplers) throw. ZDT6's sampler starts at the truncated floor 0.280775, slightly below the minimized f1. +- ZDT3: known f₁ intervals. The second left endpoint in this library is 0.1822287280; pymoo 0.6.2 writes 0.182228780. The library literal is locked. - DTLZ1: Das–Dennis × 0.5 on simplex - DTLZ2/3/4: Das–Dennis projected to unit sphere @@ -72,7 +73,7 @@ Sample **≥ 500** points on continuous bi-objective fronts (common practice). Typical ZDT: r = (1.1, 1.1). Always document r; never compare HV across different r. -## 4. Equivalence protocol (one-to-one claim) +## 4. Equivalence protocol 1. Same problem definition (bounds, n, evaluate) 2. Same Das–Dennis partitions → identical ref set size @@ -81,8 +82,11 @@ Typical ZDT: r = (1.1, 1.1). Always document r; never compare HV across differen - ZDT2 quality A/B: pop=52, **gens=250**, `PymooCompatible` (matches unsga3-bend). gens=100 is an early-stress snapshot, not the quality bar. `RankNicheDistance` is optional, not the ZDT2 default. - DTLZ2: pop=92, **gens=150**, `PymooCompatible` 4. Fixed seed **or** 15–31 seeds → median + IQR IGD -5. Compare IGD (and HV for M=2) to pymoo `UNSGA3` -6. Shipping bar: median IGD within ~1–2% of pymoo on ZDT1/DTLZ2 (or non-inferior Wilcoxon). ZDT2 has no published C# Wilcoxon table; quality budget is 250 gens. +5. Compare IGD (and HV for M=2) to pymoo `UNSGA3` on the **same front definition**. C# uses the full non-dominated front; pymoo's harness value is `res.F` (the niche optimum). A gap between those two reporters is a set-definition gap until both sides are reduced the same way. +6. Shipping bar: the published 15-seed table is **not** “median IGD within ~1–2% of pymoo.” + ZDT1: Mann–Whitney U = 65, p = 0.0512394, median ratio 0.764107 (pymoo column is `res.F`). + DTLZ2: U = 222, p = 6.15164×10⁻⁶, median ratio 1.58638 (pymoo column is n_var=10). That test rejects equal distributions at α = 0.05. + ZDT2 has no published C# Wilcoxon table; the quality budget is 250 generations. A matched 15-seed re-run has not replaced the table. Export path: dump final `F` as CSV from both sides; compute IGD in this library. diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index c584c2b..7dc1952 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -17,7 +17,7 @@ Living plan for Unsga3. Issues track concrete work; this page is the narrative. - [ ] Multi-seed Wilcoxon results checked in / refreshed on release - [ ] `net8.0` (+ `net10.0`) multi-target for broader NuGet consumers - [ ] nuget.org publish (in addition to GitHub Packages) -- [ ] IGD+ / GD+ indicators +- [ ] GD+ indicator (IGD+ has a hand-case test; GD and IGD were already implemented) - [ ] Constrained demos (OSY / TNK) with self-tests - [ ] API docs site (DocFX or similar) diff --git a/docs/WILCOXON-RESULTS.md b/docs/WILCOXON-RESULTS.md index be2f04c..9d3c57c 100644 --- a/docs/WILCOXON-RESULTS.md +++ b/docs/WILCOXON-RESULTS.md @@ -12,6 +12,8 @@ Generated by `tools/oracle/run_multiseed_wilcoxon.py`. IGD = mean nearest Euclid ZDT2 **gens=100** is an early-stress snapshot, not the quality bar. C# never published a hard ZDT2 Wilcoxon table; the unpublished harness used gens=100 + `RankNicheDistance` (optional mating mode — do not silently switch all ZDT defaults to it). Quality protocol matches [unsga3-bend](https://github.com/AppSprout-dev/unsga3-bend) A/B honesty: gens=250 + `PymooCompatible`. ZDT1 and DTLZ2 numbers below are unchanged. Do not invent a ZDT2 IGD table here. +**DTLZ2 dimension (do not refresh this file in place).** The pymoo column below was produced with pymoo’s default `n_var=10` (k=8). C# used `Dtlz2Problem(k: 10)` (`n_var=12`). Seed 1 of that mismatched pair is C# 0.00403168 / pymoo 0.00349879. New oracle runs pass `n_var=12`. A matched seed-1 pair is in [ORACLE-RESULTS.md](ORACLE-RESULTS.md). **Leave this table as published until the 15 seeds are actually re-run.** + Hypothesis tests (α = 0.05, two-sided): - **Mann–Whitney U** (Wilcoxon rank-sum): independent samples, H₀: same IGD distribution. @@ -101,6 +103,7 @@ Lower IGD is better. ## Notes +- ZDT1 sets differ. The Unsga3 column is the full non-dominated front. The pymoo column is `res.F` (about one point per reference direction). The median ratio **0.764107** inherits that asymmetry. A shared-front seed-1 note is in [ORACLE-RESULTS.md](ORACLE-RESULTS.md). Do not replace this table until both columns use one front definition. - Not bit-identical: different RNG implementations and minor operator ordering. -- Practical equivalence: median IGD within ~1–2× and non-significant MWU is a strong claim; significant differences with small effect size (ratio ≈ 1) are still acceptable for a v0.x port. +- These tests are not a 1–2% equivalence claim. ZDT1 U = 65, p = 0.0512394, median ratio 0.764107, and its pymoo column is `res.F`. DTLZ2 U = 222, p = 6.15164e-06, median ratio 1.58638, and its pymoo column is n_var=10. DTLZ2 rejects equal distributions at α = 0.05. - Reproduce: `python tools/oracle/run_multiseed_wilcoxon.py` diff --git a/src/Unsga3/Algorithm/Unsga3Algorithm.cs b/src/Unsga3/Algorithm/Unsga3Algorithm.cs index 985c225..9d1b627 100644 --- a/src/Unsga3/Algorithm/Unsga3Algorithm.cs +++ b/src/Unsga3/Algorithm/Unsga3Algorithm.cs @@ -27,13 +27,19 @@ public sealed class Unsga3Algorithm /// Defaults to the number of reference directions. /// Defaults to SBX η=30. /// Defaults to polynomial mutation η=20. - /// Probability of applying SBX to a parent pair. + /// + /// Probability of applying SBX to a parent pair. Default 1.0, matching pymoo + /// SBX(prob=1.0). Seada & Deb section 4 uses 0.9. + /// /// Per-variable mutation probability; default 1/nVars at run time. /// Optional RNG seed for reproducibility. /// Mating tournament policy; use for oracle runs. /// /// Drop offspring whose decision vector matches an existing parent or earlier offspring - /// (pymoo eliminate_duplicates=True). Default true. + /// (pymoo eliminate_duplicates=True). Default true. The key is G12 + /// (12 significant digits), not 12 digits after the decimal. If mutation cannot + /// produce a new key, attempts are capped and the remaining slots may be duplicates + /// so the loop cannot hang. /// public Unsga3Algorithm( double[][] referenceDirections, @@ -66,7 +72,13 @@ public Unsga3Algorithm( _eliminateDuplicates = eliminateDuplicates; } - /// Convenience: build Das–Dennis directions then construct the algorithm. + /// + /// Convenience: build Das–Dennis directions then construct the algorithm. + /// One objective produces a single direction, so the default population size is 1. + /// The constructor requires N ≥ 2, and WithDasDennis(1, 1) throws. + /// Single-objective runs must pass ≥ 2 + /// (Seada & Deb recommend a multiple of four, and at least |H|). + /// public static Unsga3Algorithm WithDasDennis( int numberOfObjectives, int partitions, @@ -143,6 +155,8 @@ public OptimizationResult Run( var next = survival.Select(combined, _populationSize, rng); population = new Population(next); + // pymoo keeps the niche ids written during survival. This call normalizes the + // survivors again and overwrites AssociatedReference before the next mating. TournamentSelection.PrepareForSelection(population.Members, refs, normalization); generation++; } @@ -185,16 +199,27 @@ private List CreateOffspring( TryAddOffspring(offspring, c2, seen); } - // Fallback: mutated clones if de-dup exhausted attempts (should be rare). - while (offspring.Count < _populationSize) + // Mutation that cannot change x used to spin here: a duplicate was accepted + // only when one slot remained, and nothing incremented when two or more remained. + int fallbackAttempts = 0; + int fallbackCap = Math.Max(_populationSize * 20, 1); + while (offspring.Count < _populationSize && fallbackAttempts < fallbackCap) { + fallbackAttempts++; var extra = parents[rng.Next(parents.Count)].Clone(); _mutation.Mutate(extra, problem, rng, mutProb); - // Always accept in the hard-fallback path so we never deadlock. - if (seen is null || seen.Add(DecisionKey(extra.Variables)) || offspring.Count + 1 >= _populationSize) + if (seen is null || seen.Add(DecisionKey(extra.Variables))) offspring.Add(extra); } + // Last resort: accept duplicates so elimination cannot hang. + while (offspring.Count < _populationSize) + { + var extra = parents[rng.Next(parents.Count)].Clone(); + _mutation.Mutate(extra, problem, rng, mutProb); + offspring.Add(extra); + } + return offspring; } @@ -209,10 +234,12 @@ private static void TryAddOffspring(List offspring, Individual child offspring.Add(child); } - /// Stable decision-vector key for duplicate elimination (rounded to 12 dp). - private static string DecisionKey(double[] x) + /// + /// Decision-vector key for duplicate elimination. G12 is 12 significant digits, + /// not 12 digits after the decimal point. + /// + internal static string DecisionKey(double[] x) { - // Invariant culture, fixed decimals — enough for continuous SBX without false collisions. var sb = new System.Text.StringBuilder(x.Length * 18); for (int i = 0; i < x.Length; i++) { diff --git a/src/Unsga3/Core/Normalization.cs b/src/Unsga3/Core/Normalization.cs index dbb1f47..222eeb6 100644 --- a/src/Unsga3/Core/Normalization.cs +++ b/src/Unsga3/Core/Normalization.cs @@ -4,10 +4,24 @@ namespace Unsga3.Core; /// /// NSGA-III adaptive hyperplane normalization (Deb & Jain), aligned with pymoo -/// HyperplaneNormalization: +/// HyperplaneNormalization on the intercept path: /// persistent ideal / worst points, ASF extreme points (optionally from the ND front), /// intercept-based nadir with front/population fallbacks. /// +/// +/// Collapsed span is a known delta versus pymoo 0.6.2. Both sides fall back to the +/// worst point in the population when the nadir span is at most 1e-6. If that span +/// is still at most 1e-6, this library sets nadir = ideal + 1. pymoo stops at the +/// worst-of-population, so two nearly equal objectives stay a tiny span apart and +/// normalize to 0 and 1. Here they normalize to about 0 and 1e-8. See +/// NormalizationTests.Collapsed_span_sets_nadir_to_ideal_plus_one. +/// Ideal and worst are updated from every point in the pool, feasible or not. +/// There is no constrained benchmark in this library. A feasible (1, 1) beside an +/// infeasible (0, 0) therefore takes ideal (0, 0) from the infeasible point. +/// Extreme-point ASF uses the non-dominated index set when the caller supplies one +/// (constraint-domination puts only the feasible point on that front) and the whole +/// pool when it does not. See Infeasible_origin_sets_ideal_from_the_whole_pool. +/// public sealed class Normalization { private readonly int _m; @@ -197,7 +211,8 @@ private void UpdateNadir(IReadOnlyList population, int[] ndIdx) _nadir[j] = worstOfFront[j]; } - // Degenerate range → fall back to worst of population. + // Degenerate range → worst of this population, then ideal+1. + // pymoo 0.6.2 stops after the worst-of-population assignment. for (int j = 0; j < _m; j++) { if (_nadir[j] - _ideal[j] <= 1e-6) diff --git a/src/Unsga3/Metrics/ParetoFronts.cs b/src/Unsga3/Metrics/ParetoFronts.cs index e6a6c45..9150847 100644 --- a/src/Unsga3/Metrics/ParetoFronts.cs +++ b/src/Unsga3/Metrics/ParetoFronts.cs @@ -11,6 +11,9 @@ public static class ParetoFronts /// ZDT1: f2 = 1 - sqrt(f1), f1 ∈ [0,1]. public static double[][] Zdt1(int nPoints = 500) { + if (nPoints < 2) + throw new ArgumentOutOfRangeException(nameof(nPoints), "Need at least 2 points."); + var pf = new double[nPoints][]; for (int i = 0; i < nPoints; i++) { @@ -23,6 +26,9 @@ public static double[][] Zdt1(int nPoints = 500) /// ZDT2: f2 = 1 - f1². public static double[][] Zdt2(int nPoints = 500) { + if (nPoints < 2) + throw new ArgumentOutOfRangeException(nameof(nPoints), "Need at least 2 points."); + var pf = new double[nPoints][]; for (int i = 0; i < nPoints; i++) { @@ -37,7 +43,13 @@ public static double[][] Zdt2(int nPoints = 500) /// public static double[][] Zdt3(int pointsPerSegment = 100) { - // Known f1 intervals for ZDT3 Pareto set (Deb). + if (pointsPerSegment < 2) + throw new ArgumentOutOfRangeException(nameof(pointsPerSegment), "Need at least 2 points per segment."); + + // Known f1 intervals for the ZDT3 Pareto set. + // The second left endpoint is 0.1822287280. pymoo 0.6.2 writes 0.182228780 + // in zdt.py; the two literals differ at the eighth significant digit. + // This value stays as transcribed. MetricsTests locks it. double[][] intervals = { new[] { 0.0, 0.0830015349 }, @@ -63,10 +75,19 @@ public static double[][] Zdt3(int pointsPerSegment = 100) /// ZDT4 same geometry as ZDT1. public static double[][] Zdt4(int nPoints = 500) => Zdt1(nPoints); - /// ZDT6: f1 from ~0.280775 to 1, f2 = 1 - f1². + /// + /// ZDT6: f2 = 1 - f1², with f1 sampled from 0.280775 up to 1. + /// 0.280775 is a six-digit truncation of the minimized + /// f1(x) = 1 - exp(-4x) sin⁶(6πx). A uniform grid of 2,000,001 points on [0, 1] + /// found a minimum of 0.280775318847039 (x = 0.081458), which is 3.188×10⁻⁷ above + /// this floor, so the first sample sits slightly below the true front. + /// public static double[][] Zdt6(int nPoints = 500) { - // f1* = 1 - exp(-4x) sin^6(6πx) for x in [0,1]; min ≈ 0.280775 + if (nPoints < 2) + throw new ArgumentOutOfRangeException(nameof(nPoints), "Need at least 2 points."); + + // Truncated floor. See the summary comment. MetricsTests locks 0.280775. double f1Min = 0.280775; var pf = new double[nPoints][]; for (int i = 0; i < nPoints; i++) diff --git a/src/Unsga3/Metrics/PerformanceIndicators.cs b/src/Unsga3/Metrics/PerformanceIndicators.cs index b5b0fa6..d960e7e 100644 --- a/src/Unsga3/Metrics/PerformanceIndicators.cs +++ b/src/Unsga3/Metrics/PerformanceIndicators.cs @@ -116,6 +116,9 @@ public static double Hypervolume2D(IReadOnlyList front, double[] refer private static double ModifiedDistance(double[] a, double[] z) { + if (a.Length != z.Length) + throw new ArgumentException("Objective vectors must have the same length."); + // d+ from z toward a for minimization: Euclidean of max(a_j - z_j, 0) double s = 0; for (int k = 0; k < z.Length; k++) @@ -140,8 +143,11 @@ private static double NearestDistance(double[] point, IReadOnlyList se private static double Euclidean(double[] a, double[] b) { + if (a.Length != b.Length) + throw new ArgumentException("Objective vectors must have the same length."); + double s = 0; - int n = Math.Min(a.Length, b.Length); + int n = a.Length; for (int i = 0; i < n; i++) { double d = a[i] - b[i]; diff --git a/src/Unsga3/Metrics/ReferenceDirectionThinning.cs b/src/Unsga3/Metrics/ReferenceDirectionThinning.cs new file mode 100644 index 0000000..33a30ca --- /dev/null +++ b/src/Unsga3/Metrics/ReferenceDirectionThinning.cs @@ -0,0 +1,78 @@ +using Unsga3.Algorithm; + +namespace Unsga3.Metrics; + +/// +/// Optional comparison aid. pymoo UNSGA3 reports res.F as the survival +/// niche set (about one member per filled reference direction). This library scores +/// the full non-dominated front. keeps the raw objective +/// vector with the smallest perpendicular distance to each direction so a caller can +/// score a similar cardinality. +/// +/// +/// The result is not pymoo res.F. It does not repeat hyperplane normalization +/// or survival, and it is not evidence that the two algorithms match. +/// +public static class ReferenceDirectionThinning +{ + /// + /// Keep at most one point per reference direction: the member with the smallest + /// perpendicular distance to that direction. Directions with no assigned point are omitted. + /// + public static double[][] OnePerDirection( + IReadOnlyList front, + IReadOnlyList directions) + { + ArgumentNullException.ThrowIfNull(front); + ArgumentNullException.ThrowIfNull(directions); + if (directions.Count == 0) + throw new ArgumentException("Need at least one direction.", nameof(directions)); + + int m = directions[0].Length; + for (int r = 1; r < directions.Count; r++) + { + if (directions[r].Length != m) + throw new ArgumentException("All directions must have the same length.", nameof(directions)); + } + + var bestDist = new double[directions.Count]; + var bestIdx = new int[directions.Count]; + Array.Fill(bestDist, double.PositiveInfinity); + Array.Fill(bestIdx, -1); + + for (int i = 0; i < front.Count; i++) + { + var f = front[i]; + if (f.Length != m) + throw new ArgumentException( + "Each front point must have the same length as the reference directions.", + nameof(front)); + + int bestRef = 0; + double best = double.PositiveInfinity; + for (int r = 0; r < directions.Count; r++) + { + double d = ReferencePointManager.PerpendicularDistance(f, directions[r]); + if (d < best) + { + best = d; + bestRef = r; + } + } + + if (best < bestDist[bestRef]) + { + bestDist[bestRef] = best; + bestIdx[bestRef] = i; + } + } + + var kept = new List(); + for (int r = 0; r < bestIdx.Length; r++) + { + if (bestIdx[r] >= 0) + kept.Add(front[bestIdx[r]]); + } + return kept.ToArray(); + } +} diff --git a/src/Unsga3/Operators/Selection/TournamentMode.cs b/src/Unsga3/Operators/Selection/TournamentMode.cs index 75a6eef..ff3522b 100644 --- a/src/Unsga3/Operators/Selection/TournamentMode.cs +++ b/src/Unsga3/Operators/Selection/TournamentMode.cs @@ -4,15 +4,18 @@ namespace Unsga3.Operators.Selection; public enum TournamentMode { /// - /// Rank → niche count → perpendicular distance (default). - /// Stronger selection pressure across niches than stock pymoo. + /// Rank, then niche count, then perpendicular distance (constructor default). + /// Niche count is compared even when the two parents sit on different reference + /// directions. That is not Seada & Deb Algorithm 2, which picks at random + /// across directions. The default is intentionally unchanged. /// RankNicheDistance = 0, /// - /// Matches pymoo comp_by_rank_and_ref_line_dist: - /// CV first; if same niche then rank then dist-to-niche; else random. - /// Use for oracle / equivalence runs against pymoo. + /// Same-niche / different-niche split from Seada & Deb Algorithm 2 and from + /// pymoo comp_by_rank_and_ref_line_dist: constraint violation first; if the + /// parents share a niche then rank, then distance-to-niche; otherwise random. + /// A distance tie is a coin flip (pymoo). Algorithm 2 keeps the second parent on that tie. /// PymooCompatible = 1, } diff --git a/src/Unsga3/Operators/Selection/TournamentSelection.cs b/src/Unsga3/Operators/Selection/TournamentSelection.cs index ca478ea..7c5ad4a 100644 --- a/src/Unsga3/Operators/Selection/TournamentSelection.cs +++ b/src/Unsga3/Operators/Selection/TournamentSelection.cs @@ -5,7 +5,10 @@ namespace Unsga3.Operators.Selection; /// -/// U-NSGA-III niching-based binary tournament (Seada & Deb / pymoo variants). +/// Binary mating tournament. follows the +/// Seada & Deb Algorithm 2 split (same niche: rank then distance; different niches: random) +/// and pymoo's coin flip on equal distance. +/// is the constructor default and also prefers the smaller niche count across niches. /// public sealed class TournamentSelection { @@ -17,7 +20,8 @@ public TournamentSelection(TournamentMode mode = TournamentMode.RankNicheDistanc public TournamentMode Mode { get; } /// - /// Select parents (with replacement tournaments) from the population. + /// Select parents by independent tournaments with replacement. + /// This is not the paper's two shuffled passes of consecutive pairs. /// Population must already have Rank / niche association set via . /// public List SelectParents( @@ -112,7 +116,10 @@ internal static Individual WinnerPymoo(Individual a, Individual b, RandomProvide return rng.NextDouble() < 0.5 ? a : b; } - /// Recompute ranks + niche counts for tournament (normalize + associate). + /// + /// Recompute ranks and niche association for mating. This is a second normalization + /// of the survivors. pymoo keeps the niche ids from environmental selection. + /// public static void PrepareForSelection( IReadOnlyList population, ReferencePointManager references, diff --git a/src/Unsga3/Problems/DtlzProblems.cs b/src/Unsga3/Problems/DtlzProblems.cs index c41c53c..a35bffc 100644 --- a/src/Unsga3/Problems/DtlzProblems.cs +++ b/src/Unsga3/Problems/DtlzProblems.cs @@ -57,7 +57,12 @@ protected override void EvaluateCore(double[] x, double[] f, double[] g) } } -/// DTLZ2 — unit sphere (first orthant). +/// +/// DTLZ2 — unit sphere (first orthant). +/// Default k = 10 gives n = M + k − 1 (12 when M = 3), Deb's suggested k. +/// pymoo get_problem("dtlz2", n_obj=3) defaults to n_var = 10 (k = 8). +/// That default is a different search problem. The oracle passes n_var = 12. +/// public sealed class Dtlz2Problem : ProblemBase { public Dtlz2Problem(int nObjectives = 3, int k = 10) diff --git a/src/Unsga3/Utilities/DasDennis.cs b/src/Unsga3/Utilities/DasDennis.cs index 8250c30..568c095 100644 --- a/src/Unsga3/Utilities/DasDennis.cs +++ b/src/Unsga3/Utilities/DasDennis.cs @@ -2,6 +2,8 @@ namespace Unsga3.Utilities; /// /// Das–Dennis structured reference directions on the unit simplex (NSGA-III / U-NSGA-III). +/// Single layer only. Two-layer directions (an outer layer plus an inside layer), +/// which many-objective NSGA-III uses for larger M, are absent. /// public static class ReferenceDirections { @@ -39,11 +41,16 @@ public static int PartitionsForMinimumDirections(int numberOfObjectives, int min return p; } - /// Number of Das–Dennis points: C(p + M - 1, M - 1). + /// + /// Number of Das–Dennis points: C(p + M - 1, M - 1). + /// The combination is computed in a and checked into . + /// is thrown when the value does not fit in + /// (for example M = 11, p = 34, C(44, 10) = 2,481,256,778). + /// public static int Count(int numberOfObjectives, int partitions) { if (numberOfObjectives == 1) return 1; - return Binomial(partitions + numberOfObjectives - 1, numberOfObjectives - 1); + return checked((int)Binomial(partitions + numberOfObjectives - 1, numberOfObjectives - 1)); } private static void Recurse(List points, double[] current, int m, int p, int left, int index) @@ -62,7 +69,7 @@ private static void Recurse(List points, double[] current, int m, int } } - private static int Binomial(int n, int k) + private static long Binomial(int n, int k) { if (k < 0 || k > n) return 0; if (k == 0 || k == n) return 1; @@ -70,9 +77,9 @@ private static int Binomial(int n, int k) long result = 1; for (int i = 1; i <= k; i++) { - result *= n - k + i; + result = checked(result * (n - k + i)); result /= i; } - return (int)result; + return result; } } diff --git a/tests/Unsga3.Tests/Benchmarks/IgdSmokeTests.cs b/tests/Unsga3.Tests/Benchmarks/IgdSmokeTests.cs index 673b73f..9b5214a 100644 --- a/tests/Unsga3.Tests/Benchmarks/IgdSmokeTests.cs +++ b/tests/Unsga3.Tests/Benchmarks/IgdSmokeTests.cs @@ -97,9 +97,10 @@ public void Dtlz2_within_factor_of_pymoo_oracle() var obtained = result.NonDominatedSolutions.Select(i => (double[])i.Objectives.Clone()).ToArray(); double igd = PerformanceIndicators.InvertedGenerationalDistance(obtained, ParetoFronts.Dtlz2(3, 12)); const double pymooBaseline = 0.00350; - // ~3× still tracks parity work; was ~5× (0.017) pre-fix and ~10× (0.037) on default. - Assert.True(igd <= pymooBaseline * 3.0, - $"DTLZ2 IGD={igd} should be ≤ 3× pymoo baseline {pymooBaseline}"); + // Published mismatched-k scalar (pymoo n_var=10). C# seed 1 is about 1.15× this. + // 3× still passed a regression to about 2.9×. 2× rejects that and keeps the current run. + Assert.True(igd <= pymooBaseline * 2.0, + $"DTLZ2 IGD={igd} should be ≤ 2× pymoo baseline {pymooBaseline}"); } [Fact] diff --git a/tests/Unsga3.Tests/Equivalence/EquivalencePlaceholderTests.cs b/tests/Unsga3.Tests/Equivalence/EquivalencePlaceholderTests.cs index 0b6dc2f..a203bc3 100644 --- a/tests/Unsga3.Tests/Equivalence/EquivalencePlaceholderTests.cs +++ b/tests/Unsga3.Tests/Equivalence/EquivalencePlaceholderTests.cs @@ -1,15 +1,37 @@ +using Unsga3.Algorithm; +using Unsga3.Metrics; +using Unsga3.Operators.Selection; +using Unsga3.Problems; +using Unsga3.Utilities; + namespace Unsga3.Tests.Equivalence; /// -/// Placeholder for pymoo / MATLAB one-to-one equivalence (IGD/HV, fixed seeds). -/// Wire Python.NET or export CSV populations here once the oracle harness lands. +/// Seed-1 DTLZ2 regression guard. This used to assert true, so a broken run still passed. +/// The published mismatched-k pymoo scalar is 0.00350 and the C# seed-1 IGD is about 0.00403 +/// (~1.15× that scalar). A bar of 3× that scalar still passes a regression to about 2.9×. /// public class EquivalencePlaceholderTests { [Fact] - public void Harness_not_yet_wired() + public void Dtlz2_seed1_rejects_a_near_3x_regression() { - // Intentional: documents the planned equivalence suite (docs/EQUIVALENCE.md). - Assert.True(true); + // 0.00350 is the published pymoo seed-1 scalar at n_var=10 (k=8), not the matched + // n_var=12 problem. Passing this test is a regression guard, not a same-problem claim. + const double publishedMismatchedPymoo = 0.00350; + const double maxFactor = 2.0; + + var problem = new Dtlz2Problem(nObjectives: 3, k: 10); + var dirs = ReferenceDirections.DasDennis(3, 12); + var algo = new Unsga3Algorithm( + dirs, populationSize: 92, seed: 1, tournamentMode: TournamentMode.PymooCompatible); + var result = algo.Run(problem, maxGenerations: 150); + var obtained = result.NonDominatedSolutions.Select(i => (double[])i.Objectives.Clone()).ToArray(); + double igd = PerformanceIndicators.InvertedGenerationalDistance(obtained, ParetoFronts.Dtlz2(3, 12)); + + Assert.True( + igd <= publishedMismatchedPymoo * maxFactor, + $"DTLZ2 IGD={igd} exceeds {maxFactor}× {publishedMismatchedPymoo}. " + + "A factor of 3 would still pass a regression to about 2.9×."); } } diff --git a/tests/Unsga3.Tests/Unit/DasDennisTests.cs b/tests/Unsga3.Tests/Unit/DasDennisTests.cs index cd5ec17..654c4e7 100644 --- a/tests/Unsga3.Tests/Unit/DasDennisTests.cs +++ b/tests/Unsga3.Tests/Unit/DasDennisTests.cs @@ -29,6 +29,13 @@ public void Points_lie_on_unit_simplex() } } + [Fact] + public void Count_throws_when_the_combination_exceeds_int32() + { + // C(34 + 11 - 1, 10) = C(44, 10) = 2_481_256_778, which does not fit in Int32. + Assert.Throws(() => ReferenceDirections.Count(11, 34)); + } + [Fact] public void Single_objective_is_unit_scalar() { diff --git a/tests/Unsga3.Tests/Unit/DuplicateEliminationTests.cs b/tests/Unsga3.Tests/Unit/DuplicateEliminationTests.cs new file mode 100644 index 0000000..586c455 --- /dev/null +++ b/tests/Unsga3.Tests/Unit/DuplicateEliminationTests.cs @@ -0,0 +1,66 @@ +using System.Globalization; +using Unsga3.Algorithm; +using Unsga3.Core; +using Unsga3.Operators.Crossover; +using Unsga3.Operators.Mutation; +using Unsga3.Utilities; + +namespace Unsga3.Tests.Unit; + +public class DuplicateEliminationTests +{ + [Fact] + public void Decision_key_uses_G12_significant_digits() + { + double tiny = 1e-20; + string key = Unsga3Algorithm.DecisionKey(new[] { tiny }); + string significant = tiny.ToString("G12", CultureInfo.InvariantCulture); + string twelveDecimalPlaces = tiny.ToString("F12", CultureInfo.InvariantCulture); + + Assert.Equal(significant, key); + Assert.Equal("0.000000000000", twelveDecimalPlaces); + Assert.NotEqual(twelveDecimalPlaces, key); + } + + [Fact] + public void Elimination_returns_when_mutation_cannot_diversify() + { + var dirs = ReferenceDirections.DasDennis(1, 1); + var algo = new Unsga3Algorithm( + dirs, + populationSize: 6, + crossover: new CopyParentsCrossover(), + mutation: new NoopMutation(), + seed: 1, + eliminateDuplicates: true); + + var result = algo.Run(new IdentityProblem(), maxGenerations: 1); + Assert.Equal(6, result.FinalPopulation.Count); + } + + private sealed class IdentityProblem : ProblemBase + { + public IdentityProblem() + : base(1, 1, 0, new[] { (0.0, 1.0) }) + { + } + + protected override void EvaluateCore(double[] x, double[] f, double[] g) => f[0] = x[0]; + } + + private sealed class CopyParentsCrossover : ICrossover + { + public (Individual Child1, Individual Child2) Crossover( + Individual parent1, Individual parent2, IProblem problem, RandomProvider rng) + { + return (parent1.Clone(), parent2.Clone()); + } + } + + private sealed class NoopMutation : IMutation + { + public void Mutate(Individual individual, IProblem problem, RandomProvider rng, double probabilityPerVariable) + { + } + } +} diff --git a/tests/Unsga3.Tests/Unit/MatingRenormalizeTests.cs b/tests/Unsga3.Tests/Unit/MatingRenormalizeTests.cs new file mode 100644 index 0000000..8cec062 --- /dev/null +++ b/tests/Unsga3.Tests/Unit/MatingRenormalizeTests.cs @@ -0,0 +1,66 @@ +using Unsga3.Algorithm; +using Unsga3.Core; +using Unsga3.Operators.Selection; +using Unsga3.Operators.Survival; +using Unsga3.Utilities; + +namespace Unsga3.Tests.Unit; + +/// +/// After survival, normalizes the +/// survivors again and overwrites niche ids. pymoo keeps the ids written during survival. +/// +public class MatingRenormalizeTests +{ + [Fact] + public void PrepareForSelection_overwrites_survival_niche_ids() + { + var dirs = ReferenceDirections.DasDennis(2, 1); + var refs = new ReferencePointManager(dirs); + var norm = new Normalization(2); + var survival = new NondominatedSortingSurvival(refs, norm); + + var pool = new List + { + Point(0.0, 1.0), + Point(1.0, 0.0), + Point(0.2, 0.9), + Point(0.9, 0.2), + Point(0.4, 0.4), + }; + + // targetSize < |ND front| forces last-front niching, which writes niche ids. + var survivors = survival.Select(pool, targetSize: 3, rng: null); + int[] survivalNiches = survivors.Select(s => s.AssociatedReference).ToArray(); + Assert.All(survivalNiches, id => Assert.InRange(id, 0, dirs.Length - 1)); + + foreach (var s in survivors) + { + s.AssociatedReference = 999; + s.PerpendicularDistance = -1; + } + + TournamentSelection.PrepareForSelection(survivors, refs, norm); + int[] matingNiches = survivors.Select(s => s.AssociatedReference).ToArray(); + + Assert.DoesNotContain(999, matingNiches); + Assert.All(survivors, s => Assert.True(s.PerpendicularDistance >= 0)); + + // Locked dump for this pool (Das–Dennis p=1, target 3, no RNG). + // The ids match, and the sentinel 999 is gone, so mating rewrote the fields + // pymoo would have kept. Matching ids here is not a claim that every generation matches. + Assert.Equal(new[] { 0, 1, 0 }, survivalNiches); + Assert.Equal(new[] { 0, 1, 0 }, matingNiches); + Assert.Equal(0.0, survivors[0].PerpendicularDistance, 9); + Assert.Equal(0.0, survivors[1].PerpendicularDistance, 9); + Assert.Equal(0.2, survivors[2].PerpendicularDistance, 9); + } + + private static Individual Point(double f1, double f2) + { + var ind = new Individual(1, 2); + ind.Objectives[0] = f1; + ind.Objectives[1] = f2; + return ind; + } +} diff --git a/tests/Unsga3.Tests/Unit/MetricsTests.cs b/tests/Unsga3.Tests/Unit/MetricsTests.cs index ad5bf93..739fadc 100644 --- a/tests/Unsga3.Tests/Unit/MetricsTests.cs +++ b/tests/Unsga3.Tests/Unit/MetricsTests.cs @@ -53,6 +53,70 @@ public void HV2D_two_points() Assert.InRange(hv, 0.54, 0.56); } + [Fact] + public void IgdPlus_and_gd_match_the_hand_case() + { + // A = {(0, 1)} against Z = {(0, 0), (1, 0)}. + // IGD+ modified distance is 1 for both reference points. GD nearest distance is 1. + var obtained = new[] { new[] { 0.0, 1.0 } }; + var reference = new[] { new[] { 0.0, 0.0 }, new[] { 1.0, 0.0 } }; + double igdPlus = PerformanceIndicators.InvertedGenerationalDistancePlus(obtained, reference); + double gd = PerformanceIndicators.GenerationalDistance(obtained, reference); + Assert.Equal(1.0, igdPlus, 9); + Assert.Equal(1.0, gd, 9); + } + + [Fact] + public void Indicators_reject_mismatched_objective_lengths() + { + var wide = new[] { new[] { 0.0, 0.0 } }; + var shortVector = new[] { new[] { 0.0 } }; + Assert.Throws(() => + PerformanceIndicators.InvertedGenerationalDistance(wide, shortVector)); + Assert.Throws(() => + PerformanceIndicators.GenerationalDistance(shortVector, wide)); + Assert.Throws(() => + PerformanceIndicators.InvertedGenerationalDistancePlus(wide, shortVector)); + } + + [Fact] + public void Zdt_samplers_reject_a_single_point() + { + Assert.Throws(() => ParetoFronts.Zdt1(1)); + Assert.Throws(() => ParetoFronts.Zdt2(1)); + Assert.Throws(() => ParetoFronts.Zdt4(1)); + Assert.Throws(() => ParetoFronts.Zdt6(1)); + Assert.Throws(() => ParetoFronts.Zdt3(1)); + } + + [Fact] + public void Zdt3_second_segment_starts_at_the_library_literal() + { + // pymoo 0.6.2 uses 0.182228780. This library uses 0.1822287280. + var front = ParetoFronts.Zdt3(pointsPerSegment: 2); + Assert.Equal(0.1822287280, front[2][0], 12); + } + + [Fact] + public void Zdt6_floor_is_below_the_sampled_minimum() + { + var front = ParetoFronts.Zdt6(2); + Assert.Equal(0.280775, front[0][0], 12); + + double min = double.PositiveInfinity; + const int steps = 200_000; + for (int i = 0; i <= steps; i++) + { + double x = i / (double)steps; + double s = Math.Sin(6.0 * Math.PI * x); + double f1 = 1.0 - Math.Exp(-4.0 * x) * Math.Pow(s * s, 3); + if (f1 < min) min = f1; + } + + double gap = min - front[0][0]; + Assert.InRange(gap, 1e-7, 1e-6); + } + [Fact] public void ParetoFronts_Zdt1_on_curve() { @@ -60,6 +124,30 @@ public void ParetoFronts_Zdt1_on_curve() Assert.InRange(p[1], 1.0 - Math.Sqrt(p[0]) - 1e-9, 1.0 - Math.Sqrt(p[0]) + 1e-9); } + [Fact] + public void OnePerDirection_keeps_the_closer_point_on_a_shared_ray() + { + var directions = new[] { new[] { 1.0, 0.0 }, new[] { 0.0, 1.0 } }; + var onAxis = new[] { 1.0, 0.0 }; + var nearby = new[] { 0.9, 0.1 }; + var other = new[] { 0.0, 1.0 }; + var kept = ReferenceDirectionThinning.OnePerDirection( + new[] { nearby, onAxis, other }, + directions); + + Assert.Equal(2, kept.Length); + Assert.Contains(kept, p => p[0] == 1.0 && p[1] == 0.0); + Assert.Contains(kept, p => p[0] == 0.0 && p[1] == 1.0); + } + + [Fact] + public void OnePerDirection_rejects_a_short_objective_vector() + { + var directions = new[] { new[] { 1.0, 0.0 } }; + Assert.Throws(() => + ReferenceDirectionThinning.OnePerDirection(new[] { new[] { 1.0 } }, directions)); + } + [Fact] public void Dtlz2_front_on_unit_sphere() { diff --git a/tests/Unsga3.Tests/Unit/NormalizationTests.cs b/tests/Unsga3.Tests/Unit/NormalizationTests.cs index 4995505..b0d730c 100644 --- a/tests/Unsga3.Tests/Unit/NormalizationTests.cs +++ b/tests/Unsga3.Tests/Unit/NormalizationTests.cs @@ -95,6 +95,72 @@ public void Ideal_point_is_persistent_across_calls() Assert.Equal(0.1, norm.IdealPoint[1], 9); } + [Fact] + public void Collapsed_span_sets_nadir_to_ideal_plus_one() + { + // {2, 2+1e-8}. Span stays below 1e-6 after the worst-of-population fallback, + // so nadir becomes ideal + 1 = 3. pymoo 0.6.2 would keep the 1e-8 span. + var norm = new Normalization(1); + var pop = new List { Make1(2.0), Make1(2.0 + 1e-8) }; + var normalized = norm.Normalize(pop); + + Assert.Equal(2.0, norm.IdealPoint[0], 12); + Assert.Equal(3.0, norm.NadirPoint[0], 12); + Assert.Equal(0.0, normalized[0][0], 9); + Assert.InRange(normalized[1][0], 1e-9, 1e-7); + } + + [Fact] + public void Infeasible_origin_sets_ideal_from_the_whole_pool() + { + // Locks the current rule. Feasible (1,1) vs infeasible (0,0). + // Ideal and worst are taken from the whole pool. Constraint-domination + // puts only the feasible point on the first front, so the ND extreme + // search does not see the origin; the ideal still does. + var feasible = new Individual(1, 2, 1); + feasible.Objectives[0] = 1; + feasible.Objectives[1] = 1; + feasible.Constraints[0] = 0; + feasible.RefreshConstraintViolation(); + + var infeasible = new Individual(1, 2, 1); + infeasible.Objectives[0] = 0; + infeasible.Objectives[1] = 0; + infeasible.Constraints[0] = 1; + infeasible.RefreshConstraintViolation(); + + Assert.True(feasible.IsFeasible); + Assert.False(infeasible.IsFeasible); + + var pop = new List { feasible, infeasible }; + var fronts = NonDominatedSort.Sort(pop); + Assert.Equal(new[] { 0 }, fronts[0]); + + var norm = new Normalization(2); + var normalized = norm.Normalize(pop, fronts[0]); + + Assert.Equal(0.0, norm.IdealPoint[0], 12); + Assert.Equal(0.0, norm.IdealPoint[1], 12); + Assert.Equal(1.0, norm.NadirPoint[0], 12); + Assert.Equal(1.0, norm.NadirPoint[1], 12); + Assert.Equal(1.0, normalized[0][0], 9); + Assert.Equal(1.0, normalized[0][1], 9); + Assert.Equal(0.0, normalized[1][0], 9); + Assert.Equal(0.0, normalized[1][1], 9); + + // Unfiltered ASF scores the infeasible origin ahead of (1,1). + double[] ideal = { 0.0, 0.0 }; + Assert.True(Normalization.Asf(infeasible.Objectives, 0, ideal) + < Normalization.Asf(feasible.Objectives, 0, ideal)); + } + + private static Individual Make1(double a) + { + var ind = new Individual(1, 1); + ind.Objectives[0] = a; + return ind; + } + private static Individual Make(double a, double b, double c) { var ind = new Individual(1, 3); diff --git a/tests/Unsga3.Tests/Unit/TournamentSelectionTests.cs b/tests/Unsga3.Tests/Unit/TournamentSelectionTests.cs index c8e540c..6ac62f0 100644 --- a/tests/Unsga3.Tests/Unit/TournamentSelectionTests.cs +++ b/tests/Unsga3.Tests/Unit/TournamentSelectionTests.cs @@ -1,4 +1,5 @@ using Unsga3.Algorithm; +using Unsga3.Operators.Crossover; using Unsga3.Operators.Selection; using Unsga3.Utilities; @@ -61,4 +62,93 @@ public void Pymoo_different_niche_is_random_not_rank() } Assert.InRange(betterWins, 5, 35); // not deterministic rank dominance } + + [Fact] + public void Pymoo_same_niche_equal_rank_prefers_shorter_distance() + { + var closer = new Individual(1, 2) + { + Rank = 0, + AssociatedReference = 2, + PerpendicularDistance = 0.1, + NicheCount = 9, + }; + var farther = new Individual(1, 2) + { + Rank = 0, + AssociatedReference = 2, + PerpendicularDistance = 0.4, + NicheCount = 1, + }; + Assert.Same(closer, TournamentSelection.Winner( + closer, farther, new RandomProvider(0), TournamentMode.PymooCompatible)); + } + + [Fact] + public void Pymoo_same_niche_distance_tie_is_a_coin_flip() + { + var a = new Individual(1, 2) { Rank = 0, AssociatedReference = 4, PerpendicularDistance = 0.2 }; + var b = new Individual(1, 2) { Rank = 0, AssociatedReference = 4, PerpendicularDistance = 0.2 }; + int aWins = 0; + for (int seed = 0; seed < 40; seed++) + { + var w = TournamentSelection.Winner(a, b, new RandomProvider(seed), TournamentMode.PymooCompatible); + if (ReferenceEquals(w, a)) aWins++; + } + Assert.InRange(aWins, 5, 35); + } + + [Fact] + public void RankNiche_prefers_lower_niche_count_across_different_niches() + { + var sparse = new Individual(1, 2) + { + Rank = 0, + NicheCount = 1, + AssociatedReference = 0, + PerpendicularDistance = 0.5, + }; + var crowded = new Individual(1, 2) + { + Rank = 0, + NicheCount = 6, + AssociatedReference = 1, + PerpendicularDistance = 0.01, + }; + Assert.Same(sparse, TournamentSelection.Winner( + sparse, crowded, new RandomProvider(1), TournamentMode.RankNicheDistance)); + } + + [Fact] + public void Pymoo_ignores_niche_count_when_niches_differ() + { + var sparse = new Individual(1, 2) + { + Rank = 1, + NicheCount = 1, + AssociatedReference = 0, + PerpendicularDistance = 0.5, + }; + var crowded = new Individual(1, 2) + { + Rank = 0, + NicheCount = 6, + AssociatedReference = 1, + PerpendicularDistance = 0.01, + }; + int sparseWins = 0; + for (int seed = 0; seed < 40; seed++) + { + var w = TournamentSelection.Winner(sparse, crowded, new RandomProvider(seed), TournamentMode.PymooCompatible); + if (ReferenceEquals(w, sparse)) sparseWins++; + } + Assert.InRange(sparseWins, 5, 35); + } + + [Fact] + public void Default_sbx_probability_is_one() + { + // pymoo NSGA3/UNSGA3 uses SBX(prob=1.0). Seada & Deb section 4 uses pc = 0.9. + Assert.Equal(1.0, new SimulatedBinaryCrossover().Probability); + } } diff --git a/tests/Unsga3.Tests/Unit/WithDasDennisTests.cs b/tests/Unsga3.Tests/Unit/WithDasDennisTests.cs new file mode 100644 index 0000000..a28e2c8 --- /dev/null +++ b/tests/Unsga3.Tests/Unit/WithDasDennisTests.cs @@ -0,0 +1,24 @@ +using Unsga3.Algorithm; +using Unsga3.Problems; + +namespace Unsga3.Tests.Unit; + +public class WithDasDennisTests +{ + [Fact] + public void Single_objective_default_population_throws() + { + // Das–Dennis for M=1 is one direction. N defaults to |H| and must be at least 2. + var ex = Assert.Throws(() => Unsga3Algorithm.WithDasDennis(1, 1)); + Assert.Equal("populationSize", ex.ParamName); + } + + [Fact] + public void Single_objective_runs_when_caller_chooses_population() + { + var algo = Unsga3Algorithm.WithDasDennis(1, 1, populationSize: 8, seed: 1); + var result = algo.Run(new SphereProblem(nVariables: 2), maxGenerations: 2); + Assert.Equal(8, result.FinalPopulation.Count); + Assert.Equal(2, result.GenerationsExecuted); + } +} diff --git a/tools/OracleCompare/Program.cs b/tools/OracleCompare/Program.cs index 738e0c9..6e1a627 100644 --- a/tools/OracleCompare/Program.cs +++ b/tools/OracleCompare/Program.cs @@ -8,6 +8,7 @@ using Unsga3.Utilities; // Fixed-protocol C# side of the pymoo oracle (see tools/oracle/run_pymoo_oracle.py). +// IGD is scored on the full feasible non-dominated front. pymoo's script scores res.F. // // dotnet run --project tools/OracleCompare -- --problem zdt1 --partitions 12 --pop 52 --gens 100 --seed 1 --pymoo-mode // # ZDT2 quality protocol (matches unsga3-bend A/B): gens=250 + --pymoo-mode. @@ -76,7 +77,7 @@ var mode = pymooMode ? TournamentMode.PymooCompatible : TournamentMode.RankNicheDistance; Console.WriteLine( - $"Unsga3 | problem={problemName} M={m} refs={dirs.Length} pop={popSize} gens={gens} seed={seed} tournament={mode}"); + $"Unsga3 | problem={problemName} M={m} n_var={problem.NumberOfVariables} refs={dirs.Length} pop={popSize} gens={gens} seed={seed} tournament={mode}"); var algo = new Unsga3Algorithm(dirs, popSize, seed: seed, tournamentMode: mode); var result = algo.Run(problem, gens); @@ -93,6 +94,7 @@ ? PerformanceIndicators.Hypervolume2D(obtained, new[] { 1.1, 1.1 }) : null; +Console.WriteLine($"front=non_dominated n={obtained.Length}"); Console.WriteLine($"front_size={obtained.Length}"); Console.WriteLine($"IGD={igd.ToString("G6", CultureInfo.InvariantCulture)}"); if (hv is double h) @@ -116,12 +118,14 @@ ["tournament"] = mode.ToString(), ["problem"] = problemName, ["n_obj"] = m, + ["n_var"] = problem.NumberOfVariables, ["partitions"] = partitions, ["n_ref_dirs"] = dirs.Length, ["pop_size"] = popSize, ["n_gen"] = gens, ["seed"] = seed, ["n_solutions"] = obtained.Length, + ["front_definition"] = "non_dominated_feasible", ["igd"] = igd, ["hv2"] = hv, ["F_csv"] = Path.GetFileName(fPath), diff --git a/tools/oracle/analyze_dtlz2_gap.py b/tools/oracle/analyze_dtlz2_gap.py index 15d5777..a83acf0 100644 --- a/tools/oracle/analyze_dtlz2_gap.py +++ b/tools/oracle/analyze_dtlz2_gap.py @@ -64,7 +64,9 @@ def main() -> None: from pymoo.util.ref_dirs import get_reference_directions ref = get_reference_directions("das-dennis", 3, n_partitions=12) - pf = get_problem("dtlz2", n_obj=3).pareto_front(ref) + # Spherical PF does not depend on k. n_var=12 matches Dtlz2Problem(k=10); + # pymoo's default n_var=10 (k=8) is a different search problem, not a different PF. + pf = get_problem("dtlz2", n_obj=3, n_var=12).pareto_front(ref) print("PF size", len(pf), "refs", len(ref)) for F, lab in [(cs, "cs-pymoo"), (csd, "cs-def"), (py, "pymoo")]: diff --git a/tools/oracle/run_multiseed_wilcoxon.py b/tools/oracle/run_multiseed_wilcoxon.py index 34e85ed..52ca71e 100644 --- a/tools/oracle/run_multiseed_wilcoxon.py +++ b/tools/oracle/run_multiseed_wilcoxon.py @@ -31,6 +31,14 @@ ORACLE_COMPARE = ROOT / "tools" / "OracleCompare" +# C# Dtlz2Problem(k=10) ⇒ n_var = M + k - 1. pymoo's default for M=3 is n_var=10 (k=8). +DTLZ2_K = 10 + + +def dtlz2_n_var(n_obj: int, k: int = DTLZ2_K) -> int: + return n_obj + k - 1 + + @dataclass class Protocol: name: str @@ -39,6 +47,7 @@ class Protocol: gens: int n_obj: int csharp_pymoo_mode: bool # True → TournamentMode.PymooCompatible + n_var: int | None = None # DTLZ2 sets 12; None keeps the pymoo problem default PROTOCOLS: dict[str, Protocol] = { @@ -47,7 +56,17 @@ class Protocol: # Optional: csharp_pymoo_mode=False (RankNicheDistance) + gens=100 is the # unpublished Wilcoxon ZDT2 mating snapshot — do not treat it as the default. "zdt2": Protocol("zdt2", partitions=12, pop=52, gens=250, n_obj=2, csharp_pymoo_mode=True), - "dtlz2": Protocol("dtlz2", partitions=12, pop=92, gens=150, n_obj=3, csharp_pymoo_mode=True), + # n_var=12 matches C#. The checked-in WILCOXON-RESULTS.md pymoo column is the + # older n_var=10 run and must not be regenerated until those seeds are re-run. + "dtlz2": Protocol( + "dtlz2", + partitions=12, + pop=92, + gens=150, + n_obj=3, + csharp_pymoo_mode=True, + n_var=dtlz2_n_var(3), + ), } @@ -178,7 +197,8 @@ def run_pymoo(proto: Protocol, seed: int) -> float: import numpy as np if proto.name == "dtlz2": - problem = get_problem("dtlz2", n_obj=proto.n_obj) + n_var = proto.n_var if proto.n_var is not None else dtlz2_n_var(proto.n_obj) + problem = get_problem("dtlz2", n_obj=proto.n_obj, n_var=n_var) ref_dirs = get_reference_directions( "das-dennis", proto.n_obj, n_partitions=proto.partitions ) @@ -206,6 +226,7 @@ def run_pymoo(proto: Protocol, seed: int) -> float: "source": "pymoo", "algorithm": "UNSGA3", "problem": proto.name, + "n_var": int(problem.n_var), "pop_size": proto.pop, "n_gen": proto.gens, "seed": seed, @@ -332,8 +353,8 @@ def markdown_report(results: dict) -> str: "## Notes", "", "- Not bit-identical: different RNG implementations and minor operator ordering.", - "- Practical equivalence: median IGD within ~1–2× and non-significant MWU is a strong claim; " - "significant differences with small effect size (ratio ≈ 1) are still acceptable for a v0.x port.", + "- These tests are not a 1–2% equivalence claim. Report the Mann–Whitney result as computed. " + "DTLZ2 pymoo runs in this script use n_var=12; a table generated before that change is the n_var=10 column.", "- Reproduce: `python tools/oracle/run_multiseed_wilcoxon.py`", "", ] @@ -403,6 +424,7 @@ def main() -> int: "partitions": proto.partitions, "pop": proto.pop, "gens": proto.gens, + "n_var": proto.n_var, "csharp_pymoo_mode": proto.csharp_pymoo_mode, }, "csharp": summarize(cs_igds), diff --git a/tools/oracle/run_pymoo_oracle.py b/tools/oracle/run_pymoo_oracle.py index 955e871..c0a667a 100644 --- a/tools/oracle/run_pymoo_oracle.py +++ b/tools/oracle/run_pymoo_oracle.py @@ -6,10 +6,21 @@ SBX η=30, PM η=20 (pymoo defaults for NSGA3/UNSGA3), Das-Dennis refs, seed=1, export final F + IGD. +The exported front is pymoo res.F (the survival niche set, about one point per +filled reference direction). It is not the final population's full non-dominated +front. C# OracleCompare scores that full front. Compare those files only after +putting both sides on the same front definition and the same reference set. + +DTLZ2 decision dimension: C# Dtlz2Problem(k=10) uses n_var = M + k - 1 = 12. +pymoo get_problem("dtlz2", n_obj=3) defaults to n_var=10 (k=8). This script +passes n_var=12 unless --n-var is set. --n-var 10 reproduces the historical +mismatched column only; it is not the apples-to-apples protocol. + Usage: python run_pymoo_oracle.py python run_pymoo_oracle.py --problem zdt1 --partitions 12 --pop 52 --gens 100 --seed 1 python run_pymoo_oracle.py --problem zdt2 --partitions 12 --pop 52 --seed 1 + python run_pymoo_oracle.py --problem dtlz2 --partitions 12 --pop 92 --gens 150 --seed 1 # omitted --gens on zdt2 is 250 (quality protocol; matches unsga3-bend A/B). # --gens 100 is an early-stress snapshot, not the quality bar. """ @@ -23,6 +34,15 @@ import numpy as np +# Deb et al. suggest k=10 for DTLZ2. C# Dtlz2Problem defaults to that k. +DTLZ2_K = 10 + + +def dtlz2_n_var(n_obj: int, k: int = DTLZ2_K) -> int: + """n = M + k - 1. For M=3, k=10 this is 12, not pymoo's default 10.""" + return n_obj + k - 1 + + def default_gens(problem: str) -> int: """Quality-protocol generations when --gens is omitted. @@ -47,6 +67,16 @@ def main() -> int: help="generations (default: zdt2=250, else 100; explicit value always wins)", ) p.add_argument("--seed", type=int, default=1) + p.add_argument( + "--n-var", + type=int, + default=None, + help=( + "Decision variables. DTLZ2 default is M+k-1 with k=10 (n_var=12), " + "matching Dtlz2Problem(k:10). pymoo's own default is 10 (k=8); " + "pass --n-var 10 only to reproduce that historical mismatched column." + ), + ) p.add_argument("--out-dir", type=Path, default=Path(__file__).resolve().parent / "out") args = p.parse_args() gens = args.gens if args.gens is not None else default_gens(args.problem) @@ -62,21 +92,26 @@ def main() -> int: print(e, file=sys.stderr) return 2 + n_var: int | None if args.problem == "dtlz2": n_obj = 3 - problem = get_problem("dtlz2", n_obj=n_obj) + n_var = args.n_var if args.n_var is not None else dtlz2_n_var(n_obj) + problem = get_problem("dtlz2", n_obj=n_obj, n_var=n_var) ref_dirs = get_reference_directions("das-dennis", n_obj, n_partitions=args.partitions) pf = problem.pareto_front(ref_dirs) else: n_obj = 2 problem = get_problem(args.problem) + n_var = int(problem.n_var) if args.n_var is None else args.n_var + if args.n_var is not None: + problem = get_problem(args.problem, n_var=n_var) ref_dirs = get_reference_directions("das-dennis", n_obj, n_partitions=args.partitions) pf = problem.pareto_front() pop = args.pop if args.pop is not None else len(ref_dirs) algo = UNSGA3(ref_dirs, pop_size=pop) - print(f"pymoo UNSGA3 | problem={args.problem} M={n_obj} refs={len(ref_dirs)} " + print(f"pymoo UNSGA3 | problem={args.problem} M={n_obj} n_var={n_var} refs={len(ref_dirs)} " f"pop={pop} gens={gens} seed={args.seed}") res = minimize( @@ -101,17 +136,21 @@ def main() -> int: "algorithm": "UNSGA3", "problem": args.problem, "n_obj": n_obj, + "n_var": n_var, + "k": (n_var - n_obj + 1) if args.problem == "dtlz2" else None, "partitions": args.partitions, "n_ref_dirs": int(len(ref_dirs)), "pop_size": pop, "n_gen": gens, "seed": args.seed, "n_solutions": int(F.shape[0]), + "front_definition": "res.F", "igd": igd, "F_csv": str(f_path.name), } meta_path.write_text(json.dumps(meta, indent=2), encoding="utf-8") + print(f"front=res.F n={int(F.shape[0])} (survival optimum, not the full population ND front)") print(f"IGD={igd:.6g}") print(f"wrote {f_path}") print(f"wrote {meta_path}")