diff --git a/CHANGELOG.md b/CHANGELOG.md index 03aa133..ebae0be 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,17 +10,21 @@ and this project adheres to [Semantic Versioning](https://semver.org/). ### Added - Optional **TypeSafe / Jev System One** pass (`tools/typesafe-pareto`) for Score + Choice judgments over a small Pareto candidate sample. Additive semantic layer — does not replace NSGA-III / U-NSGA-III objectives. Live calls only when `TYPESAFE_API_KEY` is set; CI uses mocks. Metrics append to `metrics/typesafe-runs.jsonl`. +- `ReferenceDirections.TwoLayerDasDennis`: outer Das–Dennis layer plus an inside layer scaled toward the centroid (Deb & Jain, NSGA-III Part I; default scale 0.5). Duplicates are removed. Single-layer `DasDennis` / `WithDasDennis` stays the default. M≤3 oracles do not need the second layer. +- Constrained problems `OsyProblem`, `TnkProblem`, and `C1Dtlz1Problem`, with formulation and survival self-tests. No IGD table. - Opt-in `MatingPoolMode.TwoShuffledPasses`: two shuffled passes of consecutive tournaments (Seada & Deb Algorithm 1 and pymoo 0.6.2 `TournamentSelection`, pressure 2). When N parents are drawn from N members, each member is a contestant twice. The default pool stays N independent tournaments with replacement, which is the published ZDT1 Wilcoxon protocol. - Opt-in `TournamentMode.Algorithm2`: a same-niche perpendicular-distance tie keeps the second parent. `PymooCompatible` still coin-flips. `RankNicheDistance` is unchanged. -### Fixed - -- Odd population sizes no longer skip `parents[N-1]` in SBX. The leftover parent is paired with parent 0. Even-N pair order is unchanged (ZDT1 pop 52, DTLZ2 pop 92). - ### Changed +- Environmental selection follows pymoo `Survival.filter_infeasible` (Jain & Deb, NSGA-III Part II). Feasible members are niched. If fewer than N are feasible, the rest are filled by ascending constraint violation and are not assigned a niche. When any member is feasible, ideal, worst, and ASF extremes come from feasible objectives. An all-infeasible generation does not move the hyperplane. This replaces the locked whole-pool ideal fixture. - Docs and XML comments describe `initialPopulation` and hybrid loops in generic terms (domain-adapter warm-start / grid-seed). No product-repo names. +### Fixed + +- `CompareConstraintDominated` no longer falls through to Pareto when two infeasible individuals share a constraint violation. Deb's constraint-domination leaves them mutually non-dominated. Tournament comparators are unchanged. +- Odd population sizes no longer skip `parents[N-1]` in SBX. The leftover parent is paired with parent 0. Even-N pair order is unchanged (ZDT1 pop 52, DTLZ2 pop 92). + ## [0.1.5] — 2026-09-22 Forensic-audit wave (G1–G10): document known gaps against pymoo and Seada & Deb, and lock them with tests. The published 15-seed table is not rewritten. **Default tournament stays `RankNicheDistance`.** diff --git a/README.md b/README.md index 440ca57..3e95035 100644 --- a/README.md +++ b/README.md @@ -73,13 +73,14 @@ var algo = new Unsga3Algorithm(dirs, populationSize: 92, seed: 1, | `Unsga3Algorithm` | Main entry — `Run(problem, gens)` | | `Individual` | Variables / Objectives / Constraints | | `OptimizationResult` | Final population + non-dominated set | -| `ReferenceDirections.DasDennis` | Structured reference points | +| `ReferenceDirections.DasDennis` | Single-layer reference points (default) | +| `ReferenceDirections.TwoLayerDasDennis` | Outer layer plus an inside layer for larger M | | `SimulatedBinaryCrossover` / `PolynomialMutation` | Variation operators | | `PerformanceIndicators` | IGD, GD, 2-D hypervolume | | `TournamentMode` | Default rank→niche, `PymooCompatible` (distance-tie coin flip), or opt-in `Algorithm2` (keep the second parent) | | `MatingPoolMode` | Default independent tournaments; opt-in `TwoShuffledPasses` (paper / pymoo 0.6.2) | -Built-in problems: ZDT1–4/6, DTLZ1–4/7, Sphere, Ackley, Rosenbrock. +Built-in problems: ZDT1–4/6, DTLZ1–4/7, Sphere, Ackley, Rosenbrock, and constrained OSY, TNK, and C1-DTLZ1 (self-tests only; no IGD table yet). ## Build & test @@ -154,6 +155,7 @@ Algorithm papers (please cite these when publishing results): - Seada, H. & Deb, K. (2016). *A Unified Evolutionary Optimization Procedure for Single, Multiple, and Many Objectives.* IEEE Trans. Evol. Comput. - Deb, K. & Jain, H. (2014). *An Evolutionary Many-Objective Optimization Algorithm Using Reference-Point-Based Nondominated Sorting Approach (NSGA-III), Part I.* IEEE Trans. Evol. Comput. +- Jain, H. & Deb, K. (2014). *NSGA-III, Part II: Handling Constraints and Extending to an Adaptive Approach.* IEEE Trans. Evol. Comput. - Das, I. & Dennis, J. E. (1998). *Normal-Boundary Intersection.* SIAM J. Optim. ## License diff --git a/docs/EQUIVALENCE.md b/docs/EQUIVALENCE.md index 3d8febd..d36cdae 100644 --- a/docs/EQUIVALENCE.md +++ b/docs/EQUIVALENCE.md @@ -11,6 +11,7 @@ See also **[RESEARCH-STANDARDS.md](RESEARCH-STANDARDS.md)** for the literature + - COIN Report 2014022; Seada & Deb, IEEE TEVC 2016 - [pymoo `UNSGA3`](https://pymoo.org/algorithms/moo/unsga3.html) - pymoo `nsga3.py` — `HyperplaneNormalization`, `associate_to_niches`, `niching` +- pymoo `Survival.filter_infeasible` — feasible niching, then ascending-CV fill (Jain & Deb, NSGA-III Part II) - Indicators: [pymoo performance indicators](https://pymoo.org/misc/indicators.html) (GD, IGD, IGD+, HV) ## Protocol @@ -44,12 +45,13 @@ ZDT2 **gens=100** is an early-stress snapshot (collapse on Bend, C#, and pymoo), | Bi | ZDT1–4, ZDT6 | 2 | yes | | Many | DTLZ1–4, DTLZ7 | 3+ | yes | | Hard | WFG1, WFG2, WFG9 | 3–5 | planned | -| Constrained | OSY, TNK, C1-DTLZ1 | 2–3 | planned | +| Constrained | OSY, TNK, C1-DTLZ1 | 2–3 | yes (self-tests; no IGD table) | ## Unit checks - Reference association & niche counts -- Non-dominated ranks / constraint domination +- Non-dominated ranks / constraint domination (equal CV is mutual non-domination) +- Feasible niching, CV fill, and OSY / TNK / C1-DTLZ1 (`ConstrainedSurvivalTests`, `ConstrainedProblemTests`). No IGD table. - **Normalization intercepts / ASF axis extremes** (`NormalizationTests`) - Tournament pressure - `populationSize: null` ⇒ `|refs|` @@ -71,14 +73,17 @@ ZDT2 **gens=100** is an early-stress snapshot (collapse on Bend, C#, and pymoo), | Mating pool (`TwoShuffledPasses`, opt-in) | two shuffled passes of consecutive pairs. When N parents are drawn from N members, each index is a contestant twice | pymoo 0.6.2 `TournamentSelection` (pressure 2) **is** this pool: `n_perms = ceil(2N / N)` permutations, then consecutive pairs | | 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 | +| Reference layers | `DasDennis` / `WithDasDennis` stay single-layer (the M≤3 oracles). `TwoLayerDasDennis` adds Part I's inside layer, default scale 0.5, duplicates removed | multi-layer factory; NSGA-III consumes whatever `ref_dirs` it is given | | 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 | +| Constrained survival | Niche feasible members only. If fewer than N are feasible, fill by ascending CV (original index breaks ties) with no niche assignment. All-infeasible generations truncate by CV | `Survival.filter_infeasible` | +| Feasible hyperplane | Ideal, worst, and ASF use feasible objectives when any feasible member exists. An all-infeasible generation does not move the hyperplane | `HyperplaneNormalization.update` sees the feasible subset passed by `filter_infeasible` | +| Equal CV in the sort | Two infeasible individuals with the same violation are mutually non-dominated (Deb). Objectives are not compared | constraint-domination; CV ordering is separate from Pareto | +| Tournament equal CV | unchanged. `RankNicheDistance` still falls through to rank, niche, and distance. `PymooCompatible` still coin-flips | pymoo coin-flips | ## Mating pool diff --git a/docs/NOTICE.md b/docs/NOTICE.md index 24eabf2..ba290f4 100644 --- a/docs/NOTICE.md +++ b/docs/NOTICE.md @@ -12,7 +12,8 @@ It is **not** affiliated with, endorsed by, or derived from the source code of p | Work | Role in Unsga3 | |------|----------------| | Seada, H. & Deb, K. (2016). *A Unified Evolutionary Optimization Procedure for Single, Multiple, and Many Objectives.* IEEE TEVC. | U-NSGA-III algorithm | -| Deb, K. & Jain, H. (2014). *NSGA-III.* IEEE TEVC. | Reference-point survival, hyperplane normalization | +| Deb, K. & Jain, H. (2014). *An Evolutionary Many-Objective Optimization Algorithm Using Reference-Point-Based Nondominated Sorting Approach, Part I: Solving Problems With Box Constraints.* IEEE TEVC. | Reference directions (including the inside layer), hyperplane normalization, niching | +| Jain, H. & Deb, K. (2014). *An Evolutionary Many-Objective Optimization Algorithm Using Reference-Point-Based Nondominated Sorting Approach, Part II: Handling Constraints and Extending to an Adaptive Approach.* IEEE TEVC. | Feasible-only niching and constraint-violation fill | | Das, I. & Dennis, J. E. (1998). *Normal-Boundary Intersection.* SIAM J. Optim. | Das–Dennis reference directions | | Deb, K. et al. ZDT / DTLZ suites | Standard test problems | | Deb & Agrawal — SBX; Deb et al. — polynomial mutation | Variation operators | diff --git a/docs/ORACLE-RESULTS.md b/docs/ORACLE-RESULTS.md index 8351c11..da7ec2f 100644 --- a/docs/ORACLE-RESULTS.md +++ b/docs/ORACLE-RESULTS.md @@ -131,7 +131,7 @@ ZDT2 quality A/B is gens=250 + `PymooCompatible` (not the loose smoke bar). ZDT1 | IGD mean-distance | **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. | +| Infeasible points | **feasible hyperplane**. Ideal, worst, and ASF use feasible objectives when any feasible member exists. An all-infeasible generation does not move the hyperplane. Survival niches the feasible set and fills a shortfall by ascending CV. The old whole-pool fixture (feasible (1, 1) vs infeasible (0, 0) → ideal (0, 0)) is replaced. OSY, TNK, and C1-DTLZ1 have self-tests. No constrained IGD table. | | 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** | diff --git a/docs/RESEARCH-STANDARDS.md b/docs/RESEARCH-STANDARDS.md index a1445d3..3d70092 100644 --- a/docs/RESEARCH-STANDARDS.md +++ b/docs/RESEARCH-STANDARDS.md @@ -19,7 +19,7 @@ Sources consulted (2026-08): | Bi-objective | ZDT1–4, ZDT6 | Convex / non-convex / disconnected / multimodal / biased density | | Many-objective | DTLZ1–4, DTLZ7 | Linear simplex, sphere, multimodal sphere, bias, disconnected | | Hard (later) | WFG1–9 | Non-separable, scaled, deceptive | -| Constrained (later) | OSY, TNK, C1-DTLZ1 | Constraint-domination path | +| Constrained | OSY, TNK, C1-DTLZ1 | In the library, with survival self-tests. No IGD table until measured runs exist | **Seada selection paper** stresses ZDT4 (many local fronts) and DTLZ1 for M=3,5. @@ -35,7 +35,7 @@ Default dimensions (Deb / pymoo convention): |------|----------------| | 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 the unit simplex, **single layer**. Two-layer directions for larger M are absent | +| Reference set | **Das–Dennis** on the unit simplex. **Single layer is the default** (`DasDennis`, `WithDasDennis`, M≤3 oracles). `TwoLayerDasDennis` adds Part I's inside layer (default scale 1/2) for larger M | | 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) | diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index 7dc1952..d38075d 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -18,7 +18,8 @@ Living plan for Unsga3. Issues track concrete work; this page is the narrative. - [ ] `net8.0` (+ `net10.0`) multi-target for broader NuGet consumers - [ ] nuget.org publish (in addition to GitHub Packages) - [ ] GD+ indicator (IGD+ has a hand-case test; GD and IGD were already implemented) -- [ ] Constrained demos (OSY / TNK) with self-tests +- [x] Constrained demos (OSY, TNK, C1-DTLZ1) with formulation and survival self-tests +- [ ] Constrained IGD / Wilcoxon vs pymoo (no table until measured runs exist) - [ ] API docs site (DocFX or similar) ## Medium term (0.3+) diff --git a/src/Unsga3/Algorithm/Unsga3Algorithm.cs b/src/Unsga3/Algorithm/Unsga3Algorithm.cs index db1fe8f..a2fa5d2 100644 --- a/src/Unsga3/Algorithm/Unsga3Algorithm.cs +++ b/src/Unsga3/Algorithm/Unsga3Algorithm.cs @@ -85,7 +85,9 @@ public Unsga3Algorithm( } /// - /// Convenience: build Das–Dennis directions then construct the algorithm. + /// Convenience: build single-layer Das–Dennis directions then construct the algorithm. + /// Two-layer directions are , passed + /// to the constructor; M≤3 oracles stay on this single-layer factory. /// 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 diff --git a/src/Unsga3/Core/NonDominatedSort.cs b/src/Unsga3/Core/NonDominatedSort.cs index 48e7966..eeeda2f 100644 --- a/src/Unsga3/Core/NonDominatedSort.cs +++ b/src/Unsga3/Core/NonDominatedSort.cs @@ -85,7 +85,10 @@ public static int CompareConstraintDominated(Individual a, Individual b) { if (a.ConstraintViolation < b.ConstraintViolation) return -1; if (a.ConstraintViolation > b.ConstraintViolation) return 1; - // equal violation → fall through to objective comparison + // Deb constraint-domination: equal violation ⇒ neither dominates. + // Objectives are not consulted. Tournament comparators are separate + // and still break an equal-CV tie their own way. + return 0; } return ComparePareto(a.Objectives, b.Objectives); diff --git a/src/Unsga3/Core/Normalization.cs b/src/Unsga3/Core/Normalization.cs index 222eeb6..294d45a 100644 --- a/src/Unsga3/Core/Normalization.cs +++ b/src/Unsga3/Core/Normalization.cs @@ -15,12 +15,13 @@ namespace Unsga3.Core; /// 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. +/// When any member is feasible, ideal, worst, and ASF extremes are taken from the +/// feasible subset only (Jain & Deb Part II; pymoo Survival.filter_infeasible +/// passes that subset into HyperplaneNormalization). An all-infeasible +/// population does not move the persistent hyperplane. See +/// Feasible_members_set_the_hyperplane. +/// Extreme-point ASF uses the non-dominated index set when the caller supplies one, +/// intersected with the feasible subset, and the feasible members when it does not. /// public sealed class Normalization { @@ -60,9 +61,12 @@ public Normalization(int nObjectives) /// /// Update ideal / nadir from and return normalized - /// objective vectors (same order). When is - /// provided, extreme points are sought only among that subset (pymoo / NSGA-III); - /// otherwise the whole population is used. + /// objective vectors (same order). When any member is feasible, only feasible + /// objectives move the hyperplane. An all-infeasible population is projected with + /// the current hyperplane and does not update it. When + /// is provided, extreme points are sought + /// only among that subset (pymoo / NSGA-III), still restricted to feasible members; + /// otherwise the feasible members are used. /// public double[][] Normalize( IReadOnlyList population, @@ -72,10 +76,31 @@ public double[][] Normalize( int n = population.Count; if (n == 0) return Array.Empty(); - // Persistent ideal / worst over the run (pymoo: never loses the best ideal). + int feasibleCount = 0; for (int i = 0; i < n; i++) { - var f = population[i].Objectives; + if (population[i].IsFeasible) + feasibleCount++; + } + + // All-infeasible: keep the hyperplane from the last feasible generation. + if (feasibleCount == 0) + return ProjectWithFixedHyperplane(population); + + // Indices that may move ideal / worst. All-feasible is 0..n-1, in order, + // so the unconstrained update matches the previous whole-pool loop. + int[] poolIdx = new int[feasibleCount]; + int w = 0; + for (int i = 0; i < n; i++) + { + if (population[i].IsFeasible) + poolIdx[w++] = i; + } + + // Persistent ideal / worst over the run (pymoo: never loses the best ideal). + for (int p = 0; p < poolIdx.Length; p++) + { + var f = population[poolIdx[p]].Objectives; for (int j = 0; j < _m; j++) { if (f[j] < _ideal[j]) _ideal[j] = f[j]; @@ -83,22 +108,11 @@ public double[][] Normalize( } } - // Extreme-point search set: ND front when supplied, else all. - int[] ndIdx; - if (nonDominatedIndices is { Count: > 0 }) - { - ndIdx = new int[nonDominatedIndices.Count]; - for (int i = 0; i < nonDominatedIndices.Count; i++) - ndIdx[i] = nonDominatedIndices[i]; - } - else - { - ndIdx = new int[n]; - for (int i = 0; i < n; i++) ndIdx[i] = i; - } + // Extreme-point search set: ND front when supplied, else the hyperplane pool. + int[] ndIdx = ExtremeIndices(population, nonDominatedIndices, poolIdx); UpdateExtremePoints(population, ndIdx); - UpdateNadir(population, ndIdx); + UpdateNadir(population, ndIdx, poolIdx); var normalized = new double[n][]; for (int i = 0; i < n; i++) @@ -123,6 +137,95 @@ public double[][] Normalize( return normalized; } + /// + /// ND indices the caller supplied, restricted to the hyperplane pool, or the + /// pool itself when the caller did not supply a front. + /// + private static int[] ExtremeIndices( + IReadOnlyList population, + IReadOnlyList? nonDominatedIndices, + int[] poolIdx) + { + if (nonDominatedIndices is not { Count: > 0 }) + return (int[])poolIdx.Clone(); + + // All-feasible pool is every index, so a supplied front is copied as-is. + if (poolIdx.Length == population.Count) + { + var copy = new int[nonDominatedIndices.Count]; + for (int i = 0; i < nonDominatedIndices.Count; i++) + copy[i] = nonDominatedIndices[i]; + return copy; + } + + var kept = new List(nonDominatedIndices.Count); + for (int i = 0; i < nonDominatedIndices.Count; i++) + { + int idx = nonDominatedIndices[i]; + if ((uint)idx < (uint)population.Count && population[idx].IsFeasible) + kept.Add(idx); + } + + return kept.Count > 0 ? kept.ToArray() : (int[])poolIdx.Clone(); + } + + /// + /// Normalize with the persistent ideal and nadir. Does not write them. + /// Before any feasible generation, a temporary ideal / nadir is used only + /// for this projection so association stays finite. + /// + private double[][] ProjectWithFixedHyperplane(IReadOnlyList population) + { + int n = population.Count; + var origin = new double[_m]; + var nadir = new double[_m]; + bool persisted = !double.IsPositiveInfinity(_ideal[0]); + if (persisted) + { + Array.Copy(_ideal, origin, _m); + Array.Copy(_nadir, nadir, _m); + } + else + { + for (int j = 0; j < _m; j++) + { + origin[j] = double.PositiveInfinity; + nadir[j] = double.NegativeInfinity; + } + + for (int i = 0; i < n; i++) + { + var f = population[i].Objectives; + for (int j = 0; j < _m; j++) + { + if (f[j] < origin[j]) origin[j] = f[j]; + if (f[j] > nadir[j]) nadir[j] = f[j]; + } + } + + for (int j = 0; j < _m; j++) + { + if (nadir[j] - origin[j] <= 1e-6) + nadir[j] = origin[j] + 1.0; + } + } + + var normalized = new double[n][]; + for (int i = 0; i < n; i++) + { + normalized[i] = new double[_m]; + var f = population[i].Objectives; + for (int j = 0; j < _m; j++) + { + double denom = nadir[j] - origin[j]; + if (denom < 1e-12) denom = 1.0; + normalized[i][j] = (f[j] - origin[j]) / denom; + } + } + + return normalized; + } + /// Reset persistent state (tests / fresh run). Algorithm creates a new instance per run. public void Reset() { @@ -171,9 +274,9 @@ private void UpdateExtremePoints(IReadOnlyList population, int[] ndI _extremePoints = extreme; } - private void UpdateNadir(IReadOnlyList population, int[] ndIdx) + private void UpdateNadir(IReadOnlyList population, int[] ndIdx, int[] poolIdx) { - // Worst of front / population this generation. + // Worst of front / hyperplane pool this generation. var worstOfFront = new double[_m]; var worstOfPop = new double[_m]; for (int j = 0; j < _m; j++) @@ -181,9 +284,9 @@ private void UpdateNadir(IReadOnlyList population, int[] ndIdx) worstOfFront[j] = double.NegativeInfinity; worstOfPop[j] = double.NegativeInfinity; } - for (int i = 0; i < population.Count; i++) + for (int p = 0; p < poolIdx.Length; p++) { - var f = population[i].Objectives; + var f = population[poolIdx[p]].Objectives; for (int j = 0; j < _m; j++) if (f[j] > worstOfPop[j]) worstOfPop[j] = f[j]; } diff --git a/src/Unsga3/Operators/Selection/TournamentSelection.cs b/src/Unsga3/Operators/Selection/TournamentSelection.cs index 168012d..fa8daff 100644 --- a/src/Unsga3/Operators/Selection/TournamentSelection.cs +++ b/src/Unsga3/Operators/Selection/TournamentSelection.cs @@ -210,6 +210,8 @@ private List SelectTwoShuffledPasses( /// /// Recompute ranks and niche association for mating. This is a second normalization /// of the survivors. pymoo keeps the niche ids from environmental selection. + /// The normalizer keeps the feasible hyperplane: infeasible survivors do not move + /// ideal or nadir. /// public static void PrepareForSelection( IReadOnlyList population, diff --git a/src/Unsga3/Operators/Survival/NondominatedSortingSurvival.cs b/src/Unsga3/Operators/Survival/NondominatedSortingSurvival.cs index 4730ae2..21102f2 100644 --- a/src/Unsga3/Operators/Survival/NondominatedSortingSurvival.cs +++ b/src/Unsga3/Operators/Survival/NondominatedSortingSurvival.cs @@ -6,6 +6,10 @@ namespace Unsga3.Operators.Survival; /// /// NSGA-III environmental selection: fill complete fronts, then niche-preserve the last front. +/// When any member is infeasible, niching follows pymoo Survival.filter_infeasible +/// (Jain & Deb, NSGA-III Part II): niche the feasible subset, and fill a shortfall by +/// ascending constraint violation with no niche assignment. The hyperplane is updated +/// from feasible objectives whenever that subset is non-empty. /// public sealed class NondominatedSortingSurvival { @@ -28,6 +32,9 @@ public NondominatedSortingSurvival(ReferencePointManager references, Normalizati /// When a niche already has members, a random candidate in that niche is chosen (pymoo-style) /// if is provided; otherwise lowest index (deterministic). /// Normalization uses the first front only for extreme points (pymoo HyperplaneNormalization). + /// An all-feasible pool follows that path unchanged. A mixed or infeasible pool niches + /// only feasible members and fills any shortfall by ascending constraint violation + /// (original index breaks CV ties). Those fillers are not associated with a niche. /// public List Select( IReadOnlyList combined, @@ -39,11 +46,95 @@ public List Select( if (combined.Count <= targetSize) return combined.Select(i => i.Clone()).ToList(); - var fronts = NonDominatedSort.Sort(combined); + if (!AnyInfeasible(combined)) + return SelectByNiching(combined, targetSize, rng); + + var feasible = new List(); + var infeasibleIdx = new List(); + for (int i = 0; i < combined.Count; i++) + { + if (combined[i].IsFeasible) + feasible.Add(combined[i]); + else + infeasibleIdx.Add(i); + } + + // All infeasible: truncate by CV. Do not move the hyperplane. + if (feasible.Count == 0) + return TakeByAscendingCv(combined, infeasibleIdx, targetSize); + + int nFeasibleSurvive = Math.Min(feasible.Count, targetSize); + var selected = SelectByNiching(feasible, nFeasibleSurvive, rng); + + int need = targetSize - selected.Count; + if (need <= 0) + return selected; + + SortByAscendingCv(combined, infeasibleIdx); + + int take = Math.Min(need, infeasibleIdx.Count); + for (int k = 0; k < take; k++) + selected.Add(CloneWithoutNiche(combined[infeasibleIdx[k]])); + + return selected; + } + + private List SelectByNiching( + IReadOnlyList population, + int targetSize, + RandomProvider? rng) + { + var fronts = NonDominatedSort.Sort(population); // Extreme points from the ND front only — matches pymoo ReferenceDirectionSurvival. IReadOnlyList? nd = fronts.Count > 0 ? fronts[0] : null; - var normalized = _normalization.Normalize(combined, nd); - return SelectWithIndices(combined, fronts, targetSize, normalized, rng); + var normalized = _normalization.Normalize(population, nd); + if (population.Count <= targetSize) + return population.Select(i => i.Clone()).ToList(); + + return SelectWithIndices(population, fronts, targetSize, normalized, rng); + } + + private static bool AnyInfeasible(IReadOnlyList population) + { + for (int i = 0; i < population.Count; i++) + { + if (!population[i].IsFeasible) + return true; + } + + return false; + } + + private static void SortByAscendingCv(IReadOnlyList combined, List indices) + { + indices.Sort((i, j) => + { + int cmp = combined[i].ConstraintViolation.CompareTo(combined[j].ConstraintViolation); + return cmp != 0 ? cmp : i.CompareTo(j); + }); + } + + private static List TakeByAscendingCv( + IReadOnlyList combined, + List infeasibleIdx, + int targetSize) + { + SortByAscendingCv(combined, infeasibleIdx); + int take = Math.Min(targetSize, infeasibleIdx.Count); + var selected = new List(take); + for (int k = 0; k < take; k++) + selected.Add(CloneWithoutNiche(combined[infeasibleIdx[k]])); + + return selected; + } + + /// CV fill does not assign a niche (pymoo merges these after _do). + private static Individual CloneWithoutNiche(Individual source) + { + var clone = source.Clone(); + clone.AssociatedReference = -1; + clone.PerpendicularDistance = double.PositiveInfinity; + return clone; } private List SelectWithIndices( diff --git a/src/Unsga3/Problems/ConstrainedProblems.cs b/src/Unsga3/Problems/ConstrainedProblems.cs new file mode 100644 index 0000000..42aaef6 --- /dev/null +++ b/src/Unsga3/Problems/ConstrainedProblems.cs @@ -0,0 +1,110 @@ +using Unsga3.Core; + +namespace Unsga3.Problems; + +/// +/// OSY (Osyczka and Kundu). Six variables, two objectives, six inequalities. +/// Constraint scaling matches pymoo 0.6: each g ≥ 0 form is divided by its +/// constant, then negated so a positive value is a violation (g ≤ 0). +/// +public sealed class OsyProblem : ProblemBase +{ + public OsyProblem() + : base(6, 2, 6, new (double, double)[] + { + (0.0, 10.0), + (0.0, 10.0), + (1.0, 5.0), + (0.0, 6.0), + (1.0, 5.0), + (0.0, 10.0), + }) + { + } + + protected override void EvaluateCore(double[] x, double[] f, double[] g) + { + f[0] = -(25.0 * Sq(x[0] - 2.0) + + Sq(x[1] - 2.0) + + Sq(x[2] - 1.0) + + Sq(x[3] - 4.0) + + Sq(x[4] - 1.0)); + f[1] = Sq(x[0]) + Sq(x[1]) + Sq(x[2]) + Sq(x[3]) + Sq(x[4]) + Sq(x[5]); + + // pymoo stores G = -g_feas after the same divisors. + g[0] = -(x[0] + x[1] - 2.0) / 2.0; + g[1] = -(6.0 - x[0] - x[1]) / 6.0; + g[2] = -(2.0 - x[1] + x[0]) / 2.0; + g[3] = -(2.0 - x[0] + 3.0 * x[1]) / 2.0; + g[4] = -(4.0 - Sq(x[2] - 3.0) - x[3]) / 4.0; + g[5] = -(Sq(x[4] - 3.0) + x[5] - 4.0) / 4.0; + } + + private static double Sq(double v) => v * v; +} + +/// +/// TNK (Tanaka). Two variables on [0, π] × [1e-30, π], two inequalities in g ≤ 0 form. +/// The disk constraint uses pymoo's factor of two: +/// 2((x−0.5)² + (y−0.5)²) − 1 ≤ 0. +/// +public sealed class TnkProblem : ProblemBase +{ + public TnkProblem() + : base(2, 2, 2, new (double, double)[] + { + (0.0, Math.PI), + (1e-30, Math.PI), + }) + { + } + + protected override void EvaluateCore(double[] x, double[] f, double[] g) + { + f[0] = x[0]; + f[1] = x[1]; + double circle = x[0] * x[0] + x[1] * x[1] - 1.0 + - 0.1 * Math.Cos(16.0 * Math.Atan(x[0] / x[1])); + g[0] = -circle; + double dx = x[0] - 0.5; + double dy = x[1] - 0.5; + g[1] = 2.0 * (dx * dx + dy * dy) - 1.0; + } +} + +/// +/// C1-DTLZ1 (Jain & Deb, NSGA-III Part II). DTLZ1 objectives plus one inequality. +/// Feasible when 1 − f_M/0.6 − Σ_{i<M} f_i/0.5 ≥ 0. Stored in g ≤ 0 form, +/// so a positive value is the violation. Default k = 5 matches +/// (n = M + k − 1). The constraint is a barrier +/// in front of the linear front; it does not replace that front. +/// +public sealed class C1Dtlz1Problem : ProblemBase +{ + public C1Dtlz1Problem(int nObjectives = 3, int k = 5) + : base(nObjectives + k - 1, nObjectives, 1, DtlzHelper.UnitBounds(nObjectives + k - 1)) + { + if (nObjectives < 2) throw new ArgumentOutOfRangeException(nameof(nObjectives)); + if (k < 1) throw new ArgumentOutOfRangeException(nameof(k)); + } + + protected override void EvaluateCore(double[] x, double[] f, double[] g) + { + int m = NumberOfObjectives; + double gv = DtlzHelper.GDtlz1(x, m); + for (int i = 0; i < m; i++) + { + double val = 0.5 * (1.0 + gv); + for (int j = 0; j < m - i - 1; j++) + val *= x[j]; + if (i > 0) + val *= 1.0 - x[m - i - 1]; + f[i] = val; + } + + double weighted = f[m - 1] / 0.6; + for (int i = 0; i < m - 1; i++) + weighted += f[i] / 0.5; + g[0] = weighted - 1.0; + } +} diff --git a/src/Unsga3/Utilities/DasDennis.cs b/src/Unsga3/Utilities/DasDennis.cs index 568c095..2e66170 100644 --- a/src/Unsga3/Utilities/DasDennis.cs +++ b/src/Unsga3/Utilities/DasDennis.cs @@ -2,8 +2,9 @@ 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. +/// is a single layer and remains the default +/// (WithDasDennis, and the M≤3 oracles). adds +/// Deb & Jain Part I's inside layer for larger M. Survival accepts either array. /// public static class ReferenceDirections { @@ -28,6 +29,79 @@ public static double[][] DasDennis(int numberOfObjectives, int partitions) return points.ToArray(); } + /// + /// Outer Das–Dennis layer plus an inside layer scaled toward the centroid + /// (Deb & Jain, NSGA-III Part I). The inside point is + /// w' = s·w + (1−s)/M, which stays on the simplex. + /// defaults to 0.5. Points that repeat an + /// earlier direction (max-norm at most 1e-12), including the shared centroid, + /// are dropped. Single-layer stays the default; + /// current M≤3 oracles do not need a second layer. + /// + /// Simplex dimension M. + /// Divisions on the boundary layer. + /// Divisions on the inside layer. + /// Scale of the inside layer toward the centroid, in (0, 1]. + public static double[][] TwoLayerDasDennis( + int numberOfObjectives, + int outerPartitions, + int innerPartitions, + double innerScaling = 0.5) + { + if (!(innerScaling > 0.0) || innerScaling > 1.0) + throw new ArgumentOutOfRangeException( + nameof(innerScaling), innerScaling, "Inner scaling must be in (0, 1]."); + + var outer = DasDennis(numberOfObjectives, outerPartitions); + var inner = DasDennis(numberOfObjectives, innerPartitions); + var scaled = new double[inner.Length][]; + for (int i = 0; i < inner.Length; i++) + scaled[i] = ScaleTowardCentroid(inner[i], innerScaling); + + var kept = new List(outer.Length + scaled.Length); + kept.AddRange(outer); + for (int i = 0; i < scaled.Length; i++) + { + if (!IsDuplicate(scaled[i], kept)) + kept.Add(scaled[i]); + } + + return kept.ToArray(); + } + + /// w' = s·w + (1−s)/M. Sum is unchanged when is on the simplex. + internal static double[] ScaleTowardCentroid(double[] w, double scaling) + { + int m = w.Length; + double shift = (1.0 - scaling) / m; + var scaled = new double[m]; + for (int j = 0; j < m; j++) + scaled[j] = w[j] * scaling + shift; + return scaled; + } + + private static bool IsDuplicate(double[] point, List kept) + { + for (int i = 0; i < kept.Count; i++) + { + var other = kept[i]; + bool near = true; + for (int j = 0; j < point.Length; j++) + { + if (Math.Abs(point[j] - other[j]) > 1e-12) + { + near = false; + break; + } + } + + if (near) + return true; + } + + return false; + } + /// /// Choose partitions so the number of directions is at least /// (useful when population size is chosen first). diff --git a/tests/Unsga3.Tests/Unit/ConstrainedProblemTests.cs b/tests/Unsga3.Tests/Unit/ConstrainedProblemTests.cs new file mode 100644 index 0000000..d22f062 --- /dev/null +++ b/tests/Unsga3.Tests/Unit/ConstrainedProblemTests.cs @@ -0,0 +1,181 @@ +using Unsga3.Algorithm; +using Unsga3.Core; +using Unsga3.Operators.Survival; +using Unsga3.Problems; +using Unsga3.Utilities; + +namespace Unsga3.Tests.Unit; + +/// +/// Formulation checks and short runs for OSY, TNK, and C1-DTLZ1. +/// No IGD numbers: those wait on measured oracle runs. +/// +public class ConstrainedProblemTests +{ + [Fact] + public void Osy_matches_the_pymoo_constraint_scaling() + { + var feasible = Eval(new OsyProblem(), 5, 1, 5, 0, 5, 0); + Assert.Equal(-274.0, feasible.Objectives[0], 12); + Assert.Equal(76.0, feasible.Objectives[1], 12); + Assert.Equal(new[] { -2.0, 0.0, -3.0, 0.0, 0.0, 0.0 }, feasible.Constraints.ToList()); + Assert.Equal(0.0, feasible.ConstraintViolation); + Assert.True(feasible.IsFeasible); + + var infeasible = Eval(new OsyProblem(), 0, 0, 1, 0, 1, 0); + Assert.Equal(-120.0, infeasible.Objectives[0], 12); + Assert.Equal(2.0, infeasible.Objectives[1], 12); + Assert.Equal(1.0, infeasible.Constraints[0], 12); + Assert.Equal(1.0, infeasible.ConstraintViolation, 12); + Assert.False(infeasible.IsFeasible); + } + + [Fact] + public void Tnk_feasible_corner_and_interior_violation() + { + var feasible = Eval(new TnkProblem(), 1, 1); + Assert.Equal(1.0, feasible.Objectives[0], 12); + Assert.Equal(1.0, feasible.Objectives[1], 12); + Assert.Equal(-0.9, feasible.Constraints[0], 12); + Assert.Equal(0.0, feasible.Constraints[1], 12); + Assert.True(feasible.IsFeasible); + + var infeasible = Eval(new TnkProblem(), 0.1, 0.1); + Assert.Equal(1.08, infeasible.Constraints[0], 12); + Assert.True(infeasible.Constraints[1] < 0); + Assert.Equal(1.08, infeasible.ConstraintViolation, 12); + Assert.False(infeasible.IsFeasible); + } + + [Fact] + public void C1_dtlz1_center_is_feasible_and_the_origin_is_not() + { + var problem = new C1Dtlz1Problem(nObjectives: 3, k: 5); + Assert.Equal(7, problem.NumberOfVariables); + Assert.Equal(1, problem.NumberOfConstraints); + + var center = Eval(problem, Enumerable.Repeat(0.5, 7).ToArray()); + Assert.Equal(0.125, center.Objectives[0], 12); + Assert.Equal(0.125, center.Objectives[1], 12); + Assert.Equal(0.25, center.Objectives[2], 12); + Assert.Equal(-1.0 / 12.0, center.Constraints[0], 12); + Assert.True(center.IsFeasible); + + var origin = Eval(problem, new double[7]); + Assert.False(origin.IsFeasible); + Assert.True(origin.ConstraintViolation > 1); + } + + [Fact] + public void Osy_survival_uses_the_feasible_hyperplane_and_cv_fill() + { + var problem = new OsyProblem(); + var feasible = Eval(problem, 5, 1, 5, 0, 5, 0); + var lowerCv = Eval(problem, 0, 0, 1, 0, 1, 0); + var higherCv = Eval(problem, 0, 0, 1, 6, 1, 0); + Assert.True(lowerCv.ConstraintViolation < higherCv.ConstraintViolation); + // Both infeasible points have a smaller f2 than the feasible point. + Assert.True(lowerCv.Objectives[1] < feasible.Objectives[1]); + + var norm = new Normalization(2); + var survival = new NondominatedSortingSurvival( + new ReferencePointManager(ReferenceDirections.DasDennis(2, 1)), norm); + higherCv.AssociatedReference = 4; + + var selected = survival.Select( + new List { higherCv, feasible, lowerCv }, + targetSize: 2, + rng: null); + + Assert.True(selected[0].IsFeasible); + Assert.Equal(lowerCv.ConstraintViolation, selected[1].ConstraintViolation, 12); + Assert.Equal(-1, selected[1].AssociatedReference); + Assert.Equal(feasible.Objectives[0], norm.IdealPoint[0], 9); + Assert.Equal(feasible.Objectives[1], norm.IdealPoint[1], 9); + } + + [Fact] + public void Tnk_survival_does_not_take_the_infeasible_ideal() + { + var problem = new TnkProblem(); + var feasible = Eval(problem, 1, 1); + var closerToOrigin = Eval(problem, 0.1, 0.1); + var far = Eval(problem, 3, 3); + Assert.True(closerToOrigin.ConstraintViolation < far.ConstraintViolation); + + var norm = new Normalization(2); + var survival = new NondominatedSortingSurvival( + new ReferencePointManager(ReferenceDirections.DasDennis(2, 1)), norm); + var selected = survival.Select( + new List { far, closerToOrigin, feasible }, + targetSize: 2, + rng: null); + + Assert.True(selected[0].IsFeasible); + Assert.Equal(closerToOrigin.ConstraintViolation, selected[1].ConstraintViolation, 9); + Assert.Equal(1.0, norm.IdealPoint[0], 9); + Assert.Equal(1.0, norm.IdealPoint[1], 9); + } + + [Fact] + public void C1_dtlz1_survival_fills_from_the_smaller_violation() + { + var problem = new C1Dtlz1Problem(); + var feasible = Eval(problem, Enumerable.Repeat(0.5, 7).ToArray()); + var mild = Eval(problem, Enumerable.Repeat(0.2, 7).ToArray()); + var severe = Eval(problem, new double[7]); + Assert.True(mild.ConstraintViolation < severe.ConstraintViolation); + Assert.True(severe.Objectives[0] < feasible.Objectives[0]); + + var norm = new Normalization(3); + var survival = new NondominatedSortingSurvival( + new ReferencePointManager(ReferenceDirections.DasDennis(3, 1)), norm); + var selected = survival.Select( + new List { severe, feasible, mild }, + targetSize: 2, + rng: null); + + Assert.True(selected[0].IsFeasible); + Assert.Equal(mild.ConstraintViolation, selected[1].ConstraintViolation, 8); + Assert.Equal(feasible.Objectives[0], norm.IdealPoint[0], 9); + Assert.Equal(feasible.Objectives[1], norm.IdealPoint[1], 9); + Assert.Equal(feasible.Objectives[2], norm.IdealPoint[2], 9); + } + + [Theory] + [InlineData("osy")] + [InlineData("tnk")] + [InlineData("c1")] + public void Seeded_run_keeps_a_feasible_member(string name) + { + IProblem problem = name switch + { + "osy" => new OsyProblem(), + "tnk" => new TnkProblem(), + "c1" => new C1Dtlz1Problem(), + _ => throw new ArgumentOutOfRangeException(nameof(name)), + }; + Individual seed = name switch + { + "osy" => new Individual(new[] { 5.0, 1, 5, 0, 5, 0 }, problem.NumberOfObjectives, problem.NumberOfConstraints), + "tnk" => new Individual(new[] { 1.0, 1.0 }, problem.NumberOfObjectives, problem.NumberOfConstraints), + "c1" => new Individual(Enumerable.Repeat(0.5, problem.NumberOfVariables).ToArray(), problem.NumberOfObjectives, problem.NumberOfConstraints), + _ => throw new ArgumentOutOfRangeException(nameof(name)), + }; + + var dirs = ReferenceDirections.DasDennis(problem.NumberOfObjectives, problem.NumberOfObjectives == 2 ? 4 : 2); + var algo = new Unsga3Algorithm(dirs, populationSize: 12, seed: 1); + var result = algo.Run(problem, maxGenerations: 3, initialPopulation: new[] { seed }); + + Assert.Equal(12, result.FinalPopulation.Count); + Assert.Contains(result.FinalPopulation, ind => ind.IsFeasible); + Assert.All(result.FinalPopulation, ind => Assert.Equal(problem.NumberOfObjectives, ind.Objectives.Length)); + } + + private static Individual Eval(IProblem problem, params double[] x) + { + var ind = new Individual(x, problem.NumberOfObjectives, problem.NumberOfConstraints); + problem.Evaluate(ind); + return ind; + } +} diff --git a/tests/Unsga3.Tests/Unit/ConstrainedSurvivalTests.cs b/tests/Unsga3.Tests/Unit/ConstrainedSurvivalTests.cs new file mode 100644 index 0000000..c8f278f --- /dev/null +++ b/tests/Unsga3.Tests/Unit/ConstrainedSurvivalTests.cs @@ -0,0 +1,144 @@ +using Unsga3.Algorithm; +using Unsga3.Core; +using Unsga3.Operators.Selection; +using Unsga3.Operators.Survival; +using Unsga3.Utilities; + +namespace Unsga3.Tests.Unit; + +/// +/// pymoo Survival.filter_infeasible: niche feasible members, fill a shortfall +/// by ascending constraint violation, and build the hyperplane from feasible objectives. +/// +public class ConstrainedSurvivalTests +{ + [Fact] + public void Niches_only_feasible_members_when_they_fill_the_population() + { + var norm = new Normalization(2); + var survival = Survival(norm); + var infeasible = Mark(0.0, 0.0, cv: 0.25, tag: 9); + infeasible.AssociatedReference = 7; + + var selected = survival.Select( + new List + { + Mark(0.2, 0.9, tag: 1), + Mark(0.9, 0.2, tag: 2), + Mark(0.4, 0.4, tag: 3), + infeasible, + }, + targetSize: 2, + rng: null); + + Assert.Equal(2, selected.Count); + Assert.All(selected, ind => Assert.True(ind.IsFeasible)); + Assert.DoesNotContain(selected, ind => ind.Variables[0] == 9); + Assert.Equal(0.2, norm.IdealPoint[0], 12); + Assert.Equal(0.2, norm.IdealPoint[1], 12); + } + + [Fact] + public void Shortfall_is_filled_by_ascending_cv_without_a_niche() + { + var norm = new Normalization(2); + var survival = Survival(norm); + var low = Mark(9, 9, cv: 0.2, tag: 2); + var mid = Mark(0, 0, cv: 0.5, tag: 3); + var high = Mark(1, 8, cv: 1.5, tag: 4); + low.AssociatedReference = 3; + mid.AssociatedReference = 4; + high.AssociatedReference = 5; + + var selected = survival.Select( + new List { high, Mark(0.2, 0.8, tag: 1), low, mid }, + targetSize: 3, + rng: null); + + Assert.Equal(new[] { 1.0, 2.0, 3.0 }, selected.Select(s => s.Variables[0])); + Assert.True(selected[0].IsFeasible); + Assert.Equal(-1, selected[1].AssociatedReference); + Assert.Equal(-1, selected[2].AssociatedReference); + Assert.Equal(double.PositiveInfinity, selected[1].PerpendicularDistance); + // Infeasible (0, 0) must not set the ideal. The only feasible point is (0.2, 0.8). + Assert.Equal(0.2, norm.IdealPoint[0], 12); + Assert.Equal(0.8, norm.IdealPoint[1], 12); + } + + [Fact] + public void Equal_cv_fill_keeps_the_earlier_index() + { + var survival = Survival(new Normalization(2)); + var earlier = Mark(4, 4, cv: 1, tag: 8); + var later = Mark(1, 1, cv: 1, tag: 9); + var selected = survival.Select( + new List { earlier, Mark(0.3, 0.7, tag: 1), later }, + targetSize: 2, + rng: null); + + Assert.Equal(1.0, selected[0].Variables[0]); + Assert.Equal(8.0, selected[1].Variables[0]); + } + + [Fact] + public void All_infeasible_truncates_by_cv_and_leaves_the_hyperplane() + { + var norm = new Normalization(2); + norm.Normalize(new List { Mark(1, 3), Mark(3, 1) }); + double ideal0 = norm.IdealPoint[0]; + double nadir0 = norm.NadirPoint[0]; + + var survival = Survival(norm); + var selected = survival.Select( + new List + { + Mark(0, 0, cv: 3, tag: 3), + Mark(4, 1, cv: 0.4, tag: 1), + Mark(2, 2, cv: 1.2, tag: 2), + }, + targetSize: 2, + rng: null); + + Assert.Equal(new[] { 1.0, 2.0 }, selected.Select(s => s.Variables[0])); + Assert.All(selected, ind => Assert.Equal(-1, ind.AssociatedReference)); + Assert.Equal(ideal0, norm.IdealPoint[0], 12); + Assert.Equal(nadir0, norm.NadirPoint[0], 12); + Assert.Equal(1.0, norm.IdealPoint[1], 12); + } + + [Fact] + public void PrepareForSelection_on_an_all_infeasible_pop_does_not_set_ideal() + { + var norm = new Normalization(2); + var refs = new ReferencePointManager(ReferenceDirections.DasDennis(2, 1)); + var pop = new List + { + Mark(1.0, 2.0, cv: 1.0), + Mark(3.0, 0.5, cv: 2.0), + }; + + TournamentSelection.PrepareForSelection(pop, refs, norm); + + Assert.True(double.IsPositiveInfinity(norm.IdealPoint[0])); + Assert.True(double.IsPositiveInfinity(norm.IdealPoint[1])); + Assert.All(pop, ind => Assert.False(double.IsNaN(ind.PerpendicularDistance))); + Assert.All(pop, ind => Assert.InRange(ind.AssociatedReference, 0, refs.Count - 1)); + } + + private static NondominatedSortingSurvival Survival(Normalization norm) + { + var refs = new ReferencePointManager(ReferenceDirections.DasDennis(2, 1)); + return new NondominatedSortingSurvival(refs, norm); + } + + private static Individual Mark(double f1, double f2, double cv = 0, double tag = 0) + { + var ind = new Individual(1, 2, 1); + ind.Variables[0] = tag; + ind.Objectives[0] = f1; + ind.Objectives[1] = f2; + ind.Constraints[0] = cv; + ind.RefreshConstraintViolation(); + return ind; + } +} diff --git a/tests/Unsga3.Tests/Unit/DasDennisTests.cs b/tests/Unsga3.Tests/Unit/DasDennisTests.cs index 654c4e7..d2db3af 100644 --- a/tests/Unsga3.Tests/Unit/DasDennisTests.cs +++ b/tests/Unsga3.Tests/Unit/DasDennisTests.cs @@ -1,3 +1,5 @@ +using Unsga3.Algorithm; +using Unsga3.Problems; using Unsga3.Utilities; namespace Unsga3.Tests.Unit; @@ -43,4 +45,76 @@ public void Single_objective_is_unit_scalar() Assert.Single(pts); Assert.Equal(1.0, pts[0][0]); } + + [Fact] + public void Two_layer_prefixes_the_outer_simplex_and_scales_the_inside_layer() + { + var outer = ReferenceDirections.DasDennis(3, 2); + var inner = ReferenceDirections.DasDennis(3, 1); + var layers = ReferenceDirections.TwoLayerDasDennis(3, 2, 1); + + Assert.Equal(outer.Length + inner.Length, layers.Length); + for (int i = 0; i < outer.Length; i++) + AssertSamePoint(outer[i], layers[i]); + + // First inner direction (0, 0, 1), scale 0.5 → (1/6, 1/6, 2/3). + var expected = ReferenceDirections.ScaleTowardCentroid(inner[0], 0.5); + AssertSamePoint(expected, layers[outer.Length]); + Assert.All(layers, w => + { + Assert.Equal(1.0, w.Sum(), 9); + Assert.All(w, v => Assert.InRange(v, -1e-12, 1.0 + 1e-12)); + }); + } + + [Fact] + public void Two_layer_drops_the_shared_centroid() + { + // p = 3 includes (1/3, 1/3, 1/3). Scaling that point by 0.5 leaves it in place. + int outerCount = ReferenceDirections.Count(3, 3); + var layers = ReferenceDirections.TwoLayerDasDennis(3, 3, 3, innerScaling: 0.5); + Assert.Equal(outerCount * 2 - 1, layers.Length); + + int centroids = layers.Count(w => + Math.Abs(w[0] - 1.0 / 3.0) < 1e-12 + && Math.Abs(w[1] - 1.0 / 3.0) < 1e-12 + && Math.Abs(w[2] - 1.0 / 3.0) < 1e-12); + Assert.Equal(1, centroids); + } + + [Fact] + public void Two_layer_scaling_of_one_collapses_to_the_outer_layer() + { + var outer = ReferenceDirections.DasDennis(3, 4); + var layers = ReferenceDirections.TwoLayerDasDennis(3, 4, 4, innerScaling: 1.0); + Assert.Equal(outer.Length, layers.Length); + for (int i = 0; i < outer.Length; i++) + AssertSamePoint(outer[i], layers[i]); + } + + [Fact] + public void Two_layer_directions_are_accepted_by_survival() + { + // M≤3 oracles keep single-layer DasDennis. This only checks the factory + // output is a legal direction array for one short unconstrained run. + var dirs = ReferenceDirections.TwoLayerDasDennis(3, 3, 2); + var algo = new Unsga3Algorithm(dirs, populationSize: 16, seed: 1); + var result = algo.Run(new Dtlz2Problem(nObjectives: 3, k: 5), maxGenerations: 2); + Assert.Equal(16, result.FinalPopulation.Count); + Assert.NotEmpty(result.NonDominatedSolutions); + } + + [Fact] + public void Two_layer_rejects_a_non_positive_scale() + { + Assert.Throws( + () => ReferenceDirections.TwoLayerDasDennis(3, 2, 1, innerScaling: 0.0)); + } + + private static void AssertSamePoint(double[] expected, double[] actual) + { + Assert.Equal(expected.Length, actual.Length); + for (int j = 0; j < expected.Length; j++) + Assert.Equal(expected[j], actual[j], 12); + } } diff --git a/tests/Unsga3.Tests/Unit/NonDominatedSortTests.cs b/tests/Unsga3.Tests/Unit/NonDominatedSortTests.cs index 481f850..454ba7b 100644 --- a/tests/Unsga3.Tests/Unit/NonDominatedSortTests.cs +++ b/tests/Unsga3.Tests/Unit/NonDominatedSortTests.cs @@ -44,6 +44,35 @@ public void Constraint_domination_prefers_feasible() Assert.Equal(1, NonDominatedSort.CompareConstraintDominated(infeas, feas)); } + [Fact] + public void Equal_constraint_violation_is_mutually_non_dominated() + { + // (0, 0) Pareto-dominates (1, 1), but Deb's constraint-domination stops + // when the violations are equal. Both land on the first front. + var better = Make(0, 0); + var worse = Make(1, 1); + better.Constraints[0] = 1.5; + worse.Constraints[0] = 1.5; + better.RefreshConstraintViolation(); + worse.RefreshConstraintViolation(); + + Assert.Equal(0, NonDominatedSort.CompareConstraintDominated(better, worse)); + Assert.Equal(0, NonDominatedSort.CompareConstraintDominated(worse, better)); + Assert.False(better.ConstraintDominates(worse)); + Assert.False(worse.ConstraintDominates(better)); + + var fronts = NonDominatedSort.Sort(new List { better, worse }); + Assert.Single(fronts); + Assert.Equal(2, fronts[0].Count); + Assert.Equal(0, better.Rank); + Assert.Equal(0, worse.Rank); + + // A strictly smaller violation still dominates, objectives aside. + worse.Constraints[0] = 1.5 + 1e-6; + worse.RefreshConstraintViolation(); + Assert.Equal(-1, NonDominatedSort.CompareConstraintDominated(better, worse)); + } + private static Individual Make(double f1, double f2) { var ind = new Individual(1, 2, 1); diff --git a/tests/Unsga3.Tests/Unit/NormalizationTests.cs b/tests/Unsga3.Tests/Unit/NormalizationTests.cs index b0d730c..c97d19e 100644 --- a/tests/Unsga3.Tests/Unit/NormalizationTests.cs +++ b/tests/Unsga3.Tests/Unit/NormalizationTests.cs @@ -111,12 +111,11 @@ public void Collapsed_span_sets_nadir_to_ideal_plus_one() } [Fact] - public void Infeasible_origin_sets_ideal_from_the_whole_pool() + public void Feasible_members_set_the_hyperplane() { - // 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. + // Feasible (1, 1) beside infeasible (0, 0). Ideal and worst come from the + // feasible point only (Jain & Deb Part II / pymoo filter_infeasible). + // One feasible point collapses, so nadir = ideal + 1. var feasible = new Individual(1, 2, 1); feasible.Objectives[0] = 1; feasible.Objectives[1] = 1; @@ -139,19 +138,43 @@ public void Infeasible_origin_sets_ideal_from_the_whole_pool() 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)); + Assert.Equal(1.0, norm.IdealPoint[0], 12); + Assert.Equal(1.0, norm.IdealPoint[1], 12); + Assert.Equal(2.0, norm.NadirPoint[0], 12); + Assert.Equal(2.0, norm.NadirPoint[1], 12); + Assert.Equal(0.0, normalized[0][0], 9); + Assert.Equal(0.0, normalized[0][1], 9); + Assert.Equal(-1.0, normalized[1][0], 9); + Assert.Equal(-1.0, normalized[1][1], 9); + } + + [Fact] + public void All_infeasible_population_does_not_move_the_hyperplane() + { + var norm = new Normalization(2); + var feasible = new List { Make2(1, 3), Make2(3, 1) }; + norm.Normalize(feasible); + Assert.Equal(1.0, norm.IdealPoint[0], 12); + Assert.Equal(1.0, norm.IdealPoint[1], 12); + double nadir0 = norm.NadirPoint[0]; + double nadir1 = norm.NadirPoint[1]; + + var infeasible = new Individual(1, 2, 1); + infeasible.Objectives[0] = 0; + infeasible.Objectives[1] = 0; + infeasible.Constraints[0] = 2; + infeasible.RefreshConstraintViolation(); + var worse = new Individual(1, 2, 1); + worse.Objectives[0] = 9; + worse.Objectives[1] = 9; + worse.Constraints[0] = 4; + worse.RefreshConstraintViolation(); + + norm.Normalize(new List { infeasible, worse }); + Assert.Equal(1.0, norm.IdealPoint[0], 12); + Assert.Equal(1.0, norm.IdealPoint[1], 12); + Assert.Equal(nadir0, norm.NadirPoint[0], 12); + Assert.Equal(nadir1, norm.NadirPoint[1], 12); } private static Individual Make1(double a)