feat(moead): MOEA/D-DE, differential evolution in place of the crossover (Li and Zhang 2009)

AGENTS.md

@@ -760,10 +760,11 @@ fn main() -> genoxide::Result<()> {
 | `Nsga2` | `Nsga2::builder(real, objectives)` | Non-dominated sorting, crowding distance |
 | `Spea2`, `SmsEmoa` | same as `Nsga2` | SPEA2's truncation spreads the front evenly; SMS-EMOA keeps the largest hypervolume contributions; both cost more per generation |
 | `Moead` | `Moead::builder(real, objectives, multi::das_dennis::<M>(divisions))` | A subproblem per weight vector; Tchebycheff, or `multi::Decomposition::Pbi { theta: 5.0 }` for 3+ objectives; cheap |
+| MOEA/D-DE | `Moead` with `.crossover(multi::DifferentialEvolutionCrossover::new(0.5, 1.0)?)` | Li and Zhang's: each child is its subproblem's solution plus `F (a − b)`; for complicated Pareto sets (genes linked to one another) and constrained problems (with `.neighbor_mating(0.2)`), where SBX stalls; slower than SBX on ZDT and DTLZ |
 | `Nsga3` | `Nsga3::builder(real, objectives, multi::das_dennis::<3>(12))` | 3+ objectives; population defaults to the number of directions (91); SBX η 30 |
 
 - Objective counts are typed: `[f64; 3]` for 2 objectives doesn't compile.
-- Constrained dominance: feasible first, then the smaller violation (`multi::dominates`). `Stop::stagnation` counts generations whose front has no solution that the previous front didn't dominate or equal: with a front larger than the population, it may never fire, so add `Stop::generations`; `Stop::target` isn't available.
+- Constrained dominance: feasible first, then the smaller violation (`multi::dominates`); MOEA/D compares a child with a solution the same way, then by the subproblem's value. `Stop::stagnation` counts generations whose front has no solution that the previous front didn't dominate or equal: with a front larger than the population, it may never fire, so add `Stop::generations`; `Stop::target` isn't available.
 - `multi::non_dominated_sort`, `multi::crowding_distance`; `multi::ParetoArchive::new(objectives)` with `.on_generation(|snapshot| archive.update(snapshot))` keeps every non-dominated solution.
 - Test problems: `multi::problems::{Zdt1, Zdt2, Zdt3, Zdt4, Zdt6}::new(n)`, `Zdt5` (a `Binary` genome of 80 bits, not in `all()`), `{Dtlz1, …, Dtlz7}::<M>::new(n)` (the technical report's numbering), `{Wfg1, …, Wfg9}::<M>::new(k, l)` (k position parameters, a multiple of M − 1; l distance parameters, even for WFG2 and WFG3; `default()`: k = 4 for 2 objectives, 2(M − 1) for more, l = 20), `{FonsecaFleming, Kursawe}::new(n)`, `Schaffer1`, `Schaffer2`, `Poloni`, `Viennet1`/`2`/`3` (3 objectives), and the constrained `Bnh`, `Srn`, `Tnk`, `Osy`, `Constr` (fitness `([f64; 2], violation)`), all minimized, for `MultiEngine::new(algorithm, problem)`. The `multi::problems::MultiProblem<M>` trait gives `representation()`, `optimal_front(points)` (`Option`: `None` for KUR, POL, VNT2, VNT3, DTLZ5, DTLZ6 with 4 or more objectives, and WFG3 with 3 or more), `ideal_point()`, `nadir_point()`, `constraints(&x)` (`g <= 0`), `reference()`; `multi::problems::all::<M>()` lists them as `Box<dyn DynMultiProblem<M>>`.
 - Constrained test problems of tunable difficulty: `multi::problems::{Ctp1, …, Ctp8}` (two objectives, two variables, disconnected fronts, fronts of separate points, infeasible bands and tunnels) and `{C1Dtlz1, C1Dtlz3, C2Dtlz2, ConvexC2Dtlz2, C3Dtlz1, C3Dtlz4}::<M>::new(n)` (Jain and Deb's constrained DTLZ; `C1Dtlz3` and `ConvexC2Dtlz2` take the paper's radius for 3, 5, 8, 10 and 15 objectives, else `with_radius(n, r)`), all with `optimal_front`, and in `all()`.

docs/features.md

@@ -78,7 +78,7 @@ What genoxide has on main; [docs.rs](https://docs.rs/genoxide) documents the lat
 
 ## Multi-objective
 
-- **Algorithms:** NSGA-II, NSGA-III, SPEA2, MOEA/D, SMS-EMOA.
+- **Algorithms:** NSGA-II, NSGA-III, SPEA2, MOEA/D, SMS-EMOA; MOEA/D-DE (Li and Zhang, 2009), differential evolution in place of MOEA/D's crossover, for Pareto sets whose genes are linked and for constrained problems, checked against the paper's F2.
 - **With:** constraints, duplicate elimination, a Pareto archive.
 - **Indicators:** hypervolume, IGD, IGD+, GD, spread.
 - **Test problems** (`multi::problems`): ZDT1-6 (ZDT5 on bit strings), DTLZ1-7, WFG1-9 (any number of objectives, checked against the authors' toolkit), Schaffer's two, Fonseca and Fleming's, Kursawe's, Poloni's and Viennet's three, and the constrained BNH, SRN, TNK, OSY and CONSTR, each with its optimal front where it's known and its reference, in Rust and Python.

docs/optimization-plan.md

@@ -232,7 +232,7 @@ tested, and orders the work in batches.
 | Derivatives | none | supplied gradients, finite differences, gradient-based methods |
 | Constraints | Deb's rules and penalties over an aggregate violation; per-constraint values in `problems::Problem::constraints` | methods that use each constraint's value and gradient (SQP, augmented Lagrangian, interior point) |
 | Expensive functions | `AsyncEngine`, `Batch`, parallel evaluation | surrogate models (Gaussian processes) and Bayesian optimization |
-| Multi-objective | NSGA-II, NSGA-III, SPEA2, MOEA/D, SMS-EMOA | ParEGO, EHVI |
+| Multi-objective | NSGA-II, NSGA-III, SPEA2, MOEA/D (with SBX or differential evolution, MOEA/D-DE), SMS-EMOA | ParEGO, EHVI |
 | Linear algebra | a symmetric eigendecomposition inside CMA-ES (tred2 and tql2 of JAMA) | Cholesky, QR, triangular solves, a dense QP solver |
 
 ### 1.2 Local derivative-free methods

python/README.md

@@ -169,7 +169,7 @@ A solution that couldn't be scored has NaN objective values and a NaN violation.
 
 - `Nsga2` spreads the front by crowding distance, which works poorly beyond 2 or 3 objectives.
 - `Nsga3` spreads it along reference directions instead.
-- `Moead` solves a single-objective subproblem per weight vector.
+- `Moead` solves a single-objective subproblem per weight vector. With `crossover=gx.DifferentialEvolutionCrossover()` it is MOEA/D-DE (Li and Zhang, 2009): each child is its subproblem's solution moved by the difference of two parents, for Pareto sets whose genes are linked. Constraint violations are compared first (Deb's rules); on constrained problems, a lower `neighbor_mating` such as 0.2 keeps the population spread.
 - `das_dennis(objectives, divisions)` gives evenly spread directions or weights for `Nsga3` and `Moead`, a row each: 91 for 3 objectives and 12 divisions.
 - `Spea2` keeps an archive of the best solutions, the non-dominated ones first, truncated by the distance to their nearest neighbors.
 - `Nsga2`, `Nsga3`, `Spea2` and `SmsEmoa` drop a child that equals a member of the population or an earlier child, and breed another. `eliminate_duplicates=False` keeps copies.
@@ -220,7 +220,7 @@ Operators:
 - **Selection:** `Tournament(size)`, `Rank(pressure)`, `Roulette()`, `StochasticUniversalSampling()`, `Truncation(fraction)`, `RandomSelection()`; against bloat (trees that grow without getting better), selections that see a genome's size, a tree's nodes: `DoubleTournament(fitness_size, parsimony)` (7 and 1.4, Luke and Panait's best), `LexicographicTournament(size, bucket_ratio=None)` (of equal fitness, the smaller wins), `Tarpeian(select, rate)` (genomes larger than the mean count as invalid with probability `rate`)
 - **Crossover:**
   - any list genome: `UniformCrossover()`, `PointCrossover(points)`, `NoCrossover()`
-  - real genomes: `SimulatedBinaryCrossover(eta)`, `BlendCrossover(alpha)`, `ArithmeticCrossover()`
+  - real genomes: `SimulatedBinaryCrossover(eta)`, `BlendCrossover(alpha)`, `ArithmeticCrossover()`; with `Moead` only, `DifferentialEvolutionCrossover(f, cr, repair)` (0.5, 1 and `"bounce"` by default: MOEA/D-DE, for Pareto sets whose genes are linked)
   - permutations: `OrderCrossover()`, `PartiallyMappedCrossover()`, `CycleCrossover()`, `EdgeRecombinationCrossover()`
 - **Mutation:**
   - binary genomes: `BitFlip(rate=... | count=...)`
@@ -805,6 +805,7 @@ Some names differ:
 | `De(control={"f": 0.5, "cr": 0.9})` | `.control(de::Control::Fixed { f: 0.5, cr: 0.9 })`; `{"min_f", "max_f", "cr"}` for `Dither`, `{"c"}` for `Jade`, `{"memory"}` for `Shade` |
 | `De(restarts="never")`, `De(restarts={"tolerance": 1e-12, "patience": 200})` | `.restarts(de::Restarts::Never)`, `.restarts(de::Restarts::OnStagnation { tolerance: 1e-12, patience: 200 })` |
 | `Pbi(theta)` | `Decomposition::Pbi { theta }` |
+| `Moead(crossover=gx.DifferentialEvolutionCrossover(f, cr, repair="random"))` | `.crossover(DifferentialEvolutionCrossover::new(f, cr)?.with_repair(moead::Repair::Random))` |
 | `gx.gp.Gp(primitives, init=gx.gp.Full((2, 6)))` | `Gp::builder(set).init(Init::Full { depths: 2..=6 }).build()?` |
 | `gx.gp.PointMutation(rate=...)`, `gx.gp.ConstantMutation(sigma)` | `PointMutation::per_node(rate)`, `ConstantMutation::gaussian(sigma)` |
 | `gx.gp.Mutations([(0.5, gx.gp.SubtreeMutation()), (0.5, gx.gp.HoistMutation())])` | `Mutations::builder().subtree(0.5).hoist(0.5).build()?` |

python/genoxide/__init__.py

@@ -101,6 +101,7 @@
     "PartiallyMappedCrossover",
     "CycleCrossover",
     "EdgeRecombinationCrossover",
+    "DifferentialEvolutionCrossover",
     # mutation
     "BitFlip",
     "UniformMutation",
@@ -644,6 +645,38 @@ def _describe(self) -> dict[str, Any]:
         return {"type": "edge_recombination"}
 
 
+@dataclass(frozen=True)
+class DifferentialEvolutionCrossover:
+    """MOEA/D-DE's differential evolution (Li and Zhang, 2009), in place of a crossover: for
+    :class:`Moead` on real genomes only. A subproblem's child is its own solution ``x`` moved by
+    the difference of the two parents: each gene becomes ``x + f * (a - b)`` with probability
+    ``cr``. A child whose parents came from the whole population may replace solutions anywhere
+    in it. Suited to Pareto sets whose genes are linked, where SBX stalls.
+
+    ``f`` is greater than 0 and at most 2, 0.5 by default; ``cr`` is 0 to 1, 1 by default (Li and
+    Zhang's settings, with polynomial mutation with eta 20 at a rate of 1 / n). ``repair`` brings
+    back a gene that leaves its bounds: ``"bounce"`` (the default), a random value between ``x``
+    and the bound it crossed, which reaches genes on their bounds; or ``"random"``, a random
+    value anywhere within the bounds, the repair the paper's text describes."""
+
+    f: float = 0.5
+    cr: float = 1.0
+    repair: Literal["bounce", "random"] = "bounce"
+
+    def _describe(self) -> dict[str, Any]:
+        if self.repair not in ("bounce", "random"):
+            raise ValueError(
+                f'DifferentialEvolutionCrossover.repair is "bounce" or "random", not '
+                f"{self.repair!r}"
+            )
+        return {
+            "type": "differential_evolution",
+            "f": _number("DifferentialEvolutionCrossover.f", self.f),
+            "cr": _number("DifferentialEvolutionCrossover.cr", self.cr),
+            "repair": self.repair,
+        }
+
+
 Crossover = Union[
     UniformCrossover,
     PointCrossover,
@@ -655,6 +688,7 @@ def _describe(self) -> dict[str, Any]:
     PartiallyMappedCrossover,
     CycleCrossover,
     EdgeRecombinationCrossover,
+    DifferentialEvolutionCrossover,
     "gp.SubtreeCrossover",
     "gp.OnePointCrossover",
 ]
@@ -5555,9 +5589,11 @@ class Moead(_MultiObjective):
     :class:`Tchebycheff` (the default) or :class:`Pbi`, which spreads fronts of 3 or more
     objectives well. Each generation, every subproblem gets a child, one of the crossover's two.
 
-    For real genomes, SBX with eta 20 is the usual crossover. The fitness function is as for
-    :class:`Nsga2`. MOEA/D has no ``eliminate_duplicates``: it replaces its neighbors one child at
-    a time.
+    For real genomes, SBX with eta 20 is the usual crossover; :class:`DifferentialEvolutionCrossover`
+    makes it MOEA/D-DE (Li and Zhang, 2009), for Pareto sets whose genes are linked. The fitness
+    function is as for :class:`Nsga2`; with a constraint violation, a feasible solution beats an
+    infeasible one and the smaller violation wins (Deb's rules). MOEA/D has no
+    ``eliminate_duplicates``: it replaces its neighbors one child at a time.
 
     Parameters
     ----------
@@ -5569,7 +5605,8 @@ class Moead(_MultiObjective):
         At least 2 weight vectors, a row each with a value per objective: finite, non-negative
         and not all 0. Their number is the population size.
     crossover : a crossover
-        How pairs of parents are combined. It must fit the genome.
+        How pairs of parents are combined. It must fit the genome. For real genomes,
+        :class:`DifferentialEvolutionCrossover` too.
     mutation : a mutation
         How children are changed. It must fit the genome.
     neighbors : int, default 20

python/src/config.rs

@@ -499,6 +499,21 @@ pub enum Crossover {
     },
     /// Trees: subtrees exchanged at a point of the common region.
     OnePoint {},
+    /// Real genomes, MOEA/D only: the subproblem's solution moved by `f` times the difference of
+    /// two parents in each gene with probability `cr` (MOEA/D-DE).
+    DifferentialEvolution {
+        f: f64,
+        cr: f64,
+        repair: Option<DeRepair>,
+    },
+}
+
+/// How MOEA/D-DE brings back a gene that leaves its bounds.
+#[derive(Clone, Copy, Debug, Deserialize)]
+#[serde(rename_all = "snake_case")]
+pub enum DeRepair {
+    Bounce,
+    Random,
 }
 
 /// A mutation: `rate` changes each gene with that probability, `count` exactly that many genes.

python/src/operators.rs

@@ -4,6 +4,8 @@
 use crate::config;
 use crate::errors::{setting, setting_named};
 use genoxide::genome::{AdaptiveReal, Binary, Genome, Integer, Permutation, Real, Representation};
+use genoxide::multi::DifferentialEvolutionCrossover;
+use genoxide::multi::moead::Repair;
 use genoxide::operator::{
     ArithmeticCrossover, BitFlip, BlendCrossover, Crossover, CycleCrossover, DoubleTournament,
     EdgeRecombinationCrossover, GaussianMutation, InsertionMutation, InversionMutation,
@@ -160,6 +162,7 @@ fn crossover_name(crossover: config::Crossover) -> &'static str {
         config::Crossover::EdgeRecombination {} => "EdgeRecombinationCrossover",
         config::Crossover::Subtree { .. } => "SubtreeCrossover",
         config::Crossover::OnePoint {} => "OnePointCrossover",
+        config::Crossover::DifferentialEvolution { .. } => "DifferentialEvolutionCrossover",
     }
 }
 
@@ -290,6 +293,10 @@ impl RealCrossover {
                 Ok(Self::Blend(setting(BlendCrossover::new(alpha))?))
             }
             config::Crossover::Arithmetic {} => Ok(Self::Arithmetic(ArithmeticCrossover::new())),
+            config::Crossover::DifferentialEvolution { .. } => Err(
+                "DifferentialEvolutionCrossover works with Moead only; use SimulatedBinaryCrossover"
+                    .to_string(),
+            ),
             _ => Err(wrong_crossover(
                 crossover,
                 "real",
@@ -322,6 +329,26 @@ impl Crossover<Real> for RealCrossover {
     }
 }
 
+/// MOEA/D-DE's differential evolution, from its settings.
+pub fn differential_evolution(
+    crossover: config::Crossover,
+) -> Result<DifferentialEvolutionCrossover> {
+    match crossover {
+        config::Crossover::DifferentialEvolution { f, cr, repair } => {
+            let names = [
+                ("f", "DifferentialEvolutionCrossover.f"),
+                ("cr", "DifferentialEvolutionCrossover.cr"),
+            ];
+            let crossover = setting_named(DifferentialEvolutionCrossover::new(f, cr), &names)?;
+            Ok(match repair {
+                None | Some(config::DeRepair::Bounce) => crossover.with_repair(Repair::Bounce),
+                Some(config::DeRepair::Random) => crossover.with_repair(Repair::Random),
+            })
+        }
+        _ => Err("not a differential evolution crossover".to_string()),
+    }
+}
+
 /// A crossover of permutations.
 #[derive(Clone, Debug, serde::Serialize, serde::Deserialize)]
 pub enum OrderCrossovers {