diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml index 7ee6898a5..77d4e406f 100644 --- a/.github/workflows/python.yml +++ b/.github/workflows/python.yml @@ -144,7 +144,7 @@ jobs: - name: Cache intersphinx inventories uses: actions/cache@v5 with: - path: python/doc/_build/doctrees/ + path: python/doc/_build/doctrees/__intersphinx_cache__/ key: intersphinx-${{ hashFiles('python/doc/source/conf.py') }} restore-keys: intersphinx- diff --git a/bibkeys.txt b/bibkeys.txt index 0a0ff5285..74fa1abbc 100644 --- a/bibkeys.txt +++ b/bibkeys.txt @@ -1,5 +1,6 @@ AugBadBroZit2009gecco BeuFonLopPaqVah09:tec +BezLopStu2015emode BezLopStu2017assessment BezLopStu2017emo BinGinRou2015gaupar @@ -14,24 +15,30 @@ Deb02nsga2 DenZha2019approxhv DiaLop2020ejor DubLopStu2011amai +EibSmi2003 EmmFon2011emo +EshSch1992 FanWan1994numtheory FonGueLopPaq2011emo FonPaqLop06:hypervolume Gla2017fast +Goldberg89 GruFon2009:emaa Grunert01 GueFon2017hv4d GueFonPaq2021hv HanJas1998 HerWer1987tabucol +Hol75 IshHeSha2019regular IshMasTanNoj2015igd Jen03 JohAraMcGSch1991 KnoCor2002cec +KukLam2005:gde3 KunLucPre1975jacm LacKlaFon2017box +LiZha2009:moead-de Lop2026hvapprox LopPaqStu09emaa LopVerDreDoe2025 @@ -39,9 +46,12 @@ LopezIbanezPhD Molchanov2005theory Muller1959sphere Rob2018unreasonable +RobFil05demo RubMel1998simulation SchEsqLarCoe2012tec SchKer2025r2v2 +StoPri1997:de +Syswerda89 VelLam1998gp WuAza2001metrics ZhoZhaJin2009igdx diff --git a/c/Makefile b/c/Makefile index 81922b8d6..919cc075b 100644 --- a/c/Makefile +++ b/c/Makefile @@ -56,6 +56,7 @@ BINDIR:=$(abspath $(BINDIR)) ## Define source files SRCS = avl.c \ cmdline.c \ + crossover.c \ dominatedsets.c \ eaf3d.c \ eaf.c \ @@ -90,6 +91,7 @@ HDRS = avl.h \ cmdline.h \ common.h \ config.h \ + crossover.h \ cvector.h \ eaf.h \ epsilon.h \ @@ -208,6 +210,7 @@ $(EXE_FILES): cmdline.o io.o avl.o: avl.h cmdline.o: cmdline.h io.h +crossover.o: crossover.h rng.h mt19937/mt19937.h dominatedsets.o: cmdline.h io.h nondominated.h epsilon.h eaf.o: eaf.h io.h bit_array.h cvector.h eaf3d.o: eaf.h io.h bit_array.h cvector.h avl.h diff --git a/c/NEWS.md b/c/NEWS.md index 317b10c1d..ebcb229c8 100644 --- a/c/NEWS.md +++ b/c/NEWS.md @@ -3,6 +3,9 @@ ## 0.20 * radixsort.h, insort.h: New. + * crossover.c, crossover.h: New crossover operators (uniform, single-point, + blend-alpha, single arithmetic, local, bounded exponential, differential + evolution and parent-centric). * main-nondominated.c: Renamed from nondominated.c. * hv and hvapprox executables now use -U instead of -u for --union to be consistent with nondominated. diff --git a/c/crossover.c b/c/crossover.c new file mode 100644 index 000000000..2ab92074a --- /dev/null +++ b/c/crossover.c @@ -0,0 +1,439 @@ +/************************************************************************* + + crossover.c + + Crossover operators for real-valued (and, for uniform and single-point + crossover, any numerical) decision vectors. See crossover.h for the + conventions shared by all functions. + + --------------------------------------------------------------------- + + Copyright (C) 2026 + Bhupinder Saini + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + + ---------------------------------------------------------------------- + + +*************************************************************************/ +#include +#include // memcpy(), memmove() +#include +#include "common.h" +#include "crossover.h" +#include "rng.h" + +// parent_centric_crossover() treats an index parent this close to the +// centroid of its three parents as degenerate. +#define PCX_DEGENERATE 1e-12 + +static inline double +clamp_to_bounds(double x, double lower, double upper) +{ + return fmin(fmax(x, lower), upper); +} + +/* Rows I and I + 1 of PARENTS are crossed over and their offspring go to rows + I and I + 1 of OFFSPRING. If I + 1 == N, row I is crossed over with row 0 + and its second offspring goes to SCRATCH, so it is discarded. */ +static inline void +get_pair(const double * parents, double * offspring, double * scratch, + size_t i, size_t n, size_t nvar, + const double ** p1, const double ** p2, double ** off1, double ** off2) +{ + const bool last = (i + 1 == n); + *p1 = parents + i * nvar; + *p2 = parents + (last ? 0 : i + 1) * nvar; + *off1 = offspring + i * nvar; + *off2 = last ? scratch : offspring + (i + 1) * nvar; +} + +static inline double * +scratch_new(size_t n, size_t nvar) +{ + return (n % 2 == 1) ? MOOCORE_MALLOC(nvar, double) : NULL; +} + +// Offspring of a pair that is not crossed over. +static inline void +copy_parents(double * off1, double * off2, const double * p1, const double * p2, + size_t nvar) +{ + memcpy(off1, p1, nvar * sizeof(*off1)); + memcpy(off2, p2, nvar * sizeof(*off2)); +} + +/* Draw K distinct indices in [0, N) that also differ from the NEXCLUDED + indices in EXCLUDED. EXCLUDED must be sorted in ascending order and have + room for K more values. Each draw picks among the indices not yet excluded + and shifts the result past the excluded ones, so there is no rejection + loop. */ +static void +sample_distinct(rng_state * rng, size_t n, size_t * excluded, size_t nexcluded, + size_t k, size_t * chosen) +{ + for (size_t c = 0; c < k; c++, nexcluded++) { + ASSUME(n > nexcluded); + size_t d = rng_integer(rng, n - nexcluded); + size_t pos = 0; + while (pos < nexcluded && excluded[pos] <= d) { + d++; + pos++; + } + // Insert D at POS to keep EXCLUDED sorted. + memmove(excluded + pos + 1, excluded + pos, + (nexcluded - pos) * sizeof(*excluded)); + excluded[pos] = d; + chosen[c] = d; + } +} + +void +uniform_crossover(double * restrict offspring, const double * restrict parents, + size_t n, size_t nvar, double swap_prob, double pair_prob, + uint32_t seed) +{ + ASSUME(n >= 2 && nvar >= 1); + rng_state * rng = rng_new(seed); + double * scratch = scratch_new(n, nvar); + for (size_t i = 0; i < n; i += 2) { + const double *p1, *p2; + double *off1, *off2; + get_pair(parents, offspring, scratch, i, n, nvar, + &p1, &p2, &off1, &off2); + const bool recombine = rng_random(rng) < pair_prob; + for (size_t j = 0; j < nvar; j++) { + const bool swap = recombine && rng_random(rng) < swap_prob; + off1[j] = swap ? p2[j] : p1[j]; + off2[j] = swap ? p1[j] : p2[j]; + } + } + free(scratch); + rng_free(rng); +} + +void +single_point_crossover(double * restrict offspring, + const double * restrict parents, + size_t n, size_t nvar, double pair_prob, uint32_t seed) +{ + ASSUME(n >= 2 && nvar >= 2); + rng_state * rng = rng_new(seed); + double * scratch = scratch_new(n, nvar); + for (size_t i = 0; i < n; i += 2) { + const double *p1, *p2; + double *off1, *off2; + get_pair(parents, offspring, scratch, i, n, nvar, + &p1, &p2, &off1, &off2); + if (!(rng_random(rng) < pair_prob)) { + copy_parents(off1, off2, p1, p2, nvar); + continue; + } + // Cut point in [1, nvar - 1], so each offspring takes at least one + // variable from each parent. + const size_t cut = 1 + rng_integer(rng, nvar - 1); + memcpy(off1, p1, cut * sizeof(*off1)); + memcpy(off1 + cut, p2 + cut, (nvar - cut) * sizeof(*off1)); + memcpy(off2, p2, cut * sizeof(*off2)); + memcpy(off2 + cut, p1 + cut, (nvar - cut) * sizeof(*off2)); + } + free(scratch); + rng_free(rng); +} + +void +blend_alpha_crossover(double * restrict offspring, + const double * restrict parents, + size_t n, size_t nvar, const double * restrict lower, + const double * restrict upper, double alpha, + double pair_prob, uint32_t seed) +{ + ASSUME(n >= 2 && nvar >= 1); + ASSUME(alpha >= 0); + rng_state * rng = rng_new(seed); + for (size_t i = 0; i < n; i += 2) { + const double * p1 = parents + i * nvar; + const double * p2 = parents + ((i + 1 < n) ? i + 1 : 0) * nvar; + double * off = offspring + (i / 2) * nvar; + // A pair that is not crossed over yields a copy of its first parent. + if (!(rng_random(rng) < pair_prob)) { + memcpy(off, p1, nvar * sizeof(*off)); + continue; + } + for (size_t j = 0; j < nvar; j++) { + const double lo = fmin(p1[j], p2[j]); + const double hi = fmax(p1[j], p2[j]); + const double span = hi - lo; + // Widen [lo, hi] by ALPHA * SPAN on each side, clip it to the + // bounds and sample uniformly from it. + const double a = fmax(lo - alpha * span, lower[j]); + const double b = fmin(hi + alpha * span, upper[j]); + off[j] = a + rng_random(rng) * (b - a); + } + } + rng_free(rng); +} + +void +single_arithmetic_crossover(double * restrict offspring, + const double * restrict parents, + size_t n, size_t nvar, double pair_prob, + uint32_t seed) +{ + ASSUME(n >= 2 && nvar >= 1); + rng_state * rng = rng_new(seed); + double * scratch = scratch_new(n, nvar); + for (size_t i = 0; i < n; i += 2) { + const double *p1, *p2; + double *off1, *off2; + get_pair(parents, offspring, scratch, i, n, nvar, + &p1, &p2, &off1, &off2); + const bool recombine = rng_random(rng) < pair_prob; + const size_t k = rng_integer(rng, nvar); + memcpy(off1, p1, nvar * sizeof(*off1)); + memcpy(off2, p2, nvar * sizeof(*off2)); + if (recombine) + off1[k] = off2[k] = 0.5 * (p1[k] + p2[k]); + } + free(scratch); + rng_free(rng); +} + +void +local_crossover(double * restrict offspring, const double * restrict parents, + size_t n, size_t nvar, double pair_prob, uint32_t seed) +{ + ASSUME(n >= 2 && nvar >= 1); + rng_state * rng = rng_new(seed); + double * scratch = scratch_new(n, nvar); + for (size_t i = 0; i < n; i += 2) { + const double *p1, *p2; + double *off1, *off2; + get_pair(parents, offspring, scratch, i, n, nvar, + &p1, &p2, &off1, &off2); + if (!(rng_random(rng) < pair_prob)) { + copy_parents(off1, off2, p1, p2, nvar); + continue; + } + for (size_t j = 0; j < nvar; j++) { + const double a = rng_random(rng); + off1[j] = a * p1[j] + (1 - a) * p2[j]; + off2[j] = (1 - a) * p1[j] + a * p2[j]; + } + } + free(scratch); + rng_free(rng); +} + +/* Offspring value of parent value X. It moves X by BETA * SPAN. BETA is + negative or positive with probability 1/2 each. On each side, it follows + an exponential distribution with scale SCALE, truncated at the bound on + that side, so the offspring stays within [LOWER, UPPER]. U and R are + uniform in [0, 1). */ +static inline double +bex_offspring(double x, double span, double lower, double upper, double scale, + double u, double r) +{ + double beta; + if (r < 0.5) { + const double e = exp((lower - x) / (scale * span)); + // 1 - U is in (0, 1], which avoids log(0) when LOWER is -INFINITY. + beta = scale * log(e + (1 - u) * (1 - e)); + } else { + const double e = exp((x - upper) / (scale * span)); + beta = -scale * log(1 - u * (1 - e)); + } + // Rounding may step just outside the bounds. + return clamp_to_bounds(x + beta * span, lower, upper); +} + +void +bounded_exponential_crossover(double * restrict offspring, + const double * restrict parents, + size_t n, size_t nvar, + const double * restrict lower, + const double * restrict upper, double scale, + double swap_prob, double pair_prob, + uint32_t seed) +{ + ASSUME(n >= 2 && nvar >= 1); + ASSUME(scale > 0); + rng_state * rng = rng_new(seed); + double * scratch = scratch_new(n, nvar); + for (size_t i = 0; i < n; i += 2) { + const double *p1, *p2; + double *off1, *off2; + get_pair(parents, offspring, scratch, i, n, nvar, + &p1, &p2, &off1, &off2); + if (!(rng_random(rng) < pair_prob)) { + copy_parents(off1, off2, p1, p2, nvar); + continue; + } + for (size_t j = 0; j < nvar; j++) { + // The distance between the parents scales the move of both + // offspring, whichever parent is smaller. + const double span = fabs(p2[j] - p1[j]); + if (span == 0) { + off1[j] = p1[j]; + off2[j] = p2[j]; + continue; + } + // Both offspring use the same random numbers. + const double u = rng_random(rng); + const double r = rng_random(rng); + off1[j] = bex_offspring(p1[j], span, lower[j], upper[j], + scale, u, r); + off2[j] = bex_offspring(p2[j], span, lower[j], upper[j], + scale, u, r); + if (swap_prob > 0 && rng_random(rng) < swap_prob) { + const double tmp = off1[j]; + off1[j] = off2[j]; + off2[j] = tmp; + } + } + } + free(scratch); + rng_free(rng); +} + +void +de_rand_1_bin_crossover(double * restrict offspring, + const double * restrict population, + size_t n, size_t nvar, const size_t * restrict targets, + size_t ntargets, const double * restrict lower, + const double * restrict upper, double scaling_factor, + double crossover_rate, uint32_t seed) +{ + ASSUME(n >= 4 && nvar >= 1); + rng_state * rng = rng_new(seed); + for (size_t t = 0; t < ntargets; t++) { + ASSUME(targets[t] < n); + const double * x = population + targets[t] * nvar; + size_t excluded[4] = { targets[t] }; + size_t r[3]; + sample_distinct(rng, n, excluded, 1, 3, r); + const double * x1 = population + r[0] * nvar; + const double * x2 = population + r[1] * nvar; + const double * x3 = population + r[2] * nvar; + // The offspring takes at least the variable JRAND from the mutant. + const size_t jrand = rng_integer(rng, nvar); + double * off = offspring + t * nvar; + for (size_t j = 0; j < nvar; j++) { + const double v = (j == jrand || rng_random(rng) < crossover_rate) + ? x1[j] + scaling_factor * (x2[j] - x3[j]) + : x[j]; + off[j] = clamp_to_bounds(v, lower[j], upper[j]); + } + } + rng_free(rng); +} + +void +de_current_1_crossover(double * restrict offspring, + const double * restrict population, + size_t n, size_t nvar, const size_t * restrict targets, + size_t ntargets, const double * restrict lower, + const double * restrict upper, double scaling_factor, + double crossover_rate, uint32_t seed) +{ + ASSUME(n >= 2 && nvar >= 1); + rng_state * rng = rng_new(seed); + for (size_t t = 0; t < ntargets; t++) { + ASSUME(targets[t] < n); + const double * x = population + targets[t] * nvar; + // The two donors differ from each other, but either may be the target. + size_t excluded[2]; + size_t r[2]; + sample_distinct(rng, n, excluded, 0, 2, r); + const double * x1 = population + r[0] * nvar; + const double * x2 = population + r[1] * nvar; + double * off = offspring + t * nvar; + for (size_t j = 0; j < nvar; j++) { + const double v = (rng_random(rng) < crossover_rate) + ? x[j] + scaling_factor * (x1[j] - x2[j]) + : x[j]; + off[j] = clamp_to_bounds(v, lower[j], upper[j]); + } + } + rng_free(rng); +} + +void +parent_centric_crossover(double * restrict offspring, + const double * restrict population, + size_t n, size_t nvar, + const size_t * restrict targets, size_t ntargets, + const double * restrict lower, + const double * restrict upper, + double sigma_zeta, double sigma_eta, uint32_t seed) +{ + ASSUME(n >= 3 && nvar >= 1); + rng_state * rng = rng_new(seed); + double * direction = MOOCORE_MALLOC(3 * nvar, double); + double * unit = direction + nvar; + double * orthogonal = unit + nvar; + for (size_t t = 0; t < ntargets; t++) { + ASSUME(targets[t] < n); + const double * x = population + targets[t] * nvar; + size_t excluded[3] = { targets[t] }; + size_t r[2]; + sample_distinct(rng, n, excluded, 1, 2, r); + const double * others[2] = { population + r[0] * nvar, + population + r[1] * nvar }; + + // Direction from the centroid of the three parents to the index + // parent X. + double norm = 0; + for (size_t j = 0; j < nvar; j++) { + direction[j] = x[j] - (x[j] + others[0][j] + others[1][j]) / 3; + norm += direction[j] * direction[j]; + } + norm = sqrt(norm); + const bool degenerate = norm < PCX_DEGENERATE; + const double divisor = degenerate ? 1 : norm; + for (size_t j = 0; j < nvar; j++) + unit[j] = direction[j] / divisor; + + // Mean distance from the other two parents to the line through X + // along DIRECTION. + double spread = 0; + for (int k = 0; k < 2; k++) { + const double * o = others[k]; + double along = 0; + for (size_t j = 0; j < nvar; j++) + along += (o[j] - x[j]) * unit[j]; + double dist = 0; + for (size_t j = 0; j < nvar; j++) { + const double perp = (o[j] - x[j]) - along * unit[j]; + dist += perp * perp; + } + spread += sqrt(dist); + } + spread /= 2; + + // Gaussian noise with its component along DIRECTION removed. This + // has the same distribution as summing independent normal steps + // along an orthonormal basis of the subspace orthogonal to DIRECTION. + double proj = 0; + for (size_t j = 0; j < nvar; j++) { + orthogonal[j] = sigma_eta * rng_standard_normal(rng); + proj += orthogonal[j] * unit[j]; + } + if (!degenerate) { + for (size_t j = 0; j < nvar; j++) + orthogonal[j] -= proj * unit[j]; + } + const double zeta = sigma_zeta * rng_standard_normal(rng); + + double * off = offspring + t * nvar; + for (size_t j = 0; j < nvar; j++) + off[j] = clamp_to_bounds( + x[j] + zeta * direction[j] + spread * orthogonal[j], + lower[j], upper[j]); + } + free(direction); + rng_free(rng); +} diff --git a/c/crossover.h b/c/crossover.h new file mode 100644 index 000000000..791fb8c09 --- /dev/null +++ b/c/crossover.h @@ -0,0 +1,95 @@ +/************************************************************************* + + crossover.h + + --------------------------------------------------------------------- + + Copyright (C) 2026 + Bhupinder Saini + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at https://mozilla.org/MPL/2.0/. + + ---------------------------------------------------------------------- + + +*************************************************************************/ +#ifndef CROSSOVER_H_ +#define CROSSOVER_H_ + +#include +#include // uint32_t +#include "libmoocore-config.h" + +// C++ needs to know that types and declarations are C, not C++. +BEGIN_C_DECLS + +/* + All matrices are stored by rows, one decision vector per row, with NVAR + variables each. OFFSPRING must be allocated by the caller. + + Functions that take PARENTS (N rows) cross over rows 2i and 2i+1. If N is + odd, the last row is crossed over with row 0. Operators that produce two + offspring per pair write them to rows 2i and 2i+1 of OFFSPRING (N rows), + dropping the second offspring of the last pair when N is odd. + blend_alpha_crossover() produces one offspring per pair, so OFFSPRING has + (N + 1) / 2 rows. + + Functions that take POPULATION (N rows) produce one offspring per row index in + TARGETS (NTARGETS values in [0, N)), so OFFSPRING has NTARGETS rows. They + choose the other parents at random from POPULATION. + + LOWER and UPPER give the bounds of each variable. Use -INFINITY and + INFINITY for unbounded variables. +*/ + +MOOCORE_API void uniform_crossover( + double * restrict offspring, const double * restrict parents, + size_t n, size_t nvar, double swap_prob, double pair_prob, uint32_t seed); + +MOOCORE_API void single_point_crossover( + double * restrict offspring, const double * restrict parents, + size_t n, size_t nvar, double pair_prob, uint32_t seed); + +MOOCORE_API void blend_alpha_crossover( + double * restrict offspring, const double * restrict parents, + size_t n, size_t nvar, const double * restrict lower, + const double * restrict upper, double alpha, double pair_prob, + uint32_t seed); + +MOOCORE_API void single_arithmetic_crossover( + double * restrict offspring, const double * restrict parents, + size_t n, size_t nvar, double pair_prob, uint32_t seed); + +MOOCORE_API void local_crossover( + double * restrict offspring, const double * restrict parents, + size_t n, size_t nvar, double pair_prob, uint32_t seed); + +MOOCORE_API void bounded_exponential_crossover( + double * restrict offspring, const double * restrict parents, + size_t n, size_t nvar, const double * restrict lower, + const double * restrict upper, double scale, double swap_prob, + double pair_prob, uint32_t seed); + +MOOCORE_API void de_rand_1_bin_crossover( + double * restrict offspring, const double * restrict population, + size_t n, size_t nvar, const size_t * restrict targets, size_t ntargets, + const double * restrict lower, const double * restrict upper, + double scaling_factor, double crossover_rate, uint32_t seed); + +MOOCORE_API void de_current_1_crossover( + double * restrict offspring, const double * restrict population, + size_t n, size_t nvar, const size_t * restrict targets, size_t ntargets, + const double * restrict lower, const double * restrict upper, + double scaling_factor, double crossover_rate, uint32_t seed); + +MOOCORE_API void parent_centric_crossover( + double * restrict offspring, const double * restrict population, + size_t n, size_t nvar, const size_t * restrict targets, size_t ntargets, + const double * restrict lower, const double * restrict upper, + double sigma_zeta, double sigma_eta, uint32_t seed); + +END_C_DECLS + +#endif // CROSSOVER_H_ diff --git a/c/rng.h b/c/rng.h index 3b5a65cc7..892140218 100644 --- a/c/rng.h +++ b/c/rng.h @@ -37,6 +37,14 @@ rng_uniform(rng_state * rng, double low, double high) return low + (high - low) * rng_random(rng); } +/* Returns an integer between [0, n) */ +static inline size_t +rng_integer(rng_state * rng, size_t n) +{ + assert(n > 0); + return STATIC_CAST(size_t, rng_random(rng) * STATIC_CAST(double, n)); +} + double rng_standard_normal(rng_state *rng); void rng_bivariate_normal_fill(rng_state * rng, diff --git a/python/benchmarks/README.md b/python/benchmarks/README.md index 002b7186c..b9e505b63 100644 --- a/python/benchmarks/README.md +++ b/python/benchmarks/README.md @@ -4,4 +4,4 @@ Install dependencies with: python3 -m pip install moocore[benchmarks] -Benchmark results: https://multi-objective.github.io/moocore/python/#benchmarks +Benchmark results: https://multi-objective.github.io/moocore/python/benchmarks/index.html diff --git a/python/benchmarks/bench.py b/python/benchmarks/bench.py index e9e7c409c..7e68188a9 100644 --- a/python/benchmarks/bench.py +++ b/python/benchmarks/bench.py @@ -262,6 +262,7 @@ def plots( relative=False, xlabel="n", show_ci=False, + add_name=True, ): """Plots the results after running the benchmark one or more times. @@ -366,7 +367,8 @@ def single_plot(results, ylabel, file_suffix): ax.fill_between(x, y1, y2, where=valid, alpha=0.2) plt.title(f"({self.cpu_model})", fontsize=10) - plt.suptitle(f"{title} for {self.name}", fontsize=12) + suptitle = f"{title} for {self.name}" if add_name else title + plt.suptitle(suptitle, fontsize=12) save2png(f"{file_prefix}_bench-{self.name}-{file_suffix}.png") # CPU times. diff --git a/python/benchmarks/bench_crossover.py b/python/benchmarks/bench_crossover.py new file mode 100644 index 000000000..86a6c9b24 --- /dev/null +++ b/python/benchmarks/bench_crossover.py @@ -0,0 +1,449 @@ +"""Crossover Benchmarks +==================== + +This example benchmarks the crossover operators in ``moocore`` against the +equivalent operators in other packages, when they exist. + +Each package receives the parents in its own data structure, which is created +before the timing starts. When a package expects the caller to choose the +other parents of DE or PCX, choosing them is timed, because ``moocore`` does it +inside the call. + +""" + +import random + +from bench import Bench + +import numpy as np +import polars as pl +import matplotlib.pyplot as plt +import moocore + +from deap import tools as deap_tools +from desdeo.emo.operators import crossover as desdeo_crossover +from desdeo.problem.schema import ( + Objective, + Problem as DesdeoProblem, + Variable, + VariableTypeEnum, +) +from desdeo.tools.patterns import Publisher +from jmetal.core.solution import FloatSolution +from jmetal.operator.crossover import ( + BLXAlphaCrossover, + DifferentialEvolutionCrossover, +) +from pymoo.algorithms.moo.nsde import VariantDE +from pymoo.core.population import Population +from pymoo.core.problem import Problem as PymooProblem +from pymoo.operators.crossover.dex import DEX +from pymoo.operators.crossover.pcx import PCX +from pymoo.operators.crossover.pntx import SinglePointCrossover +from pymoo.operators.crossover.ux import UniformCrossover + +LOWER, UPPER = 0.0, 1.0 +# Number of parents for a fixed number of variables, and vice versa. +N_PARENTS = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000] +N_VARS = [2, 5, 10, 20, 50, 100, 200, 500, 1000] +FIXED_NVAR = 30 +FIXED_NPARENTS = 100 + +rng = np.random.default_rng(42) +# DEAP and jMetalPy draw from the random module. +random.seed(42) + + +def distinct_indices(targets, n, k, exclude_targets=True): + """Draw k distinct row indices in [0, n) for each target.""" + m = len(targets) + if exclude_targets: + excluded = targets[:, None].copy() + else: + excluded = np.empty((m, 0), dtype=int) + chosen = np.empty((m, k), dtype=int) + for j in range(k): + drawn = rng.integers(0, n - excluded.shape[1], size=m) + for col in range(excluded.shape[1]): + drawn += drawn >= excluded[:, col] + chosen[:, j] = drawn + excluded = np.sort(np.column_stack([excluded, drawn]), axis=1) + return chosen + + +def check_offspring(a, b, what, n, name): + """Check that each package returns as many offspring as moocore.""" + if hasattr(b, "to_numpy"): # DESDEO + b = b.to_numpy() + elif len(b) and hasattr(b[0], "variables"): # jMetalPy + b = [s.variables for s in b] + b = np.asarray(b, dtype=float) + assert b.shape == a.shape, ( + f"In {name}, n={n}, {what} returned shape {b.shape}, but moocore returned {a.shape}" + ) + + +# pymoo: a Population and an unbounded Problem, so that pymoo does not repair +# the offspring. Offspring are clipped as in moocore when needed. +def pymoo_setup(x): + return PymooProblem(n_var=x.shape[1]), Population.new("X", x) + + +def pymoo_bounded_setup(x): + problem = PymooProblem(n_var=x.shape[1], xl=LOWER, xu=UPPER) + return problem, Population.new("X", x) + + +def pymoo_pairs(op): + def run(problem, pop): + parents = np.arange(len(pop)).reshape(-1, 2) + return op.do(problem, pop, parents, random_state=rng).get("X") + + return run + + +pymoo_dex = DEX(F=0.5, CR=1.0, at_least_once=False) +pymoo_pcx = PCX(eta=0.1, zeta=0.1, prob=1.0) +pymoo_variant_de = VariantDE( + variant="DE/rand/1/bin", CR=0.9, F=0.5, gamma=None, de_repair="to-bounds" +) + + +def pymoo_de_rand_1_bin(problem, pop): + return pymoo_variant_de.do(problem, pop, len(pop), random_state=rng).get( + "X" + ) + + +def pymoo_de_current_1(problem, pop): + targets = np.arange(len(pop)) + donors = distinct_indices(targets, len(pop), 2, exclude_targets=False) + parents = np.column_stack([targets, donors]) + off = pymoo_dex.do(problem, pop, parents, random_state=rng) + return np.clip(off.get("X"), LOWER, UPPER) + + +def pymoo_parent_centric(problem, pop): + targets = np.arange(len(pop)) + parents = np.column_stack([targets, distinct_indices(targets, len(pop), 2)]) + off = pymoo_pcx.do(problem, pop, parents, random_state=rng) + return np.clip(off.get("X"), LOWER, UPPER) + + +# DEAP: lists of floats. The operators change both individuals in place. +def deap_setup(x): + return x.tolist() + + +def deap_pairs(fun, *args): + def run(ind): + for i in range(0, len(ind) - 1, 2): + ind[i], ind[i + 1] = fun(ind[i], ind[i + 1], *args) + return ind + + return run + + +# jMetalPy: lists of FloatSolution. +def jmetalpy_setup(x): + lower, upper = [LOWER] * x.shape[1], [UPPER] * x.shape[1] + solutions = [] + for row in x: + s = FloatSolution(lower, upper, number_of_objectives=1) + s.variables = row.tolist() + solutions.append(s) + return solutions + + +def jmetalpy_pairs(op): + def run(solutions): + offspring = [] + for i in range(0, len(solutions) - 1, 2): + offspring.extend(op.execute([solutions[i], solutions[i + 1]])) + return offspring + + return run + + +jmetalpy_de_rand = DifferentialEvolutionCrossover(CR=0.9, F=0.5) +jmetalpy_de_current = DifferentialEvolutionCrossover(CR=1.0, F=0.5) + + +def jmetalpy_de_rand_1_bin(solutions): + # execute([a, b, c]) returns c + F (a - b), crossed over with the target. + targets = np.arange(len(solutions)) + offspring = [] + for i, (r1, r2, r3) in zip( + targets, distinct_indices(targets, len(solutions), 3) + ): + jmetalpy_de_rand.current_individual = solutions[i] + offspring.extend( + jmetalpy_de_rand.execute( + [solutions[r2], solutions[r3], solutions[r1]] + ) + ) + return offspring + + +def jmetalpy_de_current_1(solutions): + targets = np.arange(len(solutions)) + donors = distinct_indices(targets, len(solutions), 2, exclude_targets=False) + offspring = [] + for i, (r1, r2) in zip(targets, donors): + jmetalpy_de_current.current_individual = solutions[i] + offspring.extend( + jmetalpy_de_current.execute( + [solutions[r1], solutions[r2], solutions[i]] + ) + ) + return offspring + + +# DESDEO: a polars DataFrame and an operator built for the problem. +def desdeo_problem(nvar, binary=False): + if binary: + kind, lower, upper = VariableTypeEnum.binary, 0, 1 + else: + kind, lower, upper = VariableTypeEnum.real, LOWER, UPPER + variables = [ + Variable( + name=f"x_{i}", + symbol=f"x_{i}", + variable_type=kind, + lowerbound=lower, + upperbound=upper, + initial_value=lower, + ) + for i in range(1, nvar + 1) + ] + objective = Objective(name="f_1", symbol="f_1", func="x_1", maximize=False) + return DesdeoProblem( + name="bench", + description="Crossover benchmark", + variables=variables, + objectives=[objective], + ) + + +def desdeo_setup(cls, binary=False, **kwargs): + def setup(x): + problem = desdeo_problem(x.shape[1], binary=binary) + op = cls( + problem=problem, + publisher=Publisher(), + verbosity=0, + seed=42, + **kwargs, + ) + symbols = [v.symbol for v in problem.get_flattened_variables()] + population = pl.DataFrame(x, schema=symbols, orient="row") + return op, population, list(range(len(x))) + + return setup + + +def desdeo_run(op, population, rows): + return op.do(population=population, to_mate=rows) + + +def moocore_blend_alpha(x): + # moocore produces one offspring per pair. Call it twice to produce as + # many offspring as the other packages. + return np.vstack( + [ + moocore.blend_alpha_crossover(x, lower=LOWER, upper=UPPER, seed=rng) + for _ in range(2) + ] + ) + + +# For each operator: (data, {package: (setup, function)}). +real = "real" +binary = "binary" +operators = { + "uniform": ( + real, + { + "moocore": (None, lambda x: moocore.uniform_crossover(x, seed=rng)), + "pymoo": (pymoo_setup, pymoo_pairs(UniformCrossover(prob=1.0))), + "deap": (deap_setup, deap_pairs(deap_tools.cxUniform, 0.5)), + "desdeo": ( + desdeo_setup(desdeo_crossover.UniformCrossover), + desdeo_run, + ), + }, + ), + "single_point": ( + binary, + { + "moocore": ( + None, + lambda x: moocore.single_point_crossover(x, seed=rng), + ), + "pymoo": (pymoo_setup, pymoo_pairs(SinglePointCrossover(prob=1.0))), + "deap": (deap_setup, deap_pairs(deap_tools.cxOnePoint)), + "desdeo": ( + desdeo_setup( + desdeo_crossover.SinglePointBinaryCrossover, binary=True + ), + desdeo_run, + ), + }, + ), + "blend_alpha": ( + real, + { + "moocore": (None, moocore_blend_alpha), + "deap": (deap_setup, deap_pairs(deap_tools.cxBlend, 0.5)), + "jMetalPy": ( + jmetalpy_setup, + jmetalpy_pairs(BLXAlphaCrossover(probability=1.0, alpha=0.5)), + ), + "desdeo": ( + desdeo_setup(desdeo_crossover.BlendAlphaCrossover), + desdeo_run, + ), + }, + ), + "single_arithmetic": ( + real, + { + "moocore": ( + None, + lambda x: moocore.single_arithmetic_crossover(x, seed=rng), + ), + "desdeo": ( + desdeo_setup(desdeo_crossover.SingleArithmeticCrossover), + desdeo_run, + ), + }, + ), + "local": ( + real, + { + "moocore": (None, lambda x: moocore.local_crossover(x, seed=rng)), + "desdeo": ( + desdeo_setup(desdeo_crossover.LocalCrossover), + desdeo_run, + ), + }, + ), + "bounded_exponential": ( + real, + { + "moocore": ( + None, + lambda x: moocore.bounded_exponential_crossover( + x, lower=LOWER, upper=UPPER, seed=rng + ), + ), + "desdeo": ( + desdeo_setup(desdeo_crossover.BoundedExponentialCrossover), + desdeo_run, + ), + }, + ), + "de_rand_1_bin": ( + real, + { + "moocore": ( + None, + lambda x: moocore.de_rand_1_bin_crossover( + x, lower=LOWER, upper=UPPER, seed=rng + ), + ), + "pymoo": (pymoo_bounded_setup, pymoo_de_rand_1_bin), + "jMetalPy": (jmetalpy_setup, jmetalpy_de_rand_1_bin), + "desdeo": ( + desdeo_setup(desdeo_crossover.DifferentialEvolutionCrossover), + desdeo_run, + ), + }, + ), + "de_current_1": ( + real, + { + "moocore": ( + None, + lambda x: moocore.de_current_1_crossover( + x, lower=LOWER, upper=UPPER, seed=rng + ), + ), + "pymoo": (pymoo_setup, pymoo_de_current_1), + "jMetalPy": (jmetalpy_setup, jmetalpy_de_current_1), + }, + ), + "parent_centric": ( + real, + { + "moocore": ( + None, + lambda x: moocore.parent_centric_crossover( + x, lower=LOWER, upper=UPPER, seed=rng + ), + ), + "pymoo": (pymoo_setup, pymoo_parent_centric), + "desdeo": ( + desdeo_setup(desdeo_crossover.ParentCentricCrossover), + desdeo_run, + ), + }, + ), +} + + +def get_parents(nparents, nvar, data): + if data == binary: + return rng.integers(0, 2, size=(nparents, nvar)).astype(float) + return rng.uniform(LOWER, UPPER, size=(nparents, nvar)) + + +# Operator names for the plot titles. +titles = { + "uniform": "Uniform", + "single_point": "Single Point", + "blend_alpha": "Blend Alpha", + "single_arithmetic": "Single Arithmetic", + "local": "Local", + "bounded_exponential": "Bounded Exponential", + "de_rand_1_bin": "DE/rand/1/bin", + "de_current_1": "DE/current/1", + "parent_centric": "Parent-Centric", +} + +file_prefix = "crossover" +for op, (data, packages) in operators.items(): + print(f"# Running benchmark: {op}_crossover()") + sweeps = { + "nparents": ( + N_PARENTS, + "Number of parents", + lambda n: get_parents(n, FIXED_NVAR, data), + ), + "nvars": ( + N_VARS, + "Number of variables", + lambda n: get_parents(FIXED_NPARENTS, n, data), + ), + } + for sweep, (n, xlabel, get_testcase) in sweeps.items(): + bench = Bench( + name=f"{op}-{sweep}", + n=n, + bench={k: fun for k, (_, fun) in packages.items()}, + setup={k: setup for k, (setup, _) in packages.items() if setup}, + check=check_offspring, + max_time=10, + ) + bench(get_testcase) + bench.plots( + file_prefix=file_prefix, + title=f"{titles[op]} crossover performance over {xlabel.lower()}", + log="xy", + xlabel=xlabel, + add_name=False, + ) + plt.close("all") + +if "__file__" not in globals(): # Running interactively. + plt.show() diff --git a/python/benchmarks/compare_crossover.py b/python/benchmarks/compare_crossover.py new file mode 100644 index 000000000..bd974c3c0 --- /dev/null +++ b/python/benchmarks/compare_crossover.py @@ -0,0 +1,344 @@ +"""Crossover Return Values +======================= + +This script compares the offspring of the crossover operators in ``moocore`` +with the offspring of the equivalent operators in other packages. + +Each operator crosses over the same parents many times in each package. For +every pair of packages and every decision variable, a two-sample +Kolmogorov-Smirnov test compares their offspring. The script prints, and saves +to ``crossover_ks-.csv``, a matrix with the smallest p-value over the +variables for every pair of packages. The bounds are wide enough that they +never bind, so the comparison does not depend on how each package handles them. + +""" + +import random + +import numpy as np +import polars as pl +from scipy.stats import ks_2samp +import moocore + +from deap import tools as deap_tools +from desdeo.emo.operators import crossover as desdeo_crossover +from desdeo.problem.schema import ( + Objective, + Problem as DesdeoProblem, + Variable, + VariableTypeEnum, +) +from desdeo.tools.patterns import Publisher +from jmetal.core.solution import FloatSolution +from jmetal.operator.crossover import ( + BLXAlphaCrossover, + DifferentialEvolutionCrossover, +) +from pymoo.algorithms.moo.nsde import VariantDE +from pymoo.core.population import Population +from pymoo.core.problem import Problem as PymooProblem +from pymoo.operators.crossover.dex import DEX +from pymoo.operators.crossover.pcx import PCX +from pymoo.operators.crossover.pntx import SinglePointCrossover +from pymoo.operators.crossover.ux import UniformCrossover + +N = 100_000 +LOWER, UPPER = -10.0, 10.0 +# Two parents, and a population whose first row is the target. +A = np.array([0.2, 0.7, 0.95]) +B = np.array([0.6, 0.3, 0.9]) +POPULATION = np.array( + [ + [0.2, 0.7, 0.95], + [0.6, 0.3, 0.9], + [0.5, 0.5, 0.1], + [0.1, 0.9, 0.4], + [0.8, 0.2, 0.6], + [0.4, 0.4, 0.4], + [0.9, 0.1, 0.2], + [0.3, 0.6, 0.8], + ] +) + +rng = np.random.default_rng(42) +random.seed(42) + + +def pair_matrix(a, b, n=N): + return np.tile(np.array([a, b], dtype=float), (n, 1)) + + +def first_offspring(x): + """First offspring of each pair, when rows 2i and 2i+1 are a pair.""" + return np.asarray(x)[0::2] + + +# pymoo +def pymoo_pairs(op, x): + pop = Population.new("X", x) + parents = np.arange(len(x)).reshape(-1, 2) + off = op.do(PymooProblem(n_var=x.shape[1]), pop, parents, random_state=rng) + # pymoo returns all the first offspring before all the second offspring. + return off.get("X")[: len(parents)] + + +def pymoo_with_targets(op, population, others, n=N): + """Offspring of target row 0, with the other parents from `others`.""" + parents = np.column_stack([np.zeros(n, dtype=int), others(n)]) + problem = PymooProblem(n_var=population.shape[1]) + off = op.do( + problem, Population.new("X", population), parents, random_state=rng + ) + return off.get("X") + + +def random_rows(k, exclude_zero=True): + def draw(n): + start = 1 if exclude_zero else 0 + return np.array( + [ + rng.choice(np.arange(start, len(POPULATION)), k, replace=False) + for _ in range(n) + ] + ) + + return draw + + +def pymoo_de_rand_1_bin(population, n=N): + # VariantDE crosses over every row once per call, so call it repeatedly + # and keep the offspring of row 0. + op = VariantDE(variant="DE/rand/1/bin", CR=0.9, F=0.5, gamma=None) + problem = PymooProblem(n_var=population.shape[1]) + pop = Population.new("X", population) + calls = n // 4 + return np.array( + [ + op.do(problem, pop, len(pop), random_state=rng).get("X")[0] + for _ in range(calls) + ] + ) + + +# DEAP +def deap_pairs(fun, x, *args): + ind = x.tolist() + for i in range(0, len(ind) - 1, 2): + ind[i], ind[i + 1] = fun(ind[i], ind[i + 1], *args) + return first_offspring(ind) + + +# jMetalPy +def jmetalpy_solutions(x): + lower, upper = [LOWER] * x.shape[1], [UPPER] * x.shape[1] + solutions = [] + for row in x: + s = FloatSolution(lower, upper, number_of_objectives=1) + s.variables = row.tolist() + solutions.append(s) + return solutions + + +def jmetalpy_pairs(op, x): + s = jmetalpy_solutions(x) + offspring = [] + for i in range(0, len(s) - 1, 2): + offspring.extend(op.execute([s[i], s[i + 1]])) + return first_offspring([o.variables for o in offspring]) + + +def jmetalpy_de(population, current, n=N): + s = jmetalpy_solutions(population) + if current: + op = DifferentialEvolutionCrossover(CR=1.0, F=0.5) + donors = random_rows(2, exclude_zero=False)(n) + else: + op = DifferentialEvolutionCrossover(CR=0.9, F=0.5) + donors = random_rows(3)(n) + op.current_individual = s[0] + offspring = [] + for d in donors: + # execute([a, b, c]) returns c + F (a - b), crossed over with the target. + parents = ( + [s[d[0]], s[d[1]], s[0]] if current else [s[d[1]], s[d[2]], s[d[0]]] + ) + offspring.append(op.execute(parents)[0].variables) + return np.array(offspring) + + +# DESDEO +def desdeo_run(cls, x, targets=None, binary=False, **kwargs): + nvar = x.shape[1] + if binary: + kind, lower, upper = VariableTypeEnum.binary, 0, 1 + else: + kind, lower, upper = VariableTypeEnum.real, LOWER, UPPER + variables = [ + Variable( + name=f"x_{i}", + symbol=f"x_{i}", + variable_type=kind, + lowerbound=lower, + upperbound=upper, + initial_value=lower, + ) + for i in range(1, nvar + 1) + ] + objective = Objective(name="f_1", symbol="f_1", func="x_1", maximize=False) + problem = DesdeoProblem( + name="compare", + description="Crossover comparison", + variables=variables, + objectives=[objective], + ) + op = cls( + problem=problem, publisher=Publisher(), verbosity=0, seed=42, **kwargs + ) + population = pl.DataFrame( + x, schema=[v.symbol for v in variables], orient="row" + ) + if targets is None: + targets = list(range(len(x))) + return op.do(population=population, to_mate=targets).to_numpy() + + +def desdeo_pairs(cls, x, **kwargs): + # DESDEO stacks all the first offspring before all the second offspring. + return desdeo_run(cls, x, **kwargs)[: len(x) // 2] + + +x = pair_matrix(A, B) +xb = pair_matrix(np.zeros(6), np.ones(6)) +targets = np.zeros(N, dtype=int) +bounds = dict(lower=LOWER, upper=UPPER) + +comparisons = { + "uniform": ( + first_offspring(moocore.uniform_crossover(x, seed=rng)), + { + "pymoo": pymoo_pairs(UniformCrossover(prob=1.0), x), + "deap": deap_pairs(deap_tools.cxUniform, x, 0.5), + "desdeo": desdeo_pairs(desdeo_crossover.UniformCrossover, x), + }, + ), + "single_point": ( + first_offspring(moocore.single_point_crossover(xb, seed=rng)), + { + "pymoo": pymoo_pairs(SinglePointCrossover(prob=1.0), xb), + "deap": deap_pairs(deap_tools.cxOnePoint, xb), + "desdeo": desdeo_pairs( + desdeo_crossover.SinglePointBinaryCrossover, xb, binary=True + ), + }, + ), + "blend_alpha": ( + moocore.blend_alpha_crossover(x, seed=rng, **bounds), + { + "deap": deap_pairs(deap_tools.cxBlend, x, 0.5), + "jMetalPy": jmetalpy_pairs( + BLXAlphaCrossover(probability=1.0, alpha=0.5), x + ), + # DESDEO produces two independent offspring per pair. + "desdeo": desdeo_run(desdeo_crossover.BlendAlphaCrossover, x), + }, + ), + "single_arithmetic": ( + first_offspring(moocore.single_arithmetic_crossover(x, seed=rng)), + {"desdeo": desdeo_pairs(desdeo_crossover.SingleArithmeticCrossover, x)}, + ), + "local": ( + first_offspring(moocore.local_crossover(x, seed=rng)), + {"desdeo": desdeo_pairs(desdeo_crossover.LocalCrossover, x)}, + ), + "bounded_exponential": ( + first_offspring( + moocore.bounded_exponential_crossover(x, seed=rng, **bounds) + ), + { + "desdeo": desdeo_pairs( + desdeo_crossover.BoundedExponentialCrossover, x + ) + }, + ), + "de_rand_1_bin": ( + moocore.de_rand_1_bin_crossover( + POPULATION, targets=targets, seed=rng, **bounds + ), + { + "pymoo": pymoo_de_rand_1_bin(POPULATION), + "jMetalPy": jmetalpy_de(POPULATION, current=False), + "desdeo": desdeo_run( + desdeo_crossover.DifferentialEvolutionCrossover, + POPULATION, + targets=targets.tolist(), + ), + }, + ), + "de_current_1": ( + moocore.de_current_1_crossover( + POPULATION, targets=targets, seed=rng, **bounds + ), + { + "pymoo": pymoo_with_targets( + DEX(F=0.5, CR=1.0, at_least_once=False), + POPULATION, + random_rows(2, exclude_zero=False), + ), + "jMetalPy": jmetalpy_de(POPULATION, current=True), + }, + ), + "parent_centric": ( + moocore.parent_centric_crossover( + POPULATION, targets=targets, seed=rng, **bounds + ), + { + "pymoo": pymoo_with_targets( + PCX(eta=0.1, zeta=0.1, prob=1.0), POPULATION, random_rows(2) + ), + "desdeo": desdeo_run( + desdeo_crossover.ParentCentricCrossover, + POPULATION, + targets=targets.tolist(), + ), + }, + ), +} + + +def ks_matrix(samples): + """Smallest KS p-value over the variables, for every pair of packages.""" + names = list(samples) + p = np.full((len(names), len(names)), np.nan) + for i, a in enumerate(names): + for j, b in enumerate(names[i + 1 :], start=i + 1): + x, y = samples[a], samples[b] + p[i, j] = p[j, i] = min( + ks_2samp(x[:, k], y[:, k]).pvalue for k in range(x.shape[1]) + ) + return names, p + + +def fmt(p): + if np.isnan(p): + return "" + return "<0.001" if p < 0.001 else f"{p:.3f}" + + +for op, (ours, others) in comparisons.items(): + samples = {"moocore": np.asarray(ours, dtype=float)} + samples.update({k: np.asarray(v, dtype=float) for k, v in others.items()}) + names, p = ks_matrix(samples) + rows = [ + ["", *names], + *( + [a, *[fmt(p[i, j]) for j in range(len(names))]] + for i, a in enumerate(names) + ), + ] + with open(f"crossover_ks-{op}.csv", "w") as f: + for row in rows: + f.write(",".join(f'"{cell}"' for cell in row) + "\n") + print(f"# {op}") + for row in rows: + print(" | ".join(f"{cell:>8}" for cell in row)) + print() diff --git a/python/doc/source/REFERENCES.bib b/python/doc/source/REFERENCES.bib index b752aab29..1e533f89b 100644 --- a/python/doc/source/REFERENCES.bib +++ b/python/doc/source/REFERENCES.bib @@ -343,6 +343,17 @@ @article{LacKlaFon2017box doi = {10.1016/j.cor.2016.06.021} } +@article{LiZha2009:moead-de, + title = {Multiobjective Optimization Problems with Complicated + {Pareto} sets, {MOEA/D} and {NSGA-II}}, + author = {Li, Hui and Zhang, Qingfu }, + journal = {IEEE Transactions on Evolutionary Computation}, + volume = 13, + number = 2, + pages = {284--302}, + year = 2009 +} + @article{LopVerDreDoe2025, author = { Manuel L{\'o}pez-Ib{\'a}{\~n}ez and Diederick Vermetten and Johann Dreo and Carola Doerr }, title = {Using the Empirical Attainment Function for Analyzing @@ -412,6 +423,19 @@ @article{SchKer2025r2v2 doi = {10.1162/evco.a.366} } +@article{StoPri1997:de, + author = {Storn, Rainer and Price, Kenneth}, + title = {Differential Evolution -- A Simple and Efficient Heuristic + for Global Optimization over Continuous Spaces}, + journal = {Journal of Global Optimization}, + year = 1997, + volume = 11, + number = 4, + pages = {341--359}, + annote = {Proposed differential evolution}, + doi = {10.1023/A:1008202821328} +} + @article{WuAza2001metrics, author = {J. Wu and S. Azam}, title = {Metrics for Quality Assessment of a Multiobjective Design @@ -480,6 +504,29 @@ @incollection{AugBadBroZit2009gecco2 pages = {563--570} } +@misc{BezLopStu2015emoDEsupp, + author = { Leonardo C. T. Bezerra and Manuel L{\'o}pez-Ib{\'a}{\~n}ez and Thomas St{\"u}tzle }, + title = {To {DE} or Not to {DE}? {Multi}-objective Differential + Evolution Revisited from a Component-Wise Perspective: {Supplementary} material}, + howpublished = {\url{https://iridia.ulb.ac.be/supp/IridiaSupp2015-001/}}, + year = 2015 +} + +@incollection{BezLopStu2015emode, + booktitle = { Evolutionary Multi-criterion Optimization, EMO 2015 Part {I}}, + address = { Heidelberg, Germany}, + series = {Lecture Notes in Computer Science}, + volume = 9018, + year = 2015, + publisher = {Springer}, + editor = { Ant{\'o}nio Gaspar{-}Cunha and Carlos Henggeler Antunes and Carlos A. {Coello Coello} }, + author = { Leonardo C. T. Bezerra and Manuel L{\'o}pez-Ib{\'a}{\~n}ez and Thomas St{\"u}tzle }, + title = {To {DE} or Not to {DE}? {Multi}-objective Differential + Evolution Revisited from a Component-Wise Perspective}, + pages = {48--63}, + doi = {10.1007/978-3-319-15934-8_4} +} + @incollection{BezLopStu2017emo, address = { Cham, Switzerland}, series = {Lecture Notes in Computer Science}, @@ -550,6 +597,14 @@ @incollection{CoeSie2004igd keywords = {IGD} } +@book{EibSmi2003, + title = {Introduction to Evolutionary Computing}, + author = { Agoston E. Eiben and Smith, James E. }, + publisher = {Springer}, + year = 2003, + isbn = 3540401849 +} + @incollection{EmmFon2011emo, booktitle = { Evolutionary Multi-criterion Optimization, EMO 2011}, address = {Berlin~/ Heidelberg}, @@ -590,6 +645,17 @@ @incollection{EmmFon2011emo pages = {121--135} } +@incollection{EshSch1992, + isbn = {1-55860-263-1}, + year = 1993, + publisher = {Morgan Kaufmann Publishers}, + booktitle = {Foundations of Genetic Algorithms (FOGA)}, + editor = { Darrell Whitley }, + author = { Larry J. Eshelman and J. David Schaffer }, + title = {Real-Coded Genetic Algorithms and Interval-Schemata}, + pages = {187--202} +} + @book{FanWan1994numtheory, author = {Fang, K. T. and Wang, Y.}, title = {Number-Theoretic Methods in Statistics}, @@ -681,6 +747,15 @@ @incollection{Gla2017fast keywords = {archiving} } +@book{Goldberg89, + author = { David E. Goldberg }, + title = {Genetic Algorithms in Search, Optimization and + Machine Learning}, + publisher = {Addison-Wesley}, + address = { Boston, MA}, + year = 1989 +} + @incollection{GruFon2009:emaa, editor = { Thomas Bartz-Beielstein and Marco Chiarandini and Lu{\'i}s Paquete and Mike Preuss }, year = 2010, @@ -744,6 +819,13 @@ @techreport{HanJas1998 annote = {Proposed R2 indicator} } +@book{Hol75, + author = { J. Holland }, + title = {Adaptation in Natural and Artificial Systems}, + publisher = {University of Michigan Press}, + year = 1975 +} + @inproceedings{IshHeSha2019regular, month = jun, year = 2019, @@ -786,6 +868,19 @@ @inproceedings{KnoCor2002cec pages = {711--716} } +@inproceedings{KukLam2005:gde3, + address = {Piscataway, NJ}, + publisher = {IEEE Press}, + month = sep, + year = 2005, + booktitle = {Proceedings of the 2005 Congress on Evolutionary Computation (CEC 2005)}, + key = {IEEE CEC}, + author = {Kukkonen, S. and Lampinen, J.}, + title = {{GDE3}: the third evolution step of generalized differential + evolution}, + pages = {443--450} +} + @incollection{Lop2026hvapprox, location = {San Jos{\'e}, Costa Rica}, address = { New York, NY}, @@ -881,6 +976,19 @@ @misc{Rob2018unreasonable annote = {Proposed the Rphi low-discrepancy sequence} } +@incollection{RobFil05demo, + address = {Berlin~/ Heidelberg}, + publisher = {Springer}, + volume = 3410, + series = {Lecture Notes in Computer Science}, + editor = { Carlos A. {Coello Coello} and Hern{\'a}ndez Aguirre, Arturo and Eckart Zitzler }, + year = 2005, + booktitle = { Evolutionary Multi-criterion Optimization, EMO 2005}, + author = {Tea Robi{\v c} and Bogdan Filipi{\v c}}, + title = {{DEMO}: Differential Evolution for Multiobjective Optimization}, + pages = {520--533} +} + @book{RubMel1998simulation, author = {R. Y. Rubinstein and B. Melamed}, title = {Modern simulation and modeling}, @@ -890,6 +998,17 @@ @book{RubMel1998simulation notes = {Uniform sampling from the simplex.} } +@inproceedings{Syswerda89, + publisher = {Morgan Kaufmann Publishers, San Mateo, CA}, + editor = { J. David Schaffer }, + year = 1989, + booktitle = {Proceedings of the Third International Conference on Genetic Algorithms (ICGA'89)}, + author = { Gilbert Syswerda }, + title = {Uniform Crossover in Genetic Algorithms}, + pages = {2--9}, + keywords = {uniform crossover} +} + @inproceedings{VelLam1998gp, year = 1998, publisher = {Stanford University Bookstore}, diff --git a/python/doc/source/_static/bench/crossover_bench-blend_alpha-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-blend_alpha-nparents-time.png new file mode 100644 index 000000000..b95792893 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-blend_alpha-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-blend_alpha-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-blend_alpha-nvars-time.png new file mode 100644 index 000000000..a416e08f4 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-blend_alpha-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-bounded_exponential-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-bounded_exponential-nparents-time.png new file mode 100644 index 000000000..4eda58f0b Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-bounded_exponential-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-bounded_exponential-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-bounded_exponential-nvars-time.png new file mode 100644 index 000000000..9ce89b842 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-bounded_exponential-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-de_current_1-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-de_current_1-nparents-time.png new file mode 100644 index 000000000..77dda370e Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-de_current_1-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-de_current_1-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-de_current_1-nvars-time.png new file mode 100644 index 000000000..50ca05275 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-de_current_1-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-de_rand_1_bin-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-de_rand_1_bin-nparents-time.png new file mode 100644 index 000000000..5b5ac9a0f Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-de_rand_1_bin-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-de_rand_1_bin-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-de_rand_1_bin-nvars-time.png new file mode 100644 index 000000000..587f5e334 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-de_rand_1_bin-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-local-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-local-nparents-time.png new file mode 100644 index 000000000..3a438da00 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-local-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-local-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-local-nvars-time.png new file mode 100644 index 000000000..6131d5848 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-local-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-parent_centric-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-parent_centric-nparents-time.png new file mode 100644 index 000000000..e1663bdc6 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-parent_centric-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-parent_centric-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-parent_centric-nvars-time.png new file mode 100644 index 000000000..40451b363 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-parent_centric-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-single_arithmetic-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-single_arithmetic-nparents-time.png new file mode 100644 index 000000000..f999bdcb6 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-single_arithmetic-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-single_arithmetic-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-single_arithmetic-nvars-time.png new file mode 100644 index 000000000..5cc4af9b2 Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-single_arithmetic-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-single_point-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-single_point-nparents-time.png new file mode 100644 index 000000000..d2bbcd03d Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-single_point-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-single_point-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-single_point-nvars-time.png new file mode 100644 index 000000000..1179f2d4e Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-single_point-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-uniform-nparents-time.png b/python/doc/source/_static/bench/crossover_bench-uniform-nparents-time.png new file mode 100644 index 000000000..35321e75c Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-uniform-nparents-time.png differ diff --git a/python/doc/source/_static/bench/crossover_bench-uniform-nvars-time.png b/python/doc/source/_static/bench/crossover_bench-uniform-nvars-time.png new file mode 100644 index 000000000..3d938391b Binary files /dev/null and b/python/doc/source/_static/bench/crossover_bench-uniform-nvars-time.png differ diff --git a/python/doc/source/_static/bench/crossover_ks-blend_alpha.csv b/python/doc/source/_static/bench/crossover_ks-blend_alpha.csv new file mode 100644 index 000000000..72c43fdd4 --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-blend_alpha.csv @@ -0,0 +1,5 @@ +"","moocore","deap","jMetalPy","desdeo" +"moocore","","0.517","0.064","0.513" +"deap","0.517","","0.118","0.814" +"jMetalPy","0.064","0.118","","0.014" +"desdeo","0.513","0.814","0.014","" diff --git a/python/doc/source/_static/bench/crossover_ks-bounded_exponential.csv b/python/doc/source/_static/bench/crossover_ks-bounded_exponential.csv new file mode 100644 index 000000000..2f174e87c --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-bounded_exponential.csv @@ -0,0 +1,3 @@ +"","moocore","desdeo" +"moocore","","0.199" +"desdeo","0.199","" diff --git a/python/doc/source/_static/bench/crossover_ks-de_current_1.csv b/python/doc/source/_static/bench/crossover_ks-de_current_1.csv new file mode 100644 index 000000000..a5291d899 --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-de_current_1.csv @@ -0,0 +1,4 @@ +"","moocore","pymoo","jMetalPy" +"moocore","","<0.001","0.905" +"pymoo","<0.001","","<0.001" +"jMetalPy","0.905","<0.001","" diff --git a/python/doc/source/_static/bench/crossover_ks-de_rand_1_bin.csv b/python/doc/source/_static/bench/crossover_ks-de_rand_1_bin.csv new file mode 100644 index 000000000..6fd181c10 --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-de_rand_1_bin.csv @@ -0,0 +1,5 @@ +"","moocore","pymoo","jMetalPy","desdeo" +"moocore","","<0.001","0.275","0.786" +"pymoo","<0.001","","<0.001","<0.001" +"jMetalPy","0.275","<0.001","","0.549" +"desdeo","0.786","<0.001","0.549","" diff --git a/python/doc/source/_static/bench/crossover_ks-local.csv b/python/doc/source/_static/bench/crossover_ks-local.csv new file mode 100644 index 000000000..dadb780c2 --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-local.csv @@ -0,0 +1,3 @@ +"","moocore","desdeo" +"moocore","","0.496" +"desdeo","0.496","" diff --git a/python/doc/source/_static/bench/crossover_ks-parent_centric.csv b/python/doc/source/_static/bench/crossover_ks-parent_centric.csv new file mode 100644 index 000000000..b859f4f5e --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-parent_centric.csv @@ -0,0 +1,4 @@ +"","moocore","pymoo","desdeo" +"moocore","","0.233","0.736" +"pymoo","0.233","","0.343" +"desdeo","0.736","0.343","" diff --git a/python/doc/source/_static/bench/crossover_ks-single_arithmetic.csv b/python/doc/source/_static/bench/crossover_ks-single_arithmetic.csv new file mode 100644 index 000000000..1081baac3 --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-single_arithmetic.csv @@ -0,0 +1,3 @@ +"","moocore","desdeo" +"moocore","","0.999" +"desdeo","0.999","" diff --git a/python/doc/source/_static/bench/crossover_ks-single_point.csv b/python/doc/source/_static/bench/crossover_ks-single_point.csv new file mode 100644 index 000000000..4b6156d9d --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-single_point.csv @@ -0,0 +1,5 @@ +"","moocore","pymoo","deap","desdeo" +"moocore","","0.995","0.996","0.650" +"pymoo","0.995","","0.531","0.422" +"deap","0.996","0.531","","0.177" +"desdeo","0.650","0.422","0.177","" diff --git a/python/doc/source/_static/bench/crossover_ks-uniform.csv b/python/doc/source/_static/bench/crossover_ks-uniform.csv new file mode 100644 index 000000000..0cf2da533 --- /dev/null +++ b/python/doc/source/_static/bench/crossover_ks-uniform.csv @@ -0,0 +1,5 @@ +"","moocore","pymoo","deap","desdeo" +"moocore","","0.249","0.765","1.000" +"pymoo","0.249","","0.658","0.699" +"deap","0.765","0.658","","0.987" +"desdeo","1.000","0.699","0.987","" diff --git a/python/doc/source/benchmarks-links.inc.rst b/python/doc/source/benchmarks-links.inc.rst new file mode 100644 index 000000000..2b03a5714 --- /dev/null +++ b/python/doc/source/benchmarks-links.inc.rst @@ -0,0 +1,22 @@ +.. _BoFire: https://experimental-design.github.io/bofire/ +.. _BoTorch: https://botorch.org/ +.. _DEAP: https://deap.readthedocs.io/en/master/ +.. _DESDEO: https://desdeo.readthedocs.io/en/latest/ +.. _Nevergrad: https://facebookresearch.github.io/nevergrad/ +.. _Numba: https://numba.pydata.org/numba-doc/dev/index.html +.. _Optuna: https://optuna.org/ +.. _Xopt: https://xopt.xopt.org/index.html +.. _fast-pareto: https://github.com/nabenabe0928/fast-pareto +.. _jMetalPy: https://jmetal.github.io/jMetalPy/index.html +.. _dmosopt: https://dmosopt.github.io/dmosopt/ +.. _moarchiving: https://cma-es.github.io/moarchiving/ +.. _moocore: https://multi-objective.github.io/moocore/python/ +.. _paretobench: https://github.com/electronsandstuff/ParetoBench +.. _paretoset: https://github.com/tommyod/paretoset +.. _patatune: https://cms-patatrack.github.io/patatune/ +.. _platypus: https://github.com/Project-Platypus/Platypus +.. _pygmo: https://esa.github.io/pygmo2/ +.. _pymoo: https://pymoo.org/ +.. _pysamoo: https://anyoptimization.com/projects/pysamoo/ +.. _seqme: https://seqme.readthedocs.io +.. _trieste: https://secondmind-labs.github.io/trieste diff --git a/python/doc/source/benchmarks.inc.rst b/python/doc/source/benchmarks.inc.rst deleted file mode 100644 index 951d692b7..000000000 --- a/python/doc/source/benchmarks.inc.rst +++ /dev/null @@ -1,285 +0,0 @@ -Benchmarks -========== - -The following plots compare the performance of `moocore`_, `pymoo`_, `BoTorch`_, `DESDEO`_, `paretoset`_, `Nevergrad`_, `jMetalPy`_, `fast-pareto`_, `Optuna`_, `moarchiving`_, and `patatune`_. - -Other Python packages are not included in the comparison because they are based -on these packages for the functionality benchmarked, so they are **at least as -slow** as them. For example `Xopt`_ and `BoFire`_ use `BoTorch`_, `pysamoo`_ is -an extension of `pymoo`_. We do not include functions from the Bayesian -optimization toolbox `trieste`_, because it is much slower than `BoTorch`_ and -too slow to run the benchmarks in a reasonable time -(https://github.com/secondmind-labs/trieste/issues/917). Similarly, we do not -include `paretobench`_ because it is always slower than `BoTorch`_ for -identifying nondominated points and more than 2000 times slower than `moocore`_ -for computing the hypervolume in 3D -(https://github.com/electronsandstuff/ParetoBench/issues/57), which is too slow -to include in the benchmarks. - -Not all packages provide the same functionality. For example, `pymoo`_ does not -provide the :ref:`epsilon indicator `. `BoTorch`_ and -`paretobench`_ only provide the hypervolume. `paretoset`_ and `fast-pareto`_ -only identify nondominated points. `DESDEO`_, `DEAP`_, and `jMetalPy`_ use -`moocore`_ for computing the hypervolume, but other functionality, such as -filtering dominated points, is still slower than `moocore`_. Recent versions -of `seqme`_ and `pymoo`_ (≥0.6.2) already use `moocore`_ for most functionality -benchmarked here. - -We would like to benchmark `pygmo`_, however, it is currently impossible to install using ``pip`` (See https://github.com/esa/pygmo2/issues/152). - -We would also like to benchmark `platypus`_, but its API is too different from -other packages. For example, it only supports filtering dominated points by -creating a new list of points, but it cannot identify their position in the -input. It cannot calculate hypervolume without normalization. It doesn't -support :obj:`~numpy.typing.ArrayLike` inputs. - -The source code for the benchmarks below can be found at https://github.com/multi-objective/moocore/tree/main/python/benchmarks . - -.. _bench-ndom: - -Identifying nondominated points -------------------------------- - -The following plots compare the speed of finding nondominated solutions, equivalent to :func:`moocore.is_nondominated`. We test both ``keep_weakly=False`` (the default) and ``keep_weakly=True``. The former is only supported by `BoTorch`_ and `paretoset`_. Thus, for the purpose of benchmarking the other packages with ``keep_weakly=False``, we remove duplicated points from the data and exclude this time from the results. This makes all packages, except `moocore`_, `BoTorch`_ and `paretoset`_, seem faster than they actually are, because the user is forced to remove duplicated points when using them. - -For 2D and 3D, the plots show that `moocore`_ is 10 times faster than `DESDEO`_ and `paretoset`_ and 100-1000 times faster than the other packages. `fast-pareto`_ claims to implement a :math:`O(n\log n)` algorithm for 3D, but the benchmarks below indicate a :math:`O(n^2)` complexity, similar to other packages and significantly slower than `moocore`_. - -|ndom_bench-test2D-200k| |ndom_bench-ran3d-40k| - -.. |ndom_bench-test2D-200k| image:: _static/bench/ndom_bench-test2D-200k-time.png - :width: 49% - -.. |ndom_bench-ran3d-40k| image:: _static/bench/ndom_bench-ran3d-40k-time.png - :width: 49% - -|wndom_bench-test2D-200k| |wndom_bench-ran3d-40k| - -.. |wndom_bench-test2D-200k| image:: _static/bench/wndom_bench-test2D-200k-time.png - :width: 49% - -.. |wndom_bench-ran3d-40k| image:: _static/bench/wndom_bench-ran3d-40k-time.png - :width: 49% - - -For dimensions larger than 3, :func:`moocore.is_nondominated` uses the -best-known :math:`O(n\log^{m-2} n)` algorithm. `fast-pareto`_ claims to be -using a :math:`O(n\log^{m-2} n)` algorithm, but the benchmarks suggest that it -uses a :math:`O(m n^2)` algorithm, although slower than most other packages and -almost 100 times slower than `moocore`_. The following plots with -``keep_weakly=True`` show that `moocore`_ is at least 10 times faster than the -other packages, even for 10D. - -|wndom_bench-sphere-4d| |wndom_bench-convex-4d| - -|wndom_bench-sphere-5d| |wndom_bench-rmnk-10d| - -.. |wndom_bench-sphere-4d| image:: _static/bench/wndom_bench-sphere-4d-time.png - :width: 49% - -.. |wndom_bench-convex-4d| image:: _static/bench/wndom_bench-convex-4d-time.png - :width: 49% - -.. |wndom_bench-sphere-5d| image:: _static/bench/wndom_bench-sphere-5d-time.png - :width: 49% - -.. |wndom_bench-rmnk-10d| image:: _static/bench/wndom_bench-rmnk-10d-time.png - :width: 49% - - -Similar results are obtained with ``keep_weakly=False``, however, this setting -is only supported by `moocore`_, `BoTorch`_ and `paretoset`_, thus the other -packages are actually slower to use than reported here, because the user is -forced to remove duplicated points when using them. - -|ndom_bench-sphere-4d| |ndom_bench-convex-4d| - -|ndom_bench-sphere-5d| |ndom_bench-rmnk-10d| - -.. |ndom_bench-sphere-4d| image:: _static/bench/ndom_bench-sphere-4d-time.png - :width: 49% - -.. |ndom_bench-convex-4d| image:: _static/bench/ndom_bench-convex-4d-time.png - :width: 49% - -.. |ndom_bench-sphere-5d| image:: _static/bench/ndom_bench-sphere-5d-time.png - :width: 49% - -.. |ndom_bench-rmnk-10d| image:: _static/bench/ndom_bench-rmnk-10d-time.png - :width: 49% - - -.. _bench-ndsort: - -Nondominated Sorting (Pareto ranking) -------------------------------------- - -Nondominated sorting (or Pareto ranking) assigns a rank to each point according -to Pareto-optimality. It can be used to split a set of points into layers of -nondominated fronts. For more details, see :func:`moocore.pareto_rank`. This -function in `moocore`_ uses a :math:`O(n \log n)` algorithm for 2 dimensions -and a :math:`O(n^2\log n)` algorithm for 3 dimensions. As the plots show, this -makes `moocore`_ significantly faster than `paretoset`_, `pymoo`_ and -`DESDEO`_. With higher dimensions, `moocore`_ uses a :math:`O(n^2\log^{m-2}n)` -algorithm instead of the :math:`O(n^3)` algorithm in `paretoset`_ and `DESDEO`_ -(latter is compiled with `Numba`_). In these benchmarks points are sampled -uniformly at random within the hypercube. - -|ndsort_bench-ran-2d| |ndsort_bench-ran-3d| - -|ndsort_bench-ran-4d| |ndsort_bench-ran-5d| - -|ndsort_bench-ran-9d| |ndsort_bench-ran-10d| - -.. |ndsort_bench-ran-2d| image:: _static/bench/ndsort_bench-ran-2d-time.png - :width: 49% - -.. |ndsort_bench-ran-3d| image:: _static/bench/ndsort_bench-ran-3d-time.png - :width: 49% - -.. |ndsort_bench-ran-4d| image:: _static/bench/ndsort_bench-ran-4d-time.png - :width: 49% - -.. |ndsort_bench-ran-5d| image:: _static/bench/ndsort_bench-ran-5d-time.png - :width: 49% - -.. |ndsort_bench-ran-9d| image:: _static/bench/ndsort_bench-ran-9d-time.png - :width: 49% - -.. |ndsort_bench-ran-10d| image:: _static/bench/ndsort_bench-ran-10d-time.png - :width: 49% - - -.. _bench-hv: - -Exact computation of hypervolume --------------------------------- - -The following plots compare the speed of computing the :ref:`hypervolume indicator ` in 3D, 4D, 5D and 6D. As the plots show, `moocore`_ is 100 times faster than the other packages and 1000 times faster than `BoTorch`_ and, by extension, `Xopt`_ and `BoFire`_. - -For 3D, `Optuna`_ uses a :math:`O(n^2)` algorithm, while `moocore`_ uses the best-known :math:`O(n\log n)` one, so the gap between the two will get significantly larger with a larger number of points than those tested here. `moarchiving`_ uses the same algorithms as `moocore`_ for 3D and 4D, but a pure Python implementation, so it is still 10-100 times slower than `moocore`_'s C implementation. `BoTorch`_ and `fast-pareto`_ are not included for more than 3 objectives because **they are tens of thousands of times slower** than `moocore`_. - -Recent versions of `DESDEO`_, `DEAP`_ (≥1.4.4), `pymoo`_ (≥0.6.1.6) and `jMetalPy`_ (≥1.9.0) already use `moocore`_ for hypervolume computation. -The implementation in `patatune`_ appears to be wrong (https://github.com/cms-patatrack/patatune/issues/52). - -|hv_bench-DTLZLinearShape-3d| |hv_bench-DTLZLinearShape-4d| - -|hv_bench-DTLZLinearShape-5d| |hv_bench-DTLZLinearShape-6d| - -.. |hv_bench-DTLZLinearShape-3d| image:: _static/bench/hv_bench-DTLZLinearShape.3d-time.png - :width: 49% - -.. |hv_bench-DTLZLinearShape-4d| image:: _static/bench/hv_bench-DTLZLinearShape.4d-time.png - :width: 49% - -.. |hv_bench-DTLZLinearShape-5d| image:: _static/bench/hv_bench-DTLZLinearShape.5d-time.png - :width: 49% - -.. |hv_bench-DTLZLinearShape-6d| image:: _static/bench/hv_bench-DTLZLinearShape.6d-time.png - :width: 49% - -.. _bench-hvapprox: - -Approximation of the hypervolume --------------------------------- - -The following plots compare the accuracy and speed of approximating the -hypervolume with the various methods provided by :func:`moocore.hv_approx`. The -plots show that there is no clear winner, in terms of approximation error, -between methods ``Rphi-FWE+`` (default) and ``DZ2019-HW``, but both produce -consistently lower approximation errors than method ``DZ2019-MC``. However, ``Rphi-FWE+`` is as fast as ``DZ2019-MC`` and both are consistently faster than ``DZ2019-HW``, in particular with higher number of objectives. - -The exact computation of the hypervolume for less than 5D takes milliseconds -for thousands of points, whereas approximating the hypervolume is significantly -slower and, thus, not worth doing. - -Approximating the hypervolume becomes more useful for dimensions higher than 5, -where the exact computation becomes noticeably slower with hundreds of points. -For such problems, method ``DZ2019-HW`` becomes significantly slower than -``Rphi-FWE+``. - -|hvapprox_bench-DTLZLinearShape-6d-values| |hvapprox_bench-DTLZLinearShape-6d-time| - -|hvapprox_bench-DTLZSphereShape-6d-values| |hvapprox_bench-DTLZSphereShape-6d-time| - -|hvapprox_bench-ran-6d-values| |hvapprox_bench-ran-6d-time| - -|hvapprox_bench-DTLZLinearShape-9d-values| |hvapprox_bench-DTLZLinearShape-9d-time| - -|hvapprox_bench-ran-9d-values| |hvapprox_bench-ran-9d-time| - -|hvapprox_bench-DTLZSphereShape-10d-values| |hvapprox_bench-DTLZSphereShape-10d-time| - -.. |hvapprox_bench-DTLZSphereShape-6d-values| image:: _static/bench/hvapprox_bench-DTLZSphereShape.6d-values.png - :width: 49% - -.. |hvapprox_bench-DTLZSphereShape-6d-time| image:: _static/bench/hvapprox_bench-DTLZSphereShape.6d-time.png - :width: 49% - -.. |hvapprox_bench-DTLZLinearShape-6d-values| image:: _static/bench/hvapprox_bench-DTLZLinearShape.6d-values.png - :width: 49% - -.. |hvapprox_bench-DTLZLinearShape-6d-time| image:: _static/bench/hvapprox_bench-DTLZLinearShape.6d-time.png - :width: 49% - -.. |hvapprox_bench-ran-6d-values| image:: _static/bench/hvapprox_bench-ran.6d-values.png - :width: 49% - -.. |hvapprox_bench-ran-6d-time| image:: _static/bench/hvapprox_bench-ran.6d-time.png - :width: 49% - -.. |hvapprox_bench-DTLZLinearShape-9d-values| image:: _static/bench/hvapprox_bench-DTLZLinearShape.9d-values.png - :width: 49% - -.. |hvapprox_bench-DTLZLinearShape-9d-time| image:: _static/bench/hvapprox_bench-DTLZLinearShape.9d-time.png - :width: 49% - -.. |hvapprox_bench-ran-9d-values| image:: _static/bench/hvapprox_bench-ran.9d-values.png - :width: 49% - -.. |hvapprox_bench-ran-9d-time| image:: _static/bench/hvapprox_bench-ran.9d-time.png - :width: 49% - -.. |hvapprox_bench-DTLZSphereShape-10d-values| image:: _static/bench/hvapprox_bench-DTLZSphereShape.10d-values.png - :width: 49% - -.. |hvapprox_bench-DTLZSphereShape-10d-time| image:: _static/bench/hvapprox_bench-DTLZSphereShape.10d-time.png - :width: 49% - - -Epsilon and IGD+ indicators ---------------------------- - -The following plots benchmark the speed of computing the :ref:`epsilon indicator ` metric and :ref:`IGD+ indicator `. Although the algorithms for computing these metrics are relatively simple and easy to vectorize in Python, the `moocore`_ implementation is still 10 to 100 times faster. - -|fig_eps_bench| |fig_igd_plus_bench| - -.. |fig_eps_bench| image:: _static/bench/eps_bench-rmnk_10D_random_search-time.png - :width: 48% - -.. |fig_igd_plus_bench| image:: _static/bench/igd_plus_bench-ran.1000pts.3d-time.png - :width: 48% - - - -.. _BoFire: https://experimental-design.github.io/bofire/ -.. _BoTorch: https://botorch.org/ -.. _DEAP: https://deap.readthedocs.io/en/master/ -.. _DESDEO: https://desdeo.readthedocs.io/en/latest/ -.. _Nevergrad: https://facebookresearch.github.io/nevergrad/ -.. _Numba: https://numba.pydata.org/numba-doc/dev/index.html -.. _Optuna: https://optuna.org/ -.. _Xopt: https://xopt.xopt.org/index.html -.. _fast-pareto: https://github.com/nabenabe0928/fast-pareto -.. _jMetalPy: https://jmetal.github.io/jMetalPy/index.html -.. _dmosopt: https://dmosopt.github.io/dmosopt/ -.. _moarchiving: https://cma-es.github.io/moarchiving/ -.. _moocore: https://multi-objective.github.io/moocore/python/ -.. _paretobench: https://github.com/electronsandstuff/ParetoBench -.. _paretoset: https://github.com/tommyod/paretoset -.. _patatune: https://cms-patatrack.github.io/patatune/ -.. _platypus: https://github.com/Project-Platypus/Platypus -.. _pygmo: https://esa.github.io/pygmo2/ -.. _pymoo: https://pymoo.org/ -.. _pysamoo: https://anyoptimization.com/projects/pysamoo/ -.. _seqme: https://seqme.readthedocs.io -.. _trieste: https://secondmind-labs.github.io/trieste diff --git a/python/doc/source/benchmarks/crossover/blend_alpha.rst b/python/doc/source/benchmarks/crossover/blend_alpha.rst new file mode 100644 index 000000000..4ec10ea9e --- /dev/null +++ b/python/doc/source/benchmarks/crossover/blend_alpha.rst @@ -0,0 +1,52 @@ +.. _bench-crossover-blend-alpha: + +Blend crossover (BLX-α) +======================= + +:func:`~moocore.blend_alpha_crossover` samples each variable of the offspring +uniformly from the interval between the parents, widened on each side by +``alpha`` times its length. `DEAP`_ (``cxBlend``), `jMetalPy`_ +(``BLXAlphaCrossover``) and `DESDEO`_ (``BlendAlphaCrossover``) provide it. + +`moocore`_ produces one offspring per pair of parents, and the other packages +produce two. The two offspring of `DEAP`_ mirror each other, so the KS tests +use the first offspring of each pair from `DEAP`_ and `jMetalPy`_, and every +offspring from `DESDEO`_, whose offspring are independent. The timings call +`moocore`_ twice, so that every package produces the same number of offspring. + +The packages handle bounds differently. `moocore`_ and `DESDEO`_ clip the +sampling interval to the bounds, `jMetalPy`_ clips each sampled value, and +`DEAP`_ ignores the bounds. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-blend_alpha.csv + :header-rows: 1 + :stub-columns: 1 + +All implementations produce the same distribution of offspring. The p-value of +0.014 between `jMetalPy`_ and `DESDEO`_ is one of the low values expected by +chance: with two other random seeds, it was 0.07 and 0.31. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-blend_alpha-nparents-time.png + :width: 80% + +With 10 parents, `DEAP`_ is faster than the two calls to `moocore`_. From 20 +parents on, `moocore`_ is the fastest. With 10 000 parents, it takes 2.1 ms, +while `DEAP`_, `DESDEO`_ and `jMetalPy`_ take 7.8, 15 and 51 times as long. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-blend_alpha-nvars-time.png + :width: 80% + +With 2 variables, `DEAP`_ is faster than `moocore`_. With more variables, +`moocore`_ is the fastest. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/bounded_exponential.rst b/python/doc/source/benchmarks/crossover/bounded_exponential.rst new file mode 100644 index 000000000..cfc82f733 --- /dev/null +++ b/python/doc/source/benchmarks/crossover/bounded_exponential.rst @@ -0,0 +1,37 @@ +.. _bench-crossover-bounded-exponential: + +Bounded exponential crossover (BEX) +=================================== + +:func:`~moocore.bounded_exponential_crossover` moves each offspring away from +its parent by a random step that follows an exponential distribution truncated +at the bounds. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-bounded_exponential.csv + :header-rows: 1 + :stub-columns: 1 + +Both implementations produce the same distribution of offspring. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-bounded_exponential-nparents-time.png + :width: 80% + +`DESDEO`_ takes 20 times as long as `moocore`_ with 10 parents, and 1.2 times +as long with 10 000 parents. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-bounded_exponential-nvars-time.png + :width: 80% + +`DESDEO`_ takes 3.7 to 15 times as long as `moocore`_. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/de_current_1.rst b/python/doc/source/benchmarks/crossover/de_current_1.rst new file mode 100644 index 000000000..fd5afc62d --- /dev/null +++ b/python/doc/source/benchmarks/crossover/de_current_1.rst @@ -0,0 +1,47 @@ +.. _bench-crossover-de-current-1: + +Differential evolution: DE/current/1 +==================================== + +:func:`~moocore.de_current_1_crossover` adds a scaled difference of two random +members of the population to the target, as in MOEA/D-DE. `pymoo`_ provides +it as ``DEX`` with ``at_least_once=False``, when the target is given as the +first parent. `jMetalPy`_ provides it through the order of the parents of +``DifferentialEvolutionCrossover``, as its MOEA/D does. `jMetalPy`_ always +takes one random variable from the mutant, so it equals this operator only +with :math:`CR = 1`, which is the default. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-de_current_1.csv + :header-rows: 1 + :stub-columns: 1 + +`moocore`_ and `jMetalPy`_ produce the same distribution of offspring. +`pymoo`_ multiplies the scaling factor by a random factor within :math:`1 \pm +5 \cdot 10^{-5}`, so its offspring are within :math:`2 \cdot 10^{-5}` of the +values that the other packages produce, but never equal to them, which is enough +difference to be detected by the KS test. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-de_current_1-nparents-time.png + :width: 80% + +`moocore`_ is the fastest for every number of parents. With 10 000 parents, +it takes 1.2 ms, while `pymoo`_ and `jMetalPy`_ take 33 and 230 times as +long. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-de_current_1-nvars-time.png + :width: 80% + +`moocore`_ is the fastest for every number of variables. With 1000 variables, +`pymoo`_ and `jMetalPy`_ take 2.7 and 950 times as long. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/de_rand_1_bin.rst b/python/doc/source/benchmarks/crossover/de_rand_1_bin.rst new file mode 100644 index 000000000..20b9f8bc2 --- /dev/null +++ b/python/doc/source/benchmarks/crossover/de_rand_1_bin.rst @@ -0,0 +1,48 @@ +.. _bench-crossover-de-rand-1-bin: + +Differential evolution: DE/rand/1/bin +===================================== + +:func:`~moocore.de_rand_1_bin_crossover` adds a scaled difference of two random +members of the population to a third one, and takes each variable of the +offspring from this mutant or from the target. One random variable always +comes from the mutant. `pymoo`_ (``VariantDE``), `jMetalPy`_ +(``DifferentialEvolutionCrossover``) and `DESDEO`_ +(``DifferentialEvolutionCrossover``) provide this operator. `pymoo`_ forces a +variable to come from the mutant only when no variable was chosen, instead of +forcing a random variable in every offspring. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-de_rand_1_bin.csv + :header-rows: 1 + :stub-columns: 1 + +`moocore`_, `jMetalPy`_ and `DESDEO`_ produce the same distribution of +offspring. `pymoo`_ differs from all three because of its rule for the +forced variable. With 3 variables and :math:`CR = 0.9`, it keeps the value of +the target in 13% to 14% of the variables, against 9% to 11% in the other +packages. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-de_rand_1_bin-nparents-time.png + :width: 80% + +`moocore`_ is the fastest for every number of parents. With 10 000 parents, +it takes 1.9 ms, while `DESDEO`_, `pymoo`_ and `jMetalPy`_ take 2.2, 17 and 140 +times as long. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-de_rand_1_bin-nvars-time.png + :width: 80% + +`moocore`_ is the fastest for every number of variables. With 1000 variables, +`pymoo`_, `DESDEO`_ and `jMetalPy`_ take 2.6, 6.5 and 550 times as long. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/index.rst b/python/doc/source/benchmarks/crossover/index.rst new file mode 100644 index 000000000..7b1f9a31d --- /dev/null +++ b/python/doc/source/benchmarks/crossover/index.rst @@ -0,0 +1,38 @@ +.. _bench-crossover: + +Crossover operators +=================== + +These pages compare the :ref:`crossover operators ` of +`moocore`_ with the equivalent operators of `pymoo`_, `DEAP`_, `jMetalPy`_ and +`DESDEO`_, when they exist. Each page compares the offspring that the packages +produce and the time they take. + +The KS test results show whether two packages produce offspring with the same +distribution. Each package crosses over the same parents 100 000 times, with +bounds wide enough that they never bind, so the results do not depend on how +each package handles bounds. For every pair of packages and every decision +variable, a two-sample Kolmogorov-Smirnov test compares their offspring. Each +table shows the smallest p-value over the variables. A p-value below 0.001 +means that the distributions differ. + +The performance plots time each operator for an increasing number of parents +with 30 variables, and for an increasing number of variables with 100 parents. +Each package receives the parents in its own data structure, which is created +before the timing starts. When a package expects the caller to choose the +other parents of DE or PCX, choosing them is part of the time. + +.. toctree:: + :maxdepth: 1 + + uniform + single_point + blend_alpha + single_arithmetic + local + bounded_exponential + de_rand_1_bin + de_current_1 + parent_centric + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/local.rst b/python/doc/source/benchmarks/crossover/local.rst new file mode 100644 index 000000000..9a42d18ff --- /dev/null +++ b/python/doc/source/benchmarks/crossover/local.rst @@ -0,0 +1,36 @@ +.. _bench-crossover-local: + +Local crossover +=============== + +:func:`~moocore.local_crossover` draws a random weight for each variable and +sets the two offspring to the two weighted averages of the parents. Only +`DESDEO`_ (``LocalCrossover``) provides this operator. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-local.csv + :header-rows: 1 + :stub-columns: 1 + +Both implementations produce the same distribution of offspring. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-local-nparents-time.png + :width: 80% + +`DESDEO`_ takes 39 to 70 times as long as `moocore`_. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-local-nvars-time.png + :width: 80% + +`DESDEO`_ takes 28 to 71 times as long as `moocore`_. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/parent_centric.rst b/python/doc/source/benchmarks/crossover/parent_centric.rst new file mode 100644 index 000000000..49c8d6f5b --- /dev/null +++ b/python/doc/source/benchmarks/crossover/parent_centric.rst @@ -0,0 +1,40 @@ +.. _bench-crossover-parent-centric: + +Parent-centric crossover (PCX) +============================== + +:func:`~moocore.parent_centric_crossover` creates an offspring around each +target, from two more parents chosen at random. `pymoo`_ (``PCX``) and +`DESDEO`_ (``ParentCentricCrossover``) provide the same operator. `pymoo`_ +replaces values outside the bounds by random values, while `moocore`_ and +`DESDEO`_ clip them. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-parent_centric.csv + :header-rows: 1 + :stub-columns: 1 + +All implementations produce the same distribution of offspring. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-parent_centric-nparents-time.png + :width: 80% + +`moocore`_ is the fastest for every number of parents. With 10 000 parents, +it takes 4.0 ms, while `DESDEO`_ and `pymoo`_ take 3.7 and 12 times as long. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-parent_centric-nvars-time.png + :width: 80% + +`moocore`_ is the fastest for every number of variables. With 1000 variables, +`pymoo`_ and `DESDEO`_ take 2.1 and 3.8 times as long. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/single_arithmetic.rst b/python/doc/source/benchmarks/crossover/single_arithmetic.rst new file mode 100644 index 000000000..0fdcd44b8 --- /dev/null +++ b/python/doc/source/benchmarks/crossover/single_arithmetic.rst @@ -0,0 +1,36 @@ +.. _bench-crossover-single-arithmetic: + +Single arithmetic crossover +=========================== + +:func:`~moocore.single_arithmetic_crossover` sets one random variable of both +offspring to the average of the parents and copies the other variables. Only +`DESDEO`_ (``SingleArithmeticCrossover``) provides this operator. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-single_arithmetic.csv + :header-rows: 1 + :stub-columns: 1 + +Both implementations produce the same distribution of offspring. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-single_arithmetic-nparents-time.png + :width: 80% + +`DESDEO`_ takes 16 to 69 times as long as `moocore`_. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-single_arithmetic-nvars-time.png + :width: 80% + +`DESDEO`_ takes 31 to 164 times as long as `moocore`_. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/single_point.rst b/python/doc/source/benchmarks/crossover/single_point.rst new file mode 100644 index 000000000..bb1466006 --- /dev/null +++ b/python/doc/source/benchmarks/crossover/single_point.rst @@ -0,0 +1,42 @@ +.. _bench-crossover-single-point: + +Single-point crossover +====================== + +:func:`~moocore.single_point_crossover` draws a cut point and exchanges the +tails of the two parents. `pymoo`_ (``SinglePointCrossover``), `DEAP`_ +(``cxOnePoint``) and `DESDEO`_ (``SinglePointBinaryCrossover``) provide the +same operator. The operator of `DESDEO`_ only accepts binary variables, so all +packages cross over binary parents here: a parent of zeros and a parent of ones +in the KS tests, and random binary parents in the timings. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-single_point.csv + :header-rows: 1 + :stub-columns: 1 + +All implementations produce the same distribution of offspring. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-single_point-nparents-time.png + :width: 80% + +With 10 and 20 parents, `DEAP`_ is faster than `moocore`_. From 50 parents +on, `moocore`_ is the fastest. With 10 000 parents, it takes 0.14 ms, while +`DESDEO`_, `DEAP`_ and `pymoo`_ take 9.8, 18 and 270 times as long. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-single_point-nvars-time.png + :width: 80% + +`moocore`_ is the fastest for every number of variables. `DEAP`_ takes 2.5 to +5 times as long, and `pymoo`_ and `DESDEO`_ take 24 to 144 times as long. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/crossover/uniform.rst b/python/doc/source/benchmarks/crossover/uniform.rst new file mode 100644 index 000000000..145137bd6 --- /dev/null +++ b/python/doc/source/benchmarks/crossover/uniform.rst @@ -0,0 +1,41 @@ +.. _bench-crossover-uniform: + +Uniform crossover +================= + +:func:`~moocore.uniform_crossover` swaps each variable between the two +offspring with probability ``swap_prob``, which is 0.5 by default. `pymoo`_ +(``UniformCrossover``), `DEAP`_ (``cxUniform``) and `DESDEO`_ +(``UniformCrossover``) provide the same operator. + +KS test results +--------------- + +.. csv-table:: + :file: /_static/bench/crossover_ks-uniform.csv + :header-rows: 1 + :stub-columns: 1 + +All implementations produce the same distribution of offspring. + +Performance over number of parents +---------------------------------- + +.. image:: /_static/bench/crossover_bench-uniform-nparents-time.png + :width: 80% + +With 10 parents, the Python loop of `DEAP`_ over five pairs takes less time +than the fixed cost of a call to `moocore`_. From 20 parents on, `moocore`_ is +the fastest. With 10 000 parents, it takes 1.4 ms, while `DESDEO`_, `DEAP`_ +and `pymoo`_ take 1.7, 4.4 and 20 times as long. + +Performance over number of variables +------------------------------------ + +.. image:: /_static/bench/crossover_bench-uniform-nvars-time.png + :width: 80% + +`moocore`_ is the fastest for every number of variables. With 1000 variables, +`pymoo`_, `DEAP`_ and `DESDEO`_ take 2.0, 3.9 and 7.4 times as long. + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/dominance.rst b/python/doc/source/benchmarks/dominance.rst new file mode 100644 index 000000000..2041e4329 --- /dev/null +++ b/python/doc/source/benchmarks/dominance.rst @@ -0,0 +1,117 @@ +Nondominated points and sorting +=============================== + +.. _bench-ndom: + +Identifying nondominated points +------------------------------- + +The following plots compare the speed of finding nondominated solutions, equivalent to :func:`moocore.is_nondominated`. We test both ``keep_weakly=False`` (the default) and ``keep_weakly=True``. The former is only supported by `BoTorch`_ and `paretoset`_. Thus, for the purpose of benchmarking the other packages with ``keep_weakly=False``, we remove duplicated points from the data and exclude this time from the results. This makes all packages, except `moocore`_, `BoTorch`_ and `paretoset`_, seem faster than they actually are, because the user is forced to remove duplicated points when using them. + +For 2D and 3D, the plots show that `moocore`_ is 10 times faster than `DESDEO`_ and `paretoset`_ and 100-1000 times faster than the other packages. `fast-pareto`_ claims to implement a :math:`O(n\log n)` algorithm for 3D, but the benchmarks below indicate a :math:`O(n^2)` complexity, similar to other packages and significantly slower than `moocore`_. + +|ndom_bench-test2D-200k| |ndom_bench-ran3d-40k| + +.. |ndom_bench-test2D-200k| image:: /_static/bench/ndom_bench-test2D-200k-time.png + :width: 49% + +.. |ndom_bench-ran3d-40k| image:: /_static/bench/ndom_bench-ran3d-40k-time.png + :width: 49% + +|wndom_bench-test2D-200k| |wndom_bench-ran3d-40k| + +.. |wndom_bench-test2D-200k| image:: /_static/bench/wndom_bench-test2D-200k-time.png + :width: 49% + +.. |wndom_bench-ran3d-40k| image:: /_static/bench/wndom_bench-ran3d-40k-time.png + :width: 49% + + +For dimensions larger than 3, :func:`moocore.is_nondominated` uses the +best-known :math:`O(n\log^{m-2} n)` algorithm. `fast-pareto`_ claims to be +using a :math:`O(n\log^{m-2} n)` algorithm, but the benchmarks suggest that it +uses a :math:`O(m n^2)` algorithm, although slower than most other packages and +almost 100 times slower than `moocore`_. The following plots with +``keep_weakly=True`` show that `moocore`_ is at least 10 times faster than the +other packages, even for 10D. + +|wndom_bench-sphere-4d| |wndom_bench-convex-4d| + +|wndom_bench-sphere-5d| |wndom_bench-rmnk-10d| + +.. |wndom_bench-sphere-4d| image:: /_static/bench/wndom_bench-sphere-4d-time.png + :width: 49% + +.. |wndom_bench-convex-4d| image:: /_static/bench/wndom_bench-convex-4d-time.png + :width: 49% + +.. |wndom_bench-sphere-5d| image:: /_static/bench/wndom_bench-sphere-5d-time.png + :width: 49% + +.. |wndom_bench-rmnk-10d| image:: /_static/bench/wndom_bench-rmnk-10d-time.png + :width: 49% + + +Similar results are obtained with ``keep_weakly=False``, however, this setting +is only supported by `moocore`_, `BoTorch`_ and `paretoset`_, thus the other +packages are actually slower to use than reported here, because the user is +forced to remove duplicated points when using them. + +|ndom_bench-sphere-4d| |ndom_bench-convex-4d| + +|ndom_bench-sphere-5d| |ndom_bench-rmnk-10d| + +.. |ndom_bench-sphere-4d| image:: /_static/bench/ndom_bench-sphere-4d-time.png + :width: 49% + +.. |ndom_bench-convex-4d| image:: /_static/bench/ndom_bench-convex-4d-time.png + :width: 49% + +.. |ndom_bench-sphere-5d| image:: /_static/bench/ndom_bench-sphere-5d-time.png + :width: 49% + +.. |ndom_bench-rmnk-10d| image:: /_static/bench/ndom_bench-rmnk-10d-time.png + :width: 49% + + +.. _bench-ndsort: + +Nondominated Sorting (Pareto ranking) +------------------------------------- + +Nondominated sorting (or Pareto ranking) assigns a rank to each point according +to Pareto-optimality. It can be used to split a set of points into layers of +nondominated fronts. For more details, see :func:`moocore.pareto_rank`. This +function in `moocore`_ uses a :math:`O(n \log n)` algorithm for 2 dimensions +and a :math:`O(n^2\log n)` algorithm for 3 dimensions. As the plots show, this +makes `moocore`_ significantly faster than `paretoset`_, `pymoo`_ and +`DESDEO`_. With higher dimensions, `moocore`_ uses a :math:`O(n^2\log^{m-2}n)` +algorithm instead of the :math:`O(n^3)` algorithm in `paretoset`_ and `DESDEO`_ +(latter is compiled with `Numba`_). In these benchmarks points are sampled +uniformly at random within the hypercube. + +|ndsort_bench-ran-2d| |ndsort_bench-ran-3d| + +|ndsort_bench-ran-4d| |ndsort_bench-ran-5d| + +|ndsort_bench-ran-9d| |ndsort_bench-ran-10d| + +.. |ndsort_bench-ran-2d| image:: /_static/bench/ndsort_bench-ran-2d-time.png + :width: 49% + +.. |ndsort_bench-ran-3d| image:: /_static/bench/ndsort_bench-ran-3d-time.png + :width: 49% + +.. |ndsort_bench-ran-4d| image:: /_static/bench/ndsort_bench-ran-4d-time.png + :width: 49% + +.. |ndsort_bench-ran-5d| image:: /_static/bench/ndsort_bench-ran-5d-time.png + :width: 49% + +.. |ndsort_bench-ran-9d| image:: /_static/bench/ndsort_bench-ran-9d-time.png + :width: 49% + +.. |ndsort_bench-ran-10d| image:: /_static/bench/ndsort_bench-ran-10d-time.png + :width: 49% + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/index.rst b/python/doc/source/benchmarks/index.rst new file mode 100644 index 000000000..35ac32f46 --- /dev/null +++ b/python/doc/source/benchmarks/index.rst @@ -0,0 +1,47 @@ +.. _benchmarks: + +Benchmarks +========== + +The following pages compare the performance of `moocore`_, `pymoo`_, `BoTorch`_, `DEAP`_, `DESDEO`_, `paretoset`_, `Nevergrad`_, `jMetalPy`_, `fast-pareto`_, `Optuna`_, `moarchiving`_, and `patatune`_. + +Other Python packages are not included in the comparison because they are based +on these packages for the functionality benchmarked, so they are **at least as +slow** as them. For example `Xopt`_ and `BoFire`_ use `BoTorch`_, `pysamoo`_ is +an extension of `pymoo`_. We do not include functions from the Bayesian +optimization toolbox `trieste`_, because it is much slower than `BoTorch`_ and +too slow to run the benchmarks in a reasonable time +(https://github.com/secondmind-labs/trieste/issues/917). Similarly, we do not +include `paretobench`_ because it is always slower than `BoTorch`_ for +identifying nondominated points and more than 2000 times slower than `moocore`_ +for computing the hypervolume in 3D +(https://github.com/electronsandstuff/ParetoBench/issues/57), which is too slow +to include in the benchmarks. + +Not all packages provide the same functionality. For example, `pymoo`_ does not +provide the :ref:`epsilon indicator `. `BoTorch`_ and +`paretobench`_ only provide the hypervolume. `paretoset`_ and `fast-pareto`_ +only identify nondominated points. `DESDEO`_, `DEAP`_, and `jMetalPy`_ use +`moocore`_ for computing the hypervolume, but other functionality, such as +filtering dominated points, is still slower than `moocore`_. Recent versions +of `seqme`_ and `pymoo`_ (≥0.6.2) already use `moocore`_ for most functionality +benchmarked here. + +We would like to benchmark `pygmo`_, however, it is currently impossible to install using ``pip`` (See https://github.com/esa/pygmo2/issues/152). + +We would also like to benchmark `platypus`_, but its API is too different from +other packages. For example, it only supports filtering dominated points by +creating a new list of points, but it cannot identify their position in the +input. It cannot calculate hypervolume without normalization. It doesn't +support :obj:`~numpy.typing.ArrayLike` inputs. + +The source code for the benchmarks below can be found at https://github.com/multi-objective/moocore/tree/main/python/benchmarks . + +.. toctree:: + :maxdepth: 2 + + dominance + indicators + crossover/index + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/benchmarks/indicators.rst b/python/doc/source/benchmarks/indicators.rst new file mode 100644 index 000000000..587d0c929 --- /dev/null +++ b/python/doc/source/benchmarks/indicators.rst @@ -0,0 +1,114 @@ +Quality indicators +================== + +.. _bench-hv: + +Exact computation of hypervolume +-------------------------------- + +The following plots compare the speed of computing the :ref:`hypervolume indicator ` in 3D, 4D, 5D and 6D. As the plots show, `moocore`_ is 100 times faster than the other packages and 1000 times faster than `BoTorch`_ and, by extension, `Xopt`_ and `BoFire`_. + +For 3D, `Optuna`_ uses a :math:`O(n^2)` algorithm, while `moocore`_ uses the best-known :math:`O(n\log n)` one, so the gap between the two will get significantly larger with a larger number of points than those tested here. `moarchiving`_ uses the same algorithms as `moocore`_ for 3D and 4D, but a pure Python implementation, so it is still 10-100 times slower than `moocore`_'s C implementation. `BoTorch`_ and `fast-pareto`_ are not included for more than 3 objectives because **they are tens of thousands of times slower** than `moocore`_. + +Recent versions of `DESDEO`_, `DEAP`_ (≥1.4.4), `pymoo`_ (≥0.6.1.6) and `jMetalPy`_ (≥1.9.0) already use `moocore`_ for hypervolume computation. +The implementation in `patatune`_ appears to be wrong (https://github.com/cms-patatrack/patatune/issues/52). + +|hv_bench-DTLZLinearShape-3d| |hv_bench-DTLZLinearShape-4d| + +|hv_bench-DTLZLinearShape-5d| |hv_bench-DTLZLinearShape-6d| + +.. |hv_bench-DTLZLinearShape-3d| image:: /_static/bench/hv_bench-DTLZLinearShape.3d-time.png + :width: 49% + +.. |hv_bench-DTLZLinearShape-4d| image:: /_static/bench/hv_bench-DTLZLinearShape.4d-time.png + :width: 49% + +.. |hv_bench-DTLZLinearShape-5d| image:: /_static/bench/hv_bench-DTLZLinearShape.5d-time.png + :width: 49% + +.. |hv_bench-DTLZLinearShape-6d| image:: /_static/bench/hv_bench-DTLZLinearShape.6d-time.png + :width: 49% + +.. _bench-hvapprox: + +Approximation of the hypervolume +-------------------------------- + +The following plots compare the accuracy and speed of approximating the +hypervolume with the various methods provided by :func:`moocore.hv_approx`. The +plots show that there is no clear winner, in terms of approximation error, +between methods ``Rphi-FWE+`` (default) and ``DZ2019-HW``, but both produce +consistently lower approximation errors than method ``DZ2019-MC``. However, ``Rphi-FWE+`` is as fast as ``DZ2019-MC`` and both are consistently faster than ``DZ2019-HW``, in particular with higher number of objectives. + +The exact computation of the hypervolume for less than 5D takes milliseconds +for thousands of points, whereas approximating the hypervolume is significantly +slower and, thus, not worth doing. + +Approximating the hypervolume becomes more useful for dimensions higher than 5, +where the exact computation becomes noticeably slower with hundreds of points. +For such problems, method ``DZ2019-HW`` becomes significantly slower than +``Rphi-FWE+``. + +|hvapprox_bench-DTLZLinearShape-6d-values| |hvapprox_bench-DTLZLinearShape-6d-time| + +|hvapprox_bench-DTLZSphereShape-6d-values| |hvapprox_bench-DTLZSphereShape-6d-time| + +|hvapprox_bench-ran-6d-values| |hvapprox_bench-ran-6d-time| + +|hvapprox_bench-DTLZLinearShape-9d-values| |hvapprox_bench-DTLZLinearShape-9d-time| + +|hvapprox_bench-ran-9d-values| |hvapprox_bench-ran-9d-time| + +|hvapprox_bench-DTLZSphereShape-10d-values| |hvapprox_bench-DTLZSphereShape-10d-time| + +.. |hvapprox_bench-DTLZSphereShape-6d-values| image:: /_static/bench/hvapprox_bench-DTLZSphereShape.6d-values.png + :width: 49% + +.. |hvapprox_bench-DTLZSphereShape-6d-time| image:: /_static/bench/hvapprox_bench-DTLZSphereShape.6d-time.png + :width: 49% + +.. |hvapprox_bench-DTLZLinearShape-6d-values| image:: /_static/bench/hvapprox_bench-DTLZLinearShape.6d-values.png + :width: 49% + +.. |hvapprox_bench-DTLZLinearShape-6d-time| image:: /_static/bench/hvapprox_bench-DTLZLinearShape.6d-time.png + :width: 49% + +.. |hvapprox_bench-ran-6d-values| image:: /_static/bench/hvapprox_bench-ran.6d-values.png + :width: 49% + +.. |hvapprox_bench-ran-6d-time| image:: /_static/bench/hvapprox_bench-ran.6d-time.png + :width: 49% + +.. |hvapprox_bench-DTLZLinearShape-9d-values| image:: /_static/bench/hvapprox_bench-DTLZLinearShape.9d-values.png + :width: 49% + +.. |hvapprox_bench-DTLZLinearShape-9d-time| image:: /_static/bench/hvapprox_bench-DTLZLinearShape.9d-time.png + :width: 49% + +.. |hvapprox_bench-ran-9d-values| image:: /_static/bench/hvapprox_bench-ran.9d-values.png + :width: 49% + +.. |hvapprox_bench-ran-9d-time| image:: /_static/bench/hvapprox_bench-ran.9d-time.png + :width: 49% + +.. |hvapprox_bench-DTLZSphereShape-10d-values| image:: /_static/bench/hvapprox_bench-DTLZSphereShape.10d-values.png + :width: 49% + +.. |hvapprox_bench-DTLZSphereShape-10d-time| image:: /_static/bench/hvapprox_bench-DTLZSphereShape.10d-time.png + :width: 49% + + +Epsilon and IGD+ indicators +--------------------------- + +The following plots benchmark the speed of computing the :ref:`epsilon indicator ` metric and :ref:`IGD+ indicator `. Although the algorithms for computing these metrics are relatively simple and easy to vectorize in Python, the `moocore`_ implementation is still 10 to 100 times faster. + +|fig_eps_bench| |fig_igd_plus_bench| + +.. |fig_eps_bench| image:: /_static/bench/eps_bench-rmnk_10D_random_search-time.png + :width: 48% + +.. |fig_igd_plus_bench| image:: /_static/bench/igd_plus_bench-ran.1000pts.3d-time.png + :width: 48% + +.. include:: /benchmarks-links.inc.rst diff --git a/python/doc/source/index.rst b/python/doc/source/index.rst index 9e5c57014..abec66d90 100644 --- a/python/doc/source/index.rst +++ b/python/doc/source/index.rst @@ -10,6 +10,7 @@ moocore: Core Algorithms for Multi-Objective Optimization API reference Examples + Benchmarks whatsnew/index contribute/index @@ -30,6 +31,7 @@ The goal of the **moocore** project (https://github.com/multi-objective/moocore/ * :ref:`Identify and filter dominated vectors `. * :ref:`Quality metrics ` such as (weighted) hypervolume, epsilon, IGD+, (exact) R2, etc. * :ref:`Computation of the Empirical Attainment Function `. The empirical attainment function (EAF) describes the probabilistic distribution of the outcomes obtained by a stochastic algorithm in the objective space. +* :ref:`Crossover operators ` for evolutionary algorithms, such as uniform, blend-alpha, bounded exponential, differential evolution and parent-centric crossover. Most critical functionality is implemented in C, with the R and Python packages providing convenient interfaces to the C code. @@ -62,7 +64,8 @@ performance measures, performance assessment Detailed examples and tutorials. -.. include:: benchmarks.inc.rst +The :ref:`benchmarks ` compare the speed of moocore with other +Python packages. .. This is not really the index page, that is found in diff --git a/python/doc/source/my_unused_images.py b/python/doc/source/my_unused_images.py index b440e3cb1..ec5cce5c1 100644 --- a/python/doc/source/my_unused_images.py +++ b/python/doc/source/my_unused_images.py @@ -23,7 +23,10 @@ def check_unused_images(app, exception): if "://" in uri or uri.startswith("data:"): continue + # Sphinx stores the URI relative to the source directory, which is + # not the folder of the document for documents in subfolders. used.add((document.parent / uri).resolve()) + used.add((Path(app.srcdir) / uri).resolve()) for image in Path(app.srcdir).rglob("*.png"): if "auto_examples" in image.parts: diff --git a/python/doc/source/reference/functions.variation.rst b/python/doc/source/reference/functions.variation.rst new file mode 100644 index 000000000..6dbd4849d --- /dev/null +++ b/python/doc/source/reference/functions.variation.rst @@ -0,0 +1,131 @@ +.. _variation_operators: + +=================== +Variation operators +=================== + +.. currentmodule:: moocore + +Variation operators are used in evolutionary algorithms to +produce new decision vectors, called offspring, from +existing ones, called parents. +This usually happens in two steps: first, a crossover operator combines the +parents to produce one or more offspring, and then a mutation operator randomly +changes some of the offspring's decision variables. We currently provide a +number of crossover operators. +Every operator takes a 2D array with one +decision vector per row and accepts a ``seed`` for reproducible results. The +:ref:`crossover example ` shows +the behaviour of the crossover operators. + +Crossover +========= +.. autosummary:: + :toctree: generated/ + + uniform_crossover + single_point_crossover + blend_alpha_crossover + single_arithmetic_crossover + local_crossover + bounded_exponential_crossover + de_rand_1_bin_crossover + de_current_1_crossover + parent_centric_crossover + +The operators differ in how many parents one crossover operation uses and how +many offspring it produces: + +.. list-table:: + :header-rows: 1 + + * - Function + - Offspring per operation + - Parents per operation + * - :func:`uniform_crossover` + - 2 + - 2 + * - :func:`single_point_crossover` + - 2 + - 2 + * - :func:`blend_alpha_crossover` + - 1 + - 2 + * - :func:`single_arithmetic_crossover` + - 2 + - 2 + * - :func:`local_crossover` + - 2 + - 2 + * - :func:`bounded_exponential_crossover` + - 2 + - 2 + * - :func:`de_rand_1_bin_crossover` + - 1 + - 1, plus 3 other parents chosen at random + * - :func:`de_current_1_crossover` + - 1 + - 1, plus 2 other parents chosen at random + * - :func:`parent_centric_crossover` + - 1 + - 1, plus 2 other parents chosen at random + +The functions that use 2 parents per operation take a matrix ``parents`` and +cross over rows ``2i`` and ``2i + 1``, so the order of the rows decides which +parents are crossed over, for example after a selection step. If the number +of rows is odd, the last row is crossed over with the first one. +When an +operation produces 2 offspring, the output has as many rows as ``parents``, and +one offspring is dropped when the number of rows is odd. +:func:`blend_alpha_crossover` returns one row per pair. + +The functions that use 1 parent per operation take a matrix ``population`` and +the row indices ``targets`` of the parents, and return one offspring per target. +They choose the other parents at random from ``population``. + +Differential evolution (DE) comes in two variants that are common in +multi-objective optimization. :func:`de_rand_1_bin_crossover` is the +DE/rand/1/bin scheme of :footcite:t:`StoPri1997:de`, which most Pareto-based +multi-objective DE algorithms use :footcite:p:`BezLopStu2015emode`. +:func:`de_current_1_crossover` is the operator of MOEA/D-DE +:footcite:p:`LiZha2009:moead-de`. It starts from the target instead of a +random row and does not force any variable to change. Its 2 other parents may +include the target itself. `Tanabe and Ishibuchi (2019) +`__ call this mutation the current/1 +strategy, a term introduced by `Gong et al. (2017) +`__. + +Certain crossover operators tend to produce offspring that are close to one +or the other parent. +Operators that produce 2 offspring from 2 parents or 1 offspring from 1 parent +preserve the order, so that the index of an offspring corresponds to the index +of the parent it is closest to, assuming a parent-centric crossover operation. +For the operators that take ``population``, the parent of offspring ``i`` is +row ``targets[i]``. + + +Random seeds +============ + +Every operator takes a ``seed``. With an integer seed, the operator draws the +same random numbers on every call. This makes a single call reproducible, but +in an evolutionary algorithm every generation would reuse the same random +numbers. Instead, create one :class:`numpy.random.Generator` at the start of a +run and pass it as ``seed`` to every call. Each call then draws a new seed +from the generator, and the seed of the generator still reproduces the whole +run:: + + rng = np.random.default_rng(42) + for generation in range(n_generations): + parents = ... # Select and order the parents. + offspring = moocore.uniform_crossover(parents, seed=rng) + ... + +With ``seed=None``, every call uses a new random seed, so the results cannot be +reproduced. + + +Bibliography +============ + +.. footbibliography:: diff --git a/python/doc/source/reference/index.rst b/python/doc/source/reference/index.rst index c28449611..f2123fd85 100644 --- a/python/doc/source/reference/index.rst +++ b/python/doc/source/reference/index.rst @@ -29,3 +29,4 @@ indentation. functions.dominance functions.metrics functions.eaf + functions.variation diff --git a/python/doc/source/whatsnew/index.rst b/python/doc/source/whatsnew/index.rst index ea58ce81e..de39f4e4a 100644 --- a/python/doc/source/whatsnew/index.rst +++ b/python/doc/source/whatsnew/index.rst @@ -12,6 +12,7 @@ Version 0.4.0 - :func:`~moocore.is_nondominated` is up to 10x faster in some inputs thanks to a customized radixsort implementation. - New shapes ``"concave-cliff"`` and ``"convex-cliff"`` added to :func:`~moocore.generate_ndset`. - New function :func:`~moocore.generate_sequence` to generate sequences of dominated and nondominated points. +- New :ref:`crossover operators `: :func:`~moocore.uniform_crossover`, :func:`~moocore.single_point_crossover`, :func:`~moocore.blend_alpha_crossover`, :func:`~moocore.single_arithmetic_crossover`, :func:`~moocore.local_crossover`, :func:`~moocore.bounded_exponential_crossover`, :func:`~moocore.de_rand_1_bin_crossover`, :func:`~moocore.de_current_1_crossover` and :func:`~moocore.parent_centric_crossover`. Version 0.3.2 (11/07/2026) diff --git a/python/examples/plot_crossover.py b/python/examples/plot_crossover.py new file mode 100644 index 000000000..71ecac85c --- /dev/null +++ b/python/examples/plot_crossover.py @@ -0,0 +1,287 @@ +r""" +Crossover Operators +=================== + +This example shows the behaviour of the crossover operators, first with two +fixed parents and then with a whole population. Each figure crosses over the +same parents many times and plots all the offspring. + +First we define a few functions useful for plotting. +""" + +import matplotlib.pyplot as plt +import numpy as np +import moocore + +LOWER, UPPER = 0.0, 10.0 +BOUNDS = dict(lower=LOWER, upper=UPPER) + + +def cross_over(fun, a, b, n_pairs=50_000, **kwargs): + """Cross over parents ``a`` and ``b`` ``n_pairs`` times.""" + parents = np.tile(np.array([a, b], dtype=float), (n_pairs, 1)) + return fun(parents, seed=42, **kwargs) + + +def plot_density(ax, values, **kwargs): + """Plot the distribution of ``values`` within the bounds.""" + ax.hist( + values, + bins=200, + range=(LOWER, UPPER), + density=True, + histtype="step", + **kwargs, + ) + + +def mark_parents(ax, positions): + """Show the parents as dotted vertical lines.""" + for x in positions: + ax.axvline(x, color="gray", linestyle=":", linewidth=1) + + +# %% +# +# Crossover behaviour with two parents +# ------------------------------------ +# +# The parents are :math:`(3, 3)` and :math:`(7, 7)` within :math:`[0, 10]^2`. +# Each panel shows the distribution of the first variable of the offspring for +# the default parameters of each operator. Uniform and single arithmetic +# crossover only produce a few distinct values, which appear as spikes. +# The parent decision variables are shown as dotted vertical lines. + +a, b = [3.0, 3.0], [7.0, 7.0] +operators = { + "Bounded exponential": (moocore.bounded_exponential_crossover, BOUNDS), + "Blend alpha": (moocore.blend_alpha_crossover, BOUNDS), + "Local": (moocore.local_crossover, {}), + "Uniform": (moocore.uniform_crossover, {}), + "Single arithmetic": (moocore.single_arithmetic_crossover, {}), +} +fig, axes = plt.subplots(len(operators), 1, figsize=(6, 8), sharex=False) +for i, (ax, (name, (fun, kwargs))) in enumerate(zip(axes, operators.items())): + plot_density(ax, cross_over(fun, a, b, **kwargs)[:, 0], color=f"C{i}") + mark_parents(ax, [a[0], b[0]]) + ax.set_title(name, loc="left", fontsize="medium") + ax.set_yticks([]) + ax.set_xlabel("$x_1$") +fig.tight_layout() +plt.show() + +# %% +# +# Bounded exponential crossover +# -------------------------------------------------- +# +# The bounded exponential crossover operator samples from an exponential +# distribution that is truncated at the bounds. The distribution is centred on +# one parent, and the distance to the other parent scales its spread. The +# offspring moves towards either bound with equal probability. The operator is +# symmetric, so the second offspring is centred on the other parent in the same +# way. The parameter ``scale`` (:math:`\lambda`) controls the decay of the +# distribution. + +fig, ax = plt.subplots(figsize=(6, 3.5)) +for i, scale in enumerate([0.01, 0.05, 0.1, 0.5]): + offspring = cross_over( + moocore.bounded_exponential_crossover, a, b, scale=scale, **BOUNDS + ) + plot_density( + ax, + offspring[:, 0], + color=plt.cm.viridis(i / 4), + label=rf"$\lambda = {scale}$", + ) +mark_parents(ax, [a[0], b[0]]) +ax.set_xlabel("$x_1$") +ax.set_ylabel("Density") +ax.legend() +fig.tight_layout() +plt.show() + +# %% +# +# Blend crossover: the alpha parameter +# ------------------------------------ +# +# BLX-:math:`\alpha` samples uniformly from the interval between the parents, +# widened by ``alpha`` times their distance on each side. + +fig, ax = plt.subplots(figsize=(6, 3.5)) +for i, alpha in enumerate([0.0, 0.25, 0.5, 1.0]): + offspring = cross_over( + moocore.blend_alpha_crossover, a, b, alpha=alpha, **BOUNDS + ) + plot_density( + ax, + offspring[:, 0], + color=plt.cm.viridis(i / 4), + label=rf"$\alpha = {alpha}$", + ) +mark_parents(ax, [a[0], b[0]]) +ax.set_xlabel("$x_1$") +ax.set_ylabel("Density") +ax.legend() +fig.tight_layout() +plt.show() + +# %% +# +# Offspring geometry in two dimensions +# ------------------------------------ +# +# With the parents at opposite corners, :math:`(2, 8)` and :math:`(8, 2)`, the +# operators differ in which regions of the decision space their offspring +# can reach. Note that certain operators such as uniform crossover can never +# explore all possible decision vectors, requiring the use of mutation for +# proper exploration of the decision space. Moreover, the local crossover +# operator never produces offspring outside the bounding box defined by the +# parents, so it is not suitable for exploring the decision space beyond +# the current population. + +a, b = [2.0, 8.0], [8.0, 2.0] +operators = { + "Bounded exponential": (moocore.bounded_exponential_crossover, BOUNDS), + "Blend alpha": (moocore.blend_alpha_crossover, BOUNDS), + "Local": (moocore.local_crossover, {}), + "Uniform": (moocore.uniform_crossover, {}), + "Single arithmetic": (moocore.single_arithmetic_crossover, {}), + "Single point": (moocore.single_point_crossover, {}), +} +fig, axes = plt.subplots(2, 3, figsize=(9, 6), sharex=True, sharey=True) +for i, (ax, (name, (fun, kwargs))) in enumerate( + zip(axes.flat, operators.items()) +): + offspring = cross_over(fun, a, b, n_pairs=1000, **kwargs) + ax.scatter(offspring[:, 0], offspring[:, 1], s=4, alpha=0.3, color=f"C{i}") + ax.scatter(*np.array([a, b]).T, marker="D", s=40, color="black") + ax.set_title(name, fontsize="medium") + ax.set_xlim(LOWER, UPPER) + ax.set_ylim(LOWER, UPPER) + ax.set_aspect("equal") +for ax in axes[1]: + ax.set_xlabel("$x_1$") +for ax in axes[:, 0]: + ax.set_ylabel("$x_2$") +fig.tight_layout() +plt.show() + +# %% +# +# Crossover behaviour with a population +# ------------------------------------- +# +# In an evolutionary algorithm, crossover is applied to a whole population. +# Here, each operator is applied 40 times to a population of 60 decision +# vectors (circles), and all the offspring are plotted. +# Each operator gets the same population, but for each call, the population is +# randomly permuted and/or a different seed is used. +# +# The first population is a round cloud. Uniform, single-point and single +# arithmetic crossover place their offspring on lines parallel to the axes, +# because each offspring copies at least one variable from a parent. Uniform +# and single-point crossover look alike here only because, with two variables, +# the cut point always falls between them. With more variables, their +# distributions differ. Both DE variants spread their offspring beyond the +# population, whereas parent-centric and bounded exponential crossover keep +# them close to their parents. + +all_operators = { + "Uniform": (moocore.uniform_crossover, {}), + "Single point": (moocore.single_point_crossover, {}), + "Blend alpha": (moocore.blend_alpha_crossover, BOUNDS), + "Single arithmetic": (moocore.single_arithmetic_crossover, {}), + "Local": (moocore.local_crossover, {}), + "Bounded exponential": (moocore.bounded_exponential_crossover, BOUNDS), + "DE/rand/1/bin": (moocore.de_rand_1_bin_crossover, BOUNDS), + "DE/current/1": (moocore.de_current_1_crossover, BOUNDS), + "Parent-centric": (moocore.parent_centric_crossover, BOUNDS), +} +takes_population = { + moocore.de_rand_1_bin_crossover, + moocore.de_current_1_crossover, + moocore.parent_centric_crossover, +} + + +def plot_population(population, n_calls=40): + """Plot the offspring of every operator applied to ``population``.""" + rng = np.random.default_rng(42) + fig, axes = plt.subplots(3, 3, figsize=(9, 9), sharex=True, sharey=True) + for i, (ax, (name, (fun, kwargs))) in enumerate( + zip(axes.flat, all_operators.items()) + ): + offspring = [] + for _ in range(n_calls): + if fun in takes_population: + x = population + else: + x = rng.permutation(population) + offspring.append(fun(x, seed=rng, **kwargs)) + offspring = np.vstack(offspring) + ax.scatter( + offspring[:, 0], offspring[:, 1], s=3, alpha=0.3, color=f"C{i}" + ) + ax.scatter( + population[:, 0], + population[:, 1], + s=15, + facecolors="none", + edgecolors="black", + linewidths=0.8, + ) + ax.set_title(name, fontsize="medium") + ax.set_xlim(LOWER, UPPER) + ax.set_ylim(LOWER, UPPER) + ax.set_aspect("equal") + for ax in axes[2]: + ax.set_xlabel("$x_1$") + for ax in axes[:, 0]: + ax.set_ylabel("$x_2$") + fig.tight_layout() + plt.show() + + +rng = np.random.default_rng(42) +round_population = np.clip(rng.normal([4.0, 6.0], 1.2, size=(60, 2)), 0.2, 9.8) +plot_population(round_population) + +# %% +# +# The second population is stretched along a diagonal, as may happen when a +# population approaches a Pareto set. Uniform and single-point crossover +# combine the variables of different parents independently, so their offspring +# fill the box around the population and lose its diagonal shape. Blend alpha, +# single arithmetic and local crossover also treat each variable separately and +# place offspring away from the diagonal. Bounded exponential crossover keeps +# the offspring close to their parents. Both DE variants and parent-centric +# crossover follow the diagonal, because they move along differences between +# parents. DE extends the population beyond its ends, whereas parent-centric +# crossover stays close to its index parents. + +t = rng.uniform(0, 1, size=60) +diagonal_population = np.column_stack([2 + 6 * t, 8 - 6 * t]) +diagonal_population += rng.normal(0, 0.25, size=(60, 2)) +plot_population(diagonal_population) + +# %% +# +# The third population consists of three separate clusters of 20 members each, +# as when the Pareto set is disconnected. Crossing over members of different +# clusters produces offspring outside the clusters, unless the operator keeps +# the offspring close to their parents. Uniform and single-point crossover +# create new clusters at every combination of the clusters' variable +# values, and single arithmetic crossover adds clusters at the midpoints of one +# variable. Blend alpha and local crossover fill the space between the +# clusters. Both DE variants create new clusters by shifting the original ones +# by :math:`F` times the differences between clusters, and the offspring that +# are clipped to the bounds pile up on the edges of the box. Bounded +# exponential and parent-centric crossover keep most offspring near the cluster +# of their parent. + +centres = np.array([[2.0, 8.0], [8.0, 6.0], [5.0, 2.0]]) +clustered_population = np.repeat(centres, 20, axis=0) +clustered_population += rng.normal(0, 0.4, size=(60, 2)) +plot_population(clustered_population) diff --git a/python/pyproject.toml b/python/pyproject.toml index e3672ff31..1fff08067 100644 --- a/python/pyproject.toml +++ b/python/pyproject.toml @@ -40,6 +40,7 @@ dependencies = [ [project.optional-dependencies] benchmarks = [ "botorch", + "deap", "desdeo", "fast-pareto", "jmetalpy>=1.8", @@ -70,7 +71,7 @@ docs = [ "plotly", "pydata-sphinx-theme>=0.16", "seaborn", - "sphinx>=6", + "sphinx>=9", "sphinx-autodoc-typehints", "sphinx-copybutton", "sphinx-design", diff --git a/python/src/moocore/__init__.py b/python/src/moocore/__init__.py index d07f559f0..1e7aaa5ba 100644 --- a/python/src/moocore/__init__.py +++ b/python/src/moocore/__init__.py @@ -38,6 +38,18 @@ from ._generate import generate_ndset, generate_sequence +from ._crossover import ( + blend_alpha_crossover, + bounded_exponential_crossover, + de_current_1_crossover, + de_rand_1_bin_crossover, + local_crossover, + parent_centric_crossover, + single_arithmetic_crossover, + single_point_crossover, + uniform_crossover, +) + from importlib.metadata import version as _metadata_version __version__ = _metadata_version("moocore") @@ -52,6 +64,10 @@ "any_dominated", "apply_within_sets", "avg_hausdorff_dist", + "blend_alpha_crossover", + "bounded_exponential_crossover", + "de_current_1_crossover", + "de_rand_1_bin_crossover", "eaf", "eafdiff", "epsilon_additive", @@ -70,11 +86,16 @@ "is_nondominated", "is_nondominated_within_sets", "largest_eafdiff", + "local_crossover", "normalise", + "parent_centric_crossover", "pareto_rank", "r2_exact", "read_datasets", + "single_arithmetic_crossover", + "single_point_crossover", "total_whv_rect", + "uniform_crossover", "vorob_dev", "vorob_t", "whv_hype", diff --git a/python/src/moocore/_crossover.py b/python/src/moocore/_crossover.py new file mode 100644 index 000000000..627ce32dd --- /dev/null +++ b/python/src/moocore/_crossover.py @@ -0,0 +1,895 @@ +from __future__ import annotations + +from collections.abc import Callable +from typing import Any + +import numpy as np +from numpy.typing import ArrayLike + +from ._docsubstitute import DocSubstitute +from ._libmoocore import ffi, lib # ty: ignore[unresolved-import] +from ._utils import _get_seed_for_c, array_1d_of_length_n, np2d_to_double_array + + +def _as_matrix(x: ArrayLike, name: str, min_rows: int) -> np.ndarray: + x = np.asarray(x, dtype=float) + if x.ndim != 2 or x.shape[1] == 0: + raise ValueError(f"{name!r} must be a 2D array with at least 1 column") + if x.shape[0] < min_rows: + raise ValueError(f"{name!r} must have at least {min_rows} rows") + return x + + +def _check_probability(value: float, name: str) -> float: + if not 0 <= value <= 1: + raise ValueError(f"{name!r} must be within [0, 1]") + return float(value) + + +def _check_positive(value: float, name: str) -> float: + if not value > 0: + raise ValueError(f"{name!r} must be larger than 0") + return float(value) + + +def _bounds( + lower: ArrayLike | None, upper: ArrayLike | None, nvar: int +) -> tuple[np.ndarray, np.ndarray]: + if lower is None: + lower = np.full(nvar, -np.inf) + else: + lower = array_1d_of_length_n( + np.asarray(lower, dtype=float), nvar, name="lower" + ) + if upper is None: + upper = np.full(nvar, np.inf) + else: + upper = array_1d_of_length_n( + np.asarray(upper, dtype=float), nvar, name="upper" + ) + if np.any(lower > upper): + raise ValueError("'lower' must not be larger than 'upper'") + return np.ascontiguousarray(lower), np.ascontiguousarray(upper) + + +def _check_within_bounds( + parents: np.ndarray, lower: np.ndarray, upper: np.ndarray +) -> None: + # NaN fails both comparisons, so it is rejected too. + if not np.all((parents >= lower) & (parents <= upper)): + raise ValueError("'parents' must be within 'lower' and 'upper'") + + +def _targets(targets: ArrayLike | None, n: int) -> np.ndarray: + if targets is None: + return np.arange(n, dtype=np.uintp) + targets = np.atleast_1d(targets) + if targets.ndim != 1 or ( + targets.size > 0 and not np.issubdtype(targets.dtype, np.integer) + ): + raise ValueError("'targets' must be a 1D vector of row indices") + if np.any((targets < 0) | (targets >= n)): + raise ValueError(f"'targets' must be row indices within [0, {n - 1}]") + return np.ascontiguousarray(targets, dtype=np.uintp) + + +def _call( + fun: Callable[..., None], x: np.ndarray, nrows: int, *args: Any +) -> np.ndarray: + """Call the C function ``fun`` and return the offspring it writes.""" + out = np.empty((nrows, x.shape[1]), dtype=float) + x_p, n, nvar = np2d_to_double_array(x, ctype_shape=("size_t", "size_t")) + fun(ffi.from_buffer("double []", out), x_p, n, nvar, *args) + return out + + +def _keep_dtype(out: np.ndarray, dtype: np.dtype) -> np.ndarray: + if np.issubdtype(dtype, np.integer) or np.issubdtype(dtype, np.bool_): + return out.astype(dtype) + return out + + +def _call_with_targets( + fun: Callable[..., None], + population: np.ndarray, + targets: np.ndarray, + lower: np.ndarray, + upper: np.ndarray, + *args: Any, +) -> np.ndarray: + return _call( + fun, + population, + len(targets), + ffi.from_buffer("size_t []", targets), + len(targets), + ffi.from_buffer("double []", lower), + ffi.from_buffer("double []", upper), + *args, + ) + + +@DocSubstitute() +def uniform_crossover( + parents: ArrayLike, + /, + *, + swap_prob: float = 0.5, + pair_prob: float = 1.0, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Uniform crossover. + + Each crossover operation on a pair of parents produces two offspring. With + probability ``pair_prob``, each variable is swapped between the two + offspring with probability ``swap_prob``. Otherwise, the offspring are + copies of their parents. + + The offspring only take values from their parents, so this operator works + for binary, integer and real-valued variables alike. Integer and Boolean + inputs keep their dtype. + + Parameters + ---------- + parents : + ${parents} + swap_prob : + ${swap_prob} + pair_prob : + ${pair_prob} + seed : + ${random_seed} + + Returns + ------- + ${pair_offspring} + + See Also + -------- + single_point_crossover : Exchange the tails of two parents. + + Notes + ----- + Uniform crossover was proposed by :footcite:t:`Syswerda89`. See also + :footcite:t:`EibSmi2003`. + + References + ---------- + .. footbibliography:: + + Examples + -------- + >>> parents = np.array([[0, 0, 0, 0, 0, 0], [1, 1, 1, 1, 1, 1]]) + >>> moocore.uniform_crossover(parents, seed=5) + array([[0, 1, 0, 1, 0, 0], + [1, 0, 1, 0, 1, 1]]) + + """ + dtype = np.asarray(parents).dtype + parents = _as_matrix(parents, "parents", min_rows=2) + swap_prob = _check_probability(swap_prob, "swap_prob") + pair_prob = _check_probability(pair_prob, "pair_prob") + out = _call( + lib.uniform_crossover, + parents, + len(parents), + swap_prob, + pair_prob, + _get_seed_for_c(seed), + ) + return _keep_dtype(out, dtype) + + +@DocSubstitute() +def single_point_crossover( + parents: ArrayLike, + /, + *, + pair_prob: float = 1.0, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Single-point crossover. + + Each crossover operation on a pair of parents produces two offspring. With + probability ``pair_prob``, a cut point :math:`k` is drawn uniformly from + :math:`\{1, \dots, n_\text{var} - 1\}`. The first offspring takes the + first :math:`k` variables from the first parent and the rest from the + second parent, and the second offspring takes the complement. Otherwise, + the offspring are copies of their parents. + + The offspring only take values from their parents, so this operator works + for binary, integer and real-valued variables alike. Integer and Boolean + inputs keep their dtype. + + Parameters + ---------- + parents : + ${parents} It must have at least 2 columns. + pair_prob : + ${pair_prob} + seed : + ${random_seed} + + Returns + ------- + ${pair_offspring} + + See Also + -------- + uniform_crossover : Swap each variable independently. + + Notes + ----- + Single-point crossover is the classical crossover of genetic algorithms + :footcite:p:`Hol75,Goldberg89`. + + References + ---------- + .. footbibliography:: + + Examples + -------- + >>> parents = np.array([[0, 0, 0, 0, 0, 0], [1, 1, 1, 1, 1, 1]]) + >>> moocore.single_point_crossover(parents, seed=4) + array([[0, 0, 0, 1, 1, 1], + [1, 1, 1, 0, 0, 0]]) + + """ + dtype = np.asarray(parents).dtype + parents = _as_matrix(parents, "parents", min_rows=2) + if parents.shape[1] < 2: + raise ValueError("'parents' must have at least 2 columns") + pair_prob = _check_probability(pair_prob, "pair_prob") + out = _call( + lib.single_point_crossover, + parents, + len(parents), + pair_prob, + _get_seed_for_c(seed), + ) + return _keep_dtype(out, dtype) + + +@DocSubstitute() +def blend_alpha_crossover( + parents: ArrayLike, + /, + *, + lower: ArrayLike | None = None, + upper: ArrayLike | None = None, + alpha: float = 0.5, + pair_prob: float = 1.0, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Blend crossover (BLX-:math:`\alpha`). + + Each crossover operation on a pair of parents :math:`x, y` produces + **one** offspring :math:`z`. With probability ``pair_prob``, for each + variable :math:`j`, let :math:`m_j = \min(x_j, y_j)`, :math:`M_j = + \max(x_j, y_j)` and :math:`d_j = M_j - m_j`. Then :math:`z_j` is sampled + uniformly from :math:`[m_j - \alpha d_j, M_j + \alpha d_j]`, clipped to + the bounds of the variable. Otherwise, the offspring is a copy of the + first parent :math:`x`. + + Parameters + ---------- + parents : + ${parents} Every value must be within the bounds. + lower : + ${lower_bound} + upper : + ${upper_bound} + alpha : + Fraction of the distance between the parents by which the sampling + interval extends beyond them on each side. With ``alpha=0``, the + offspring lie between their parents. + pair_prob : + ${pair_prob} + seed : + ${random_seed} + + Returns + ------- + Array with :math:`\lceil n/2 \rceil` rows, where :math:`n` is the + number of rows of ``parents``. Row ``i`` is the offspring of rows + ``2i`` and ``2i + 1`` of ``parents``. + + See Also + -------- + local_crossover : Offspring inside the box spanned by the parents. + + Notes + ----- + BLX-:math:`\alpha` was proposed by :footcite:t:`EshSch1992`, who recommend + :math:`\alpha = 0.5` because it makes an offspring as likely to lie outside + its parents as between them. Each variable is sampled independently, so + the offspring fill the whole (clipped) box around the parents. Eshelman + and Schaffer do not handle bounds, and they produce two offspring per pair. + + To produce more than one offspring per pair, call the function again. If + ``seed`` is an integer, every call returns the same offspring. Pass the + same :class:`numpy.random.Generator` to every call to get different + offspring each time, as in the example below. + + References + ---------- + .. footbibliography:: + + Examples + -------- + >>> parents = np.array([[0.2, 0.2], [0.6, 0.8]]) + >>> rng = np.random.default_rng(42) + >>> moocore.blend_alpha_crossover(parents, lower=0, upper=1, seed=rng) + array([[0.4174605 , 0.81353203]]) + >>> moocore.blend_alpha_crossover(parents, lower=0, upper=1, seed=rng) + array([[0.03527987, 0.73987298]]) + + """ + parents = _as_matrix(parents, "parents", min_rows=2) + lower, upper = _bounds(lower, upper, parents.shape[1]) + _check_within_bounds(parents, lower, upper) + if not alpha >= 0: + raise ValueError("'alpha' must be non-negative") + pair_prob = _check_probability(pair_prob, "pair_prob") + return _call( + lib.blend_alpha_crossover, + parents, + (len(parents) + 1) // 2, + ffi.from_buffer("double []", lower), + ffi.from_buffer("double []", upper), + float(alpha), + pair_prob, + _get_seed_for_c(seed), + ) + + +@DocSubstitute() +def single_arithmetic_crossover( + parents: ArrayLike, + /, + *, + pair_prob: float = 1.0, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Single arithmetic crossover. + + Each crossover operation on a pair of parents :math:`x, y` produces two + offspring. With probability ``pair_prob``, one variable :math:`k` is + chosen uniformly at random and both offspring get :math:`(x_k + y_k) / 2` + in that variable. All other variables are copied from the parents. + Otherwise, the offspring are copies of their parents. + + Parameters + ---------- + parents : + ${parents} + pair_prob : + ${pair_prob} + seed : + ${random_seed} + + Returns + ------- + ${pair_offspring} + + See Also + -------- + local_crossover : Recombine every variable with a random weight. + + Notes + ----- + :footcite:t:`EibSmi2003` describe single arithmetic recombination with a + weight :math:`\alpha`, where one offspring gets :math:`\alpha y_k + (1 - + \alpha) x_k` and the other gets :math:`\alpha x_k + (1 - \alpha) y_k`. + This function is their example with :math:`\alpha = 1/2`, where both + offspring get the average. Picek, Jakobovic and Golub (2013) compare it + with other real-coded crossover operators. + + References + ---------- + .. footbibliography:: + + S. Picek, D. Jakobovic and M. Golub (2013). On the recombination operator + in the real-coded genetic algorithms. *2013 IEEE Congress on Evolutionary + Computation*, pages 3103--3110. https://doi.org/10.1109/CEC.2013.6557948 + + Examples + -------- + >>> parents = np.array([[0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 1.0]]) + >>> moocore.single_arithmetic_crossover(parents, seed=42) + array([[0. , 0. , 0. , 0.5], + [1. , 1. , 1. , 0.5]]) + + """ + parents = _as_matrix(parents, "parents", min_rows=2) + pair_prob = _check_probability(pair_prob, "pair_prob") + return _call( + lib.single_arithmetic_crossover, + parents, + len(parents), + pair_prob, + _get_seed_for_c(seed), + ) + + +@DocSubstitute() +def local_crossover( + parents: ArrayLike, + /, + *, + pair_prob: float = 1.0, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Local crossover. + + Each crossover operation on a pair of parents :math:`x, y` produces two + offspring. With probability ``pair_prob``, they are + + .. math:: + u_j = a_j x_j + (1 - a_j) y_j, \qquad v_j = (1 - a_j) x_j + a_j y_j, + + where each weight :math:`a_j` is drawn uniformly from :math:`[0, 1)`. The + offspring lie in the box spanned by their parents, and :math:`u + v = x + + y`. Otherwise, the offspring are copies of their parents. + + Parameters + ---------- + parents : + ${parents} + pair_prob : + ${pair_prob} + seed : + ${random_seed} + + Returns + ------- + ${pair_offspring} + + See Also + -------- + single_arithmetic_crossover : Average a single variable. + blend_alpha_crossover : Sample around the parents. + + Notes + ----- + This operator is described by Dumitrescu et al. (2000) and compared with + other real-coded crossover operators by Picek, Jakobovic and Golub (2013). + + References + ---------- + D. Dumitrescu, B. Lazzerini, L. C. Jain and A. Dumitrescu (2000). + *Evolutionary Computation*. CRC Press. https://doi.org/10.1201/9781482273960 + + S. Picek, D. Jakobovic and M. Golub (2013). On the recombination operator + in the real-coded genetic algorithms. *2013 IEEE Congress on Evolutionary + Computation*, pages 3103--3110. https://doi.org/10.1109/CEC.2013.6557948 + + Examples + -------- + >>> parents = np.array([[0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]) + >>> moocore.local_crossover(parents, seed=42) + array([[0.04928569, 0.26800606, 0.40134152], + [0.95071431, 0.73199394, 0.59865848]]) + + """ + parents = _as_matrix(parents, "parents", min_rows=2) + pair_prob = _check_probability(pair_prob, "pair_prob") + return _call( + lib.local_crossover, + parents, + len(parents), + pair_prob, + _get_seed_for_c(seed), + ) + + +@DocSubstitute() +def bounded_exponential_crossover( + parents: ArrayLike, + /, + *, + lower: ArrayLike | None = None, + upper: ArrayLike | None = None, + scale: float = 0.1, + swap_prob: float = 0.0, + pair_prob: float = 1.0, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Bounded exponential crossover (BEX). + + Each crossover operation on a pair of parents :math:`x, y` produces two + offspring :math:`u, v`. With probability ``pair_prob``, for each + variable :math:`j` with :math:`d_j = |x_j - y_j| > 0`, + + .. math:: + u_j = x_j + \beta_j^x d_j, \qquad v_j = y_j + \beta_j^y d_j, + + where :math:`\beta_j^x` is negative or positive with probability 1/2 each. + On each side, it follows an exponential distribution with scale + :math:`\lambda` (``scale``), truncated at the bound on that side, so that + :math:`u_j` stays within the bounds of the variable. Towards the lower + bound :math:`l_j`, :math:`\beta = \lambda \ln(e_l + r(1 - e_l))` with + :math:`e_l = \exp((l_j - x_j) / (\lambda d_j))`. Towards the upper bound + :math:`h_j`, :math:`\beta = -\lambda \ln(1 - r(1 - e_h))` with :math:`e_h = + \exp((x_j - h_j) / (\lambda d_j))`. The uniform random number :math:`r` + and the direction are shared by both offspring. Variables where the + parents are equal are copied. After that, each variable is swapped between + the two offspring with probability ``swap_prob``, as in + :func:`uniform_crossover`. + + Otherwise, with probability ``1 - pair_prob``, the offspring are copies of + their parents. + + Parameters + ---------- + parents : + ${parents} Every value must be within the bounds. + lower : + ${lower_bound} + upper : + ${upper_bound} + scale : + Scale :math:`\lambda > 0` of the exponential distribution. Smaller + values keep the offspring closer to their parents. + swap_prob : + ${swap_prob} + pair_prob : + ${pair_prob} + seed : + ${random_seed} + + Returns + ------- + ${pair_offspring} + + Notes + ----- + BEX was proposed by Thakur, Meghwani and Jalota (2014) as a bounded + version of the Laplace crossover of Deep and Thakur (2007). Without + bounds, :math:`\beta` follows a Laplace distribution. A parent that sits + on a bound is copied with probability 1/2, as in the paper. + + The swap of variables (``swap_prob``) is not part of the original operator. + Thakur, Meghwani and Jalota tuned :math:`\lambda = 0.65` and a crossover + probability of 0.5 for their genetic algorithm. The defaults here keep + the offspring closer to their parents and cross over every pair. + + References + ---------- + M. Thakur, S. S. Meghwani and H. Jalota (2014). A modified real coded + genetic algorithm for constrained optimization. *Applied Mathematics and + Computation*, 235:292--317. https://doi.org/10.1016/j.amc.2014.02.093 + + K. Deep and M. Thakur (2007). A new crossover operator for real coded + genetic algorithms. *Applied Mathematics and Computation*, + 188(1):895--911. https://doi.org/10.1016/j.amc.2006.10.047 + + Examples + -------- + >>> parents = np.array([[0.2, 0.2], [0.6, 0.8]]) + >>> moocore.bounded_exponential_crossover(parents, lower=0, upper=1, seed=42) + array([[0.32040486, 0.14833417], + [0.72036984, 0.74522359]]) + + """ + parents = _as_matrix(parents, "parents", min_rows=2) + lower, upper = _bounds(lower, upper, parents.shape[1]) + _check_within_bounds(parents, lower, upper) + scale = _check_positive(scale, "scale") + swap_prob = _check_probability(swap_prob, "swap_prob") + pair_prob = _check_probability(pair_prob, "pair_prob") + return _call( + lib.bounded_exponential_crossover, + parents, + len(parents), + ffi.from_buffer("double []", lower), + ffi.from_buffer("double []", upper), + scale, + swap_prob, + pair_prob, + _get_seed_for_c(seed), + ) + + +@DocSubstitute() +def de_rand_1_bin_crossover( + population: ArrayLike, + /, + *, + targets: ArrayLike | None = None, + lower: ArrayLike | None = None, + upper: ArrayLike | None = None, + scaling_factor: float = 0.5, + crossover_rate: float = 0.9, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Differential evolution crossover (DE/rand/1/bin). + + For each target :math:`x_i`, three rows :math:`x_{r_1}, x_{r_2}, x_{r_3}` + are drawn from ``population``, distinct from each other and from the + target. The offspring takes each variable :math:`j` from the mutant + :math:`x_{r_1} + F (x_{r_2} - x_{r_3})` with probability :math:`CR`, and + from the target otherwise. One variable chosen at random always comes + from the mutant. Finally, the offspring is clipped to the bounds. + + Parameters + ---------- + population : + ${population} It must have at least 4 rows. + targets : + ${targets} + lower : + ${lower_bound} + upper : + ${upper_bound} + scaling_factor : + ${scaling_factor} + crossover_rate : + ${crossover_rate} + seed : + ${random_seed} + + Returns + ------- + Array with one row per target, where row ``i`` is the offspring of + target ``targets[i]``. + + See Also + -------- + de_current_1_crossover : The variant used by MOEA/D-DE. + + Notes + ----- + This is the DE/rand/1/bin scheme of :footcite:t:`StoPri1997:de`. They + recommend :math:`F = 0.5` as a first choice. They call :math:`CR = 0.1` a + good first choice, but suggest first trying :math:`CR = 0.9` or 1.0 to see + whether a quick solution is possible, which gives the default here. + Multi-objective DE algorithms such as GDE3 :footcite:p:`KukLam2005:gde3` + and DEMO :footcite:p:`RobFil05demo` use this variant + :footcite:p:`BezLopStu2015emode`. + + Storn and Price do not specify how to handle values outside the bounds. + Clipping them, as done here, is a common choice. + + References + ---------- + .. footbibliography:: + + Examples + -------- + >>> population = np.array( + ... [[0.1, 0.1], [0.4, 0.9], [0.8, 0.3], [0.5, 0.5], [0.9, 0.9]] + ... ) + >>> moocore.de_rand_1_bin_crossover( + ... population, targets=[0, 3], lower=0, upper=1, seed=42 + ... ) + array([[1. , 0.5], + [0. , 0.1]]) + + """ + population = _as_matrix(population, "population", min_rows=4) + targets = _targets(targets, len(population)) + lower, upper = _bounds(lower, upper, population.shape[1]) + scaling_factor = _check_positive(scaling_factor, "scaling_factor") + crossover_rate = _check_probability(crossover_rate, "crossover_rate") + return _call_with_targets( + lib.de_rand_1_bin_crossover, + population, + targets, + lower, + upper, + scaling_factor, + crossover_rate, + _get_seed_for_c(seed), + ) + + +@DocSubstitute() +def de_current_1_crossover( + population: ArrayLike, + /, + *, + targets: ArrayLike | None = None, + lower: ArrayLike | None = None, + upper: ArrayLike | None = None, + scaling_factor: float = 0.5, + crossover_rate: float = 1.0, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Differential evolution crossover used by MOEA/D-DE (DE/current/1). + + For each target :math:`x_i`, two distinct rows :math:`x_{r_1}, x_{r_2}` + are drawn from ``population``. Either of them may be the target itself. + The offspring takes each variable :math:`j` from the mutant :math:`x_i + F + (x_{r_1} - x_{r_2})` with probability :math:`CR`, and from the target + otherwise. Finally, the offspring is clipped to the bounds. + + Unlike :func:`de_rand_1_bin_crossover`, the base vector is the target and + no variable is forced to come from the mutant, so with :math:`CR = 0` the + offspring is a copy of the target. + + Parameters + ---------- + population : + ${population} In MOEA/D-DE, this is the neighbourhood of the + subproblem or, with a small probability, the whole population, and + ``targets`` gives the row of the subproblem's current solution. + targets : + ${targets} + lower : + ${lower_bound} + upper : + ${upper_bound} + scaling_factor : + ${scaling_factor} + crossover_rate : + ${crossover_rate} + seed : + ${random_seed} + + Returns + ------- + Array with one row per target, where row ``i`` is the offspring of + target ``targets[i]``. + + See Also + -------- + de_rand_1_bin_crossover : The classical DE/rand/1/bin scheme. + + Notes + ----- + This is the DE operator of MOEA/D-DE :footcite:p:`LiZha2009:moead-de`, + which uses :math:`CR = 1` and :math:`F = 0.5`. With :math:`CR = 1`, the + operator is invariant to rotations of the decision space. MOEA/D-DE then + applies polynomial mutation, which is not part of this function. + + Li and Zhang do not name the mutation step. Tanabe and Ishibuchi (2019) + call it the current/1 strategy, a term introduced by Gong et al. (2017), + and find it in the source code of MOEA/D-DE, MOEA/D-DRA and MOEA/D-STM. + + Li and Zhang do not require the two donors to differ. This function draws + two distinct donors, either of which may be the target, as their C++ code + does according to Tanabe and Ishibuchi (2019). + + Li and Zhang reset values outside the bounds to a random value within + them. Their own implementation clips them instead, which is what this + function does, and which Tanabe and Ishibuchi (2019) found to work best + with this operator. Tanabe and Ishibuchi also explain why some + implementations of MOEA/D-DE use DE/rand/1 instead. + + References + ---------- + .. footbibliography:: + + R. Tanabe and H. Ishibuchi (2019). Review and analysis of three components + of the differential evolution mutation operator in MOEA/D-DE. *Soft + Computing*, 23(23):12843--12857. + https://doi.org/10.1007/s00500-019-03842-6 + + W. Gong, Y. Wang, Z. Cai and S. Yang (2017). A weighted biobjective + transformation technique for locating multiple optimal solutions of + nonlinear equation systems. *IEEE Transactions on Evolutionary + Computation*, 21(5):697--713. https://doi.org/10.1109/TEVC.2017.2670779 + + Examples + -------- + >>> neighbourhood = np.array([[0.5, 0.5], [0.4, 0.9], [0.8, 0.3]]) + >>> moocore.de_current_1_crossover( + ... neighbourhood, targets=[0], lower=0, upper=1, seed=42 + ... ) + array([[0.3, 0.8]]) + + """ + population = _as_matrix(population, "population", min_rows=2) + targets = _targets(targets, len(population)) + lower, upper = _bounds(lower, upper, population.shape[1]) + scaling_factor = _check_positive(scaling_factor, "scaling_factor") + crossover_rate = _check_probability(crossover_rate, "crossover_rate") + return _call_with_targets( + lib.de_current_1_crossover, + population, + targets, + lower, + upper, + scaling_factor, + crossover_rate, + _get_seed_for_c(seed), + ) + + +@DocSubstitute() +def parent_centric_crossover( + population: ArrayLike, + /, + *, + targets: ArrayLike | None = None, + lower: ArrayLike | None = None, + upper: ArrayLike | None = None, + sigma_zeta: float = 0.1, + sigma_eta: float = 0.1, + seed: int | np.random.Generator | None = None, +) -> np.ndarray: + r"""Parent-centric crossover (PCX). + + Each target :math:`x` is an index parent. Two more parents :math:`a, b` + are drawn from ``population``, distinct from each other and from the + target. Let :math:`d = x - g`, where :math:`g` is the centroid of the + three parents, and let :math:`\bar{D}` be the mean distance from :math:`a` + and :math:`b` to the line through :math:`x` along :math:`d`. The + offspring is + + .. math:: + x + \zeta d + \bar{D} \eta, + + where :math:`\zeta \sim N(0, \sigma_\zeta^2)` and :math:`\eta` is a + vector of independent :math:`N(0, \sigma_\eta^2)` values projected onto + the subspace orthogonal to :math:`d`. Finally, the offspring is clipped to + the bounds. + + Parameters + ---------- + population : + ${population} It must have at least 3 rows. + targets : + ${targets} + lower : + ${lower_bound} + upper : + ${upper_bound} + sigma_zeta : + Standard deviation :math:`\sigma_\zeta > 0` of the step along + :math:`d`. + sigma_eta : + Standard deviation :math:`\sigma_\eta > 0` of the steps orthogonal to + :math:`d`. + seed : + ${random_seed} + + Returns + ------- + Array with one row per target, where row ``i`` is the offspring of + index parent ``targets[i]``. + + Notes + ----- + PCX was proposed by Deb, Anand and Joshi (2002), who use three parents and + :math:`\sigma_\zeta = \sigma_\eta = 0.1`. Their equation (2) writes the + step orthogonal to :math:`d` as a sum of normal steps along :math:`\mu - 1` + orthonormal vectors, where :math:`\mu` is the number of parents. Like + pymoo, which bases its implementation on the authors' C code, this + function projects an isotropic normal vector onto the whole subspace + orthogonal to :math:`d`, so the step covers every direction orthogonal to + :math:`d`, not only :math:`\mu - 1` of them. If :math:`x` is (almost) the + centroid, the projection is skipped. + + References + ---------- + K. Deb, A. Anand and D. Joshi (2002). A computationally efficient + evolutionary algorithm for real-parameter optimization. *Evolutionary + Computation*, 10(4):371--395. https://doi.org/10.1162/106365602760972767 + + Examples + -------- + >>> population = np.array([[0.2, 0.2], [0.4, 0.9], [0.8, 0.3], [0.5, 0.5]]) + >>> moocore.parent_centric_crossover( + ... population, targets=[0, 3], lower=0, upper=1, seed=42 + ... ) + array([[0.18832072, 0.19642879], + [0.50136187, 0.44711562]]) + + """ + population = _as_matrix(population, "population", min_rows=3) + targets = _targets(targets, len(population)) + lower, upper = _bounds(lower, upper, population.shape[1]) + sigma_zeta = _check_positive(sigma_zeta, "sigma_zeta") + sigma_eta = _check_positive(sigma_eta, "sigma_eta") + return _call_with_targets( + lib.parent_centric_crossover, + population, + targets, + lower, + upper, + sigma_zeta, + sigma_eta, + _get_seed_for_c(seed), + ) diff --git a/python/src/moocore/_docsubstitute.py b/python/src/moocore/_docsubstitute.py index 4e1f142d6..2a35a0d93 100644 --- a/python/src/moocore/_docsubstitute.py +++ b/python/src/moocore/_docsubstitute.py @@ -27,6 +27,30 @@ # ------------------------------------------------------------------------- "whv_rectangles": """Weighted rectangles that will bias the computation of the hypervolume. Maybe generated by :func:`eafdiff()` with ``rectangles=True`` or by :func:`choose_eafdiff()`.""", + # ------------------------------------------------------------------------- + "parents": """2D array of numerical values, where each row is a decision vector. + Rows ``2i`` and ``2i + 1`` are crossed over. If the number of rows is odd, the last row is crossed over with the first one.""", + # ------------------------------------------------------------------------- + "pair_offspring": """Array with the same shape as ``parents``, where rows ``2i`` and ``2i + 1`` are the two offspring of rows ``2i`` and ``2i + 1`` of ``parents``. + If the number of rows is odd, the second offspring of the last pair is dropped, so the number of rows does not change.""", + # ------------------------------------------------------------------------- + "population": """2D array of numerical values, where each row is a decision vector.""", + # ------------------------------------------------------------------------- + "targets": """1D vector of row indices of ``population``. Each index produces one offspring. ``None`` uses every row.""", + # ------------------------------------------------------------------------- + "lower_bound": """Lower bounds of the decision variables. Either a single value, which will be used for all variables, + or one value per column. ``None`` means no lower bound.""", + # ------------------------------------------------------------------------- + "upper_bound": """Upper bounds of the decision variables. Either a single value, which will be used for all variables, + or one value per column. ``None`` means no upper bound.""", + # ------------------------------------------------------------------------- + "swap_prob": """Probability of swapping each variable between the two offspring.""", + # ------------------------------------------------------------------------- + "pair_prob": """Probability of applying the crossover operation to each pair of parents.""", + # ------------------------------------------------------------------------- + "scaling_factor": """Scaling factor :math:`F > 0` of the difference vector.""", + # ------------------------------------------------------------------------- + "crossover_rate": """Crossover rate :math:`CR`, the probability of taking each variable from the mutant.""", } diff --git a/python/src/moocore/_ffi_build.py b/python/src/moocore/_ffi_build.py index 87e32fff0..1803bbeec 100644 --- a/python/src/moocore/_ffi_build.py +++ b/python/src/moocore/_ffi_build.py @@ -30,9 +30,11 @@ #include "whv.h" #include "whv_hype.h" #include "hvapprox.h" +#include "crossover.h" """ sources = [ "avl.c", + "crossover.c", "eaf.c", "eaf3d.c", "eafdiff.c", diff --git a/python/src/moocore/libmoocore.h b/python/src/moocore/libmoocore.h index 5f0c0a45e..dfa783685 100644 --- a/python/src/moocore/libmoocore.h +++ b/python/src/moocore/libmoocore.h @@ -85,3 +85,14 @@ hype_sample_dist * hype_dist_gaussian_new(const double *mu, unsigned long int se void hype_dist_free(hype_sample_dist * d); double whv_hype_estimate(const double *points, size_t n, const double *ideal, const double *ref, hype_sample_dist * dist, size_t nsamples); */ + +// crossover.h +void uniform_crossover(double * restrict offspring, const double * restrict parents, size_t n, size_t nvar, double swap_prob, double pair_prob, uint32_t seed); +void single_point_crossover(double * restrict offspring, const double * restrict parents, size_t n, size_t nvar, double pair_prob, uint32_t seed); +void blend_alpha_crossover(double * restrict offspring, const double * restrict parents, size_t n, size_t nvar, const double * restrict lower, const double * restrict upper, double alpha, double pair_prob, uint32_t seed); +void single_arithmetic_crossover(double * restrict offspring, const double * restrict parents, size_t n, size_t nvar, double pair_prob, uint32_t seed); +void local_crossover(double * restrict offspring, const double * restrict parents, size_t n, size_t nvar, double pair_prob, uint32_t seed); +void bounded_exponential_crossover(double * restrict offspring, const double * restrict parents, size_t n, size_t nvar, const double * restrict lower, const double * restrict upper, double scale, double swap_prob, double pair_prob, uint32_t seed); +void de_rand_1_bin_crossover(double * restrict offspring, const double * restrict population, size_t n, size_t nvar, const size_t * restrict targets, size_t ntargets, const double * restrict lower, const double * restrict upper, double scaling_factor, double crossover_rate, uint32_t seed); +void de_current_1_crossover(double * restrict offspring, const double * restrict population, size_t n, size_t nvar, const size_t * restrict targets, size_t ntargets, const double * restrict lower, const double * restrict upper, double scaling_factor, double crossover_rate, uint32_t seed); +void parent_centric_crossover(double * restrict offspring, const double * restrict population, size_t n, size_t nvar, const size_t * restrict targets, size_t ntargets, const double * restrict lower, const double * restrict upper, double sigma_zeta, double sigma_eta, uint32_t seed); diff --git a/python/tests/test_crossover.py b/python/tests/test_crossover.py new file mode 100644 index 000000000..cb93be758 --- /dev/null +++ b/python/tests/test_crossover.py @@ -0,0 +1,425 @@ +# ruff: noqa: D100, D103 +import math + +import numpy as np +import pytest +from numpy.testing import assert_allclose, assert_array_equal + +import moocore + +LOWER, UPPER = 0.0, 10.0 + +TWO_OFFSPRING = { + "uniform": moocore.uniform_crossover, + "single_point": moocore.single_point_crossover, + "single_arithmetic": moocore.single_arithmetic_crossover, + "local": moocore.local_crossover, + "bounded_exponential": moocore.bounded_exponential_crossover, +} +PAIR_OPERATORS = {**TWO_OFFSPRING, "blend_alpha": moocore.blend_alpha_crossover} +POPULATION_OPERATORS = { + "de_rand_1_bin": moocore.de_rand_1_bin_crossover, + "de_current_1": moocore.de_current_1_crossover, + "parent_centric": moocore.parent_centric_crossover, +} +ALL_OPERATORS = {**PAIR_OPERATORS, **POPULATION_OPERATORS} +BOUNDED = { + "blend_alpha", + "bounded_exponential", + *POPULATION_OPERATORS, +} + + +def random_matrix(n, nvar=4, seed=0): + return np.random.default_rng(seed).uniform(LOWER, UPPER, size=(n, nvar)) + + +def pairs(x): + return x[0::2], x[1::2] + + +@pytest.mark.parametrize("name", TWO_OFFSPRING) +@pytest.mark.parametrize("n", [2, 3, 4, 5, 10]) +def test_two_offspring_per_pair_keep_row_count(name, n): + x = random_matrix(n) + assert TWO_OFFSPRING[name](x, seed=1).shape == x.shape + + +@pytest.mark.parametrize("n", [2, 3, 4, 5, 10]) +def test_blend_alpha_crossover_one_offspring_per_pair(n): + x = random_matrix(n) + out = moocore.blend_alpha_crossover(x, seed=1) + assert out.shape == (math.ceil(n / 2), x.shape[1]) + + +@pytest.mark.parametrize("name", POPULATION_OPERATORS) +def test_population_operators_one_offspring_per_target(name): + fun = POPULATION_OPERATORS[name] + x = random_matrix(6) + assert fun(x, seed=1).shape == x.shape + assert fun(x, targets=[4, 0, 4], seed=1).shape == (3, x.shape[1]) + assert fun(x, targets=2, seed=1).shape == (1, x.shape[1]) + assert fun(x, targets=[], seed=1).shape == (0, x.shape[1]) + + +@pytest.mark.parametrize("n", [3, 5, 7]) +def test_odd_number_of_parents_pairs_last_with_first(n): + x = random_matrix(n) + assert_array_equal(moocore.uniform_crossover(x, swap_prob=0, seed=1), x) + # Swapping every variable gives each offspring the other parent. + out = moocore.uniform_crossover(x, swap_prob=1, seed=1) + assert_array_equal(out[0], x[1]) + assert_array_equal(out[-1], x[0]) + + out = moocore.blend_alpha_crossover(x, alpha=0, seed=1) + assert np.all(out[-1] >= np.minimum(x[-1], x[0])) + assert np.all(out[-1] <= np.maximum(x[-1], x[0])) + + +@pytest.mark.parametrize("name", PAIR_OPERATORS) +@pytest.mark.parametrize("n", [3, 4, 5]) +def test_identical_parents_reproduce_themselves(name, n): + # The row has distinct values, so a transposed output would not match. + row = np.arange(1.0, 6.0) + x = np.tile(row, (n, 1)) + out = PAIR_OPERATORS[name](x, seed=1) + assert_allclose(out, np.tile(row, (len(out), 1))) + + +@pytest.mark.parametrize("name", ALL_OPERATORS) +def test_seed(name): + fun = ALL_OPERATORS[name] + x = random_matrix(8) + assert_array_equal(fun(x, seed=42), fun(x, seed=42)) + assert_array_equal( + fun(x, seed=np.random.default_rng(42)), + fun(x, seed=np.random.default_rng(42)), + ) + assert not np.array_equal(fun(x, seed=1), fun(x, seed=2)) + # The same Generator gives new offspring on each call. + rng = np.random.default_rng(42) + assert not np.array_equal(fun(x, seed=rng), fun(x, seed=rng)) + # seed=None still works. + assert fun(x).shape[1] == x.shape[1] + + +@pytest.mark.parametrize("name", ALL_OPERATORS) +def test_inputs_are_not_modified(name, immutable_call): + x = random_matrix(6) + kwargs = {} + if name in BOUNDED: + kwargs = {"lower": np.full(4, LOWER), "upper": np.full(4, UPPER)} + if name in POPULATION_OPERATORS: + kwargs["targets"] = np.array([0, 2]) + immutable_call(ALL_OPERATORS[name], x, seed=3, **kwargs) + + +@pytest.mark.parametrize("name", PAIR_OPERATORS) +def test_pair_prob_zero_copies_parents(name): + x = random_matrix(7) + kwargs = {"lower": LOWER, "upper": UPPER} if name in BOUNDED else {} + out = PAIR_OPERATORS[name](x, pair_prob=0, seed=1, **kwargs) + # Blend alpha crossover copies the first parent of each pair. + assert_array_equal(out, x[0::2] if name == "blend_alpha" else x) + + +@pytest.mark.parametrize("name", PAIR_OPERATORS) +def test_pair_prob_rate(name): + x = random_matrix(4000, nvar=10) + kwargs = {"lower": LOWER, "upper": UPPER} if name in BOUNDED else {} + out = PAIR_OPERATORS[name](x, pair_prob=0.3, seed=1, **kwargs) + if name == "blend_alpha": + copied = np.all(out == x[0::2], axis=1) + else: + copied = np.all((out == x).reshape(-1, 2 * x.shape[1]), axis=1) + assert abs(copied.mean() - 0.7) < 0.03 + + +def test_uniform_crossover_only_copies_parent_values(): + x = random_matrix(200, nvar=6) + p1, p2 = pairs(x) + off1, off2 = pairs(moocore.uniform_crossover(x, seed=1)) + from_first = off1 == p1 + assert np.all(from_first | (off1 == p2)) + assert_array_equal(off2, np.where(from_first, p2, p1)) + assert_array_equal(moocore.uniform_crossover(x, pair_prob=0, seed=1), x) + + +@pytest.mark.parametrize("swap_prob", [0.0, 0.25, 0.5, 0.75, 1.0]) +def test_uniform_crossover_swap_rate(swap_prob): + x = np.tile([np.zeros(10), np.ones(10)], (2000, 1)) + out = moocore.uniform_crossover(x, swap_prob=swap_prob, seed=1) + # The first offspring gets a 1 only where it was swapped. + assert abs(out[0::2].mean() - swap_prob) < 0.02 + + +@pytest.mark.parametrize( + "fun", [moocore.uniform_crossover, moocore.single_point_crossover] +) +@pytest.mark.parametrize("dtype", [int, bool, np.int8, np.float32]) +def test_value_copying_operators_keep_dtype(fun, dtype): + x = np.array([[0, 0, 0, 0], [1, 1, 1, 1], [0, 1, 0, 1]], dtype=dtype) + out = fun(x, seed=1) + if np.issubdtype(dtype, np.floating): + assert out.dtype == np.float64 + else: + assert out.dtype == x.dtype + assert np.all((out == 0) | (out == 1)) + + +def test_single_point_crossover_exchanges_tails(): + nvar = 5 + x = np.tile([np.zeros(nvar), np.ones(nvar)], (2000, 1)) + off1, off2 = pairs(moocore.single_point_crossover(x, seed=1)) + cuts = (off1 == 0).sum(axis=1) + assert_array_equal(off1, (np.arange(nvar) >= cuts[:, None]).astype(float)) + assert_array_equal(off2, 1 - off1) + assert set(cuts) == set(range(1, nvar)) + + +@pytest.mark.parametrize("alpha", [0.0, 0.5, 1.0]) +def test_blend_alpha_crossover_interval(alpha): + x = random_matrix(400, nvar=3) + p1, p2 = pairs(x) + out = moocore.blend_alpha_crossover( + x, lower=LOWER, upper=UPPER, alpha=alpha, seed=1 + ) + lo, hi = np.minimum(p1, p2), np.maximum(p1, p2) + span = hi - lo + assert np.all(out >= np.maximum(lo - alpha * span, LOWER) - 1e-12) + assert np.all(out <= np.minimum(hi + alpha * span, UPPER) + 1e-12) + + +def test_blend_alpha_crossover_samples_variables_independently(): + x = np.tile([[2.0, 8.0], [8.0, 2.0]], (5000, 1)) + out = moocore.blend_alpha_crossover(x, alpha=0, seed=1) + corr = np.corrcoef(out[:, 0], out[:, 1])[0, 1] + assert abs(corr) < 0.05 + + +def test_single_arithmetic_crossover_changes_one_variable(): + x = random_matrix(100, nvar=5) + p1, p2 = pairs(x) + off1, off2 = pairs(moocore.single_arithmetic_crossover(x, seed=1)) + changed = off1 != p1 + assert np.all(changed.sum(axis=1) == 1) + assert_array_equal(changed, off2 != p2) + mean = 0.5 * (p1 + p2) + assert_array_equal(off1[changed], mean[changed]) + assert_array_equal(off2[changed], mean[changed]) + assert_array_equal( + moocore.single_arithmetic_crossover(x, pair_prob=0, seed=1), x + ) + + +def test_local_crossover_offspring_lie_between_parents(): + x = random_matrix(100, nvar=5) + p1, p2 = pairs(x) + off1, off2 = pairs(moocore.local_crossover(x, seed=1)) + assert_allclose(off1 + off2, p1 + p2) + lo, hi = np.minimum(p1, p2) - 1e-12, np.maximum(p1, p2) + 1e-12 + for off in (off1, off2): + assert np.all((off >= lo) & (off <= hi)) + + +def test_bounded_exponential_crossover_shared_values(): + x = np.tile([[0.0, 5.0, 10.0, 3.0], [0.0, 5.0, 10.0, 7.0]], (50, 1)) + out = moocore.bounded_exponential_crossover(x, lower=0, upper=10, seed=1) + assert np.all(np.isfinite(out)) + assert_array_equal(out[:, :3], x[:, :3]) + assert np.all((out >= 0) & (out <= 10)) + + +def test_bounded_exponential_crossover_bounds(): + x = random_matrix(1000) + out = moocore.bounded_exponential_crossover( + x, lower=LOWER, upper=UPPER, scale=1.0, seed=1 + ) + assert np.all((out >= LOWER) & (out <= UPPER)) + out = moocore.bounded_exponential_crossover(x, scale=1.0, seed=1) + assert np.all(np.isfinite(out)) + assert np.any((out < LOWER) | (out > UPPER)) + assert_array_equal( + moocore.bounded_exponential_crossover(x, pair_prob=0, seed=1), x + ) + + +@pytest.mark.parametrize( + ("swap_prob", "expected"), [(0.0, 1.0), (0.5, 0.5), (1.0, 0.0)] +) +def test_bounded_exponential_crossover_swap_prob(swap_prob, expected): + x = np.tile([np.full(10, 2.0), np.full(10, 8.0)], (1000, 1)) + out = moocore.bounded_exponential_crossover( + x, lower=0, upper=10, swap_prob=swap_prob, seed=1 + ) + # Fraction of the variables of the first offspring that stay near its + # parent. + off1 = out[0::2] + kept = (np.abs(off1 - 2) < np.abs(off1 - 8)).mean() + assert abs(kept - expected) < 0.02 + + +@pytest.mark.parametrize("order", [1, -1]) +def test_bounded_exponential_crossover_offspring_near_their_parent(order): + pair = np.array([[3.0, 3.0], [7.0, 7.0]])[::order] + x = np.tile(pair, (2000, 1)) + out = moocore.bounded_exponential_crossover(x, lower=0, upper=10, seed=1) + assert_allclose(np.median(out[0::2], axis=0), pair[0], atol=0.1) + assert_allclose(np.median(out[1::2], axis=0), pair[1], atol=0.1) + + +def test_de_rand_1_bin_crossover_forces_one_mutant_variable(): + x = random_matrix(20, nvar=6) + out = moocore.de_rand_1_bin_crossover(x, crossover_rate=0, seed=1) + assert np.all((out != x).sum(axis=1) == 1) + + +def test_de_rand_1_bin_crossover_base_is_another_row(): + x = random_matrix(10, nvar=3) + out = moocore.de_rand_1_bin_crossover( + x, scaling_factor=1e-12, crossover_rate=1, seed=1 + ) + for i, offspring in enumerate(out): + dist = np.abs(x - offspring).max(axis=1) + nearest = np.argmin(dist) + assert nearest != i + assert dist[nearest] < 1e-9 + + +def test_de_current_1_crossover_base_is_the_target(): + x = random_matrix(10) + assert_array_equal( + moocore.de_current_1_crossover(x, crossover_rate=0, seed=1), x + ) + out = moocore.de_current_1_crossover(x, scaling_factor=1e-12, seed=1) + assert_allclose(out, x, atol=1e-9) + + +def test_de_current_1_crossover_two_rows(): + # The donors are the two rows, in either order. + x = random_matrix(2) + step = 0.5 * (x[0] - x[1]) + out = moocore.de_current_1_crossover(x, scaling_factor=0.5, seed=1) + for i in range(2): + assert np.allclose(out[i], x[i] + step) or np.allclose( + out[i], x[i] - step + ) + + +def test_parent_centric_crossover_tiny_sigmas_return_index_parent(): + x = random_matrix(10) + out = moocore.parent_centric_crossover( + x, sigma_zeta=1e-12, sigma_eta=1e-12, seed=1 + ) + assert_allclose(out, x, atol=1e-9) + + +def test_parent_centric_crossover_sigma_zeta_scales_step_along_direction(): + x = random_matrix(10) + step1 = ( + moocore.parent_centric_crossover( + x, sigma_zeta=0.1, sigma_eta=1e-12, seed=1 + ) + - x + ) + step2 = ( + moocore.parent_centric_crossover( + x, sigma_zeta=0.2, sigma_eta=1e-12, seed=1 + ) + - x + ) + assert_allclose(step2, 2 * step1, atol=1e-9) + + +@pytest.mark.parametrize("name", POPULATION_OPERATORS) +def test_population_operators_clip_to_bounds(name): + x = random_matrix(30) + if name == "parent_centric": + kwargs = {"sigma_zeta": 5.0, "sigma_eta": 5.0} + else: + kwargs = {"scaling_factor": 5.0} + out = POPULATION_OPERATORS[name]( + x, lower=LOWER, upper=UPPER, seed=1, **kwargs + ) + assert np.all((out >= LOWER) & (out <= UPPER)) + assert np.any((out == LOWER) | (out == UPPER)) + + +@pytest.mark.parametrize("name", ALL_OPERATORS) +def test_errors_shape(name): + fun = ALL_OPERATORS[name] + with pytest.raises(ValueError, match=r"must be a 2D array"): + fun([1.0, 2.0, 3.0]) + with pytest.raises(ValueError, match=r"must have at least \d rows"): + fun([[1.0, 2.0]]) + + +@pytest.mark.parametrize("name", sorted(BOUNDED)) +def test_errors_bounds(name): + fun = ALL_OPERATORS[name] + x = random_matrix(6) + with pytest.raises(ValueError, match=r"'lower' must not be larger"): + fun(x, lower=1, upper=0) + with pytest.raises(ValueError, match=r"'upper' must have length 4"): + fun(x, upper=[1, 2]) + + +@pytest.mark.parametrize( + ("fun", "kwargs", "match"), + [ + (moocore.uniform_crossover, {"swap_prob": 1.5}, "'swap_prob'"), + (moocore.uniform_crossover, {"pair_prob": -0.1}, "'pair_prob'"), + (moocore.single_arithmetic_crossover, {"pair_prob": 2}, "'pair_prob'"), + (moocore.single_point_crossover, {"pair_prob": 2}, "'pair_prob'"), + (moocore.blend_alpha_crossover, {"pair_prob": -1}, "'pair_prob'"), + (moocore.local_crossover, {"pair_prob": 1.5}, "'pair_prob'"), + (moocore.blend_alpha_crossover, {"alpha": -1}, "'alpha'"), + (moocore.bounded_exponential_crossover, {"scale": 0}, "'scale'"), + (moocore.bounded_exponential_crossover, {"swap_prob": 2}, "swap_prob"), + (moocore.de_rand_1_bin_crossover, {"scaling_factor": 0}, "scaling"), + (moocore.de_rand_1_bin_crossover, {"crossover_rate": 2}, "crossover"), + (moocore.de_current_1_crossover, {"scaling_factor": -1}, "scaling"), + (moocore.de_current_1_crossover, {"crossover_rate": -1}, "crossover"), + (moocore.parent_centric_crossover, {"sigma_zeta": 0}, "'sigma_zeta'"), + (moocore.parent_centric_crossover, {"sigma_eta": 0}, "'sigma_eta'"), + ], +) +def test_errors_parameters(fun, kwargs, match): + with pytest.raises(ValueError, match=match): + fun(random_matrix(6), **kwargs) + + +def test_errors_operator_specific(): + with pytest.raises(ValueError, match=r"at least 2 columns"): + moocore.single_point_crossover([[0], [1]]) + with pytest.raises(ValueError, match=r"at least 4 rows"): + moocore.de_rand_1_bin_crossover(random_matrix(3)) + with pytest.raises(ValueError, match=r"at least 3 rows"): + moocore.parent_centric_crossover(random_matrix(2)) + with pytest.raises(ValueError, match=r"within 'lower' and 'upper'"): + moocore.bounded_exponential_crossover( + [[0.5, 2.0], [0.2, 0.3]], lower=0, upper=1 + ) + with pytest.raises(ValueError, match=r"within 'lower' and 'upper'"): + moocore.bounded_exponential_crossover( + [[np.nan, 0.5], [0.2, 0.3]], lower=0, upper=1 + ) + with pytest.raises(ValueError, match=r"within 'lower' and 'upper'"): + moocore.blend_alpha_crossover( + [[1.5, 2.0], [0.2, 0.3]], lower=0, upper=1 + ) + + +@pytest.mark.parametrize("name", POPULATION_OPERATORS) +@pytest.mark.parametrize( + ("targets", "match"), + [ + ([6], r"within \[0, 5\]"), + ([-1], r"within \[0, 5\]"), + ([[0, 1]], "1D vector"), + ([0.5], "1D vector"), + ], +) +def test_errors_targets(name, targets, match): + with pytest.raises(ValueError, match=match): + POPULATION_OPERATORS[name](random_matrix(6), targets=targets) diff --git a/r/inst/REFERENCES.bib b/r/inst/REFERENCES.bib index b752aab29..1e533f89b 100644 --- a/r/inst/REFERENCES.bib +++ b/r/inst/REFERENCES.bib @@ -343,6 +343,17 @@ @article{LacKlaFon2017box doi = {10.1016/j.cor.2016.06.021} } +@article{LiZha2009:moead-de, + title = {Multiobjective Optimization Problems with Complicated + {Pareto} sets, {MOEA/D} and {NSGA-II}}, + author = {Li, Hui and Zhang, Qingfu }, + journal = {IEEE Transactions on Evolutionary Computation}, + volume = 13, + number = 2, + pages = {284--302}, + year = 2009 +} + @article{LopVerDreDoe2025, author = { Manuel L{\'o}pez-Ib{\'a}{\~n}ez and Diederick Vermetten and Johann Dreo and Carola Doerr }, title = {Using the Empirical Attainment Function for Analyzing @@ -412,6 +423,19 @@ @article{SchKer2025r2v2 doi = {10.1162/evco.a.366} } +@article{StoPri1997:de, + author = {Storn, Rainer and Price, Kenneth}, + title = {Differential Evolution -- A Simple and Efficient Heuristic + for Global Optimization over Continuous Spaces}, + journal = {Journal of Global Optimization}, + year = 1997, + volume = 11, + number = 4, + pages = {341--359}, + annote = {Proposed differential evolution}, + doi = {10.1023/A:1008202821328} +} + @article{WuAza2001metrics, author = {J. Wu and S. Azam}, title = {Metrics for Quality Assessment of a Multiobjective Design @@ -480,6 +504,29 @@ @incollection{AugBadBroZit2009gecco2 pages = {563--570} } +@misc{BezLopStu2015emoDEsupp, + author = { Leonardo C. T. Bezerra and Manuel L{\'o}pez-Ib{\'a}{\~n}ez and Thomas St{\"u}tzle }, + title = {To {DE} or Not to {DE}? {Multi}-objective Differential + Evolution Revisited from a Component-Wise Perspective: {Supplementary} material}, + howpublished = {\url{https://iridia.ulb.ac.be/supp/IridiaSupp2015-001/}}, + year = 2015 +} + +@incollection{BezLopStu2015emode, + booktitle = { Evolutionary Multi-criterion Optimization, EMO 2015 Part {I}}, + address = { Heidelberg, Germany}, + series = {Lecture Notes in Computer Science}, + volume = 9018, + year = 2015, + publisher = {Springer}, + editor = { Ant{\'o}nio Gaspar{-}Cunha and Carlos Henggeler Antunes and Carlos A. {Coello Coello} }, + author = { Leonardo C. T. Bezerra and Manuel L{\'o}pez-Ib{\'a}{\~n}ez and Thomas St{\"u}tzle }, + title = {To {DE} or Not to {DE}? {Multi}-objective Differential + Evolution Revisited from a Component-Wise Perspective}, + pages = {48--63}, + doi = {10.1007/978-3-319-15934-8_4} +} + @incollection{BezLopStu2017emo, address = { Cham, Switzerland}, series = {Lecture Notes in Computer Science}, @@ -550,6 +597,14 @@ @incollection{CoeSie2004igd keywords = {IGD} } +@book{EibSmi2003, + title = {Introduction to Evolutionary Computing}, + author = { Agoston E. Eiben and Smith, James E. }, + publisher = {Springer}, + year = 2003, + isbn = 3540401849 +} + @incollection{EmmFon2011emo, booktitle = { Evolutionary Multi-criterion Optimization, EMO 2011}, address = {Berlin~/ Heidelberg}, @@ -590,6 +645,17 @@ @incollection{EmmFon2011emo pages = {121--135} } +@incollection{EshSch1992, + isbn = {1-55860-263-1}, + year = 1993, + publisher = {Morgan Kaufmann Publishers}, + booktitle = {Foundations of Genetic Algorithms (FOGA)}, + editor = { Darrell Whitley }, + author = { Larry J. Eshelman and J. David Schaffer }, + title = {Real-Coded Genetic Algorithms and Interval-Schemata}, + pages = {187--202} +} + @book{FanWan1994numtheory, author = {Fang, K. T. and Wang, Y.}, title = {Number-Theoretic Methods in Statistics}, @@ -681,6 +747,15 @@ @incollection{Gla2017fast keywords = {archiving} } +@book{Goldberg89, + author = { David E. Goldberg }, + title = {Genetic Algorithms in Search, Optimization and + Machine Learning}, + publisher = {Addison-Wesley}, + address = { Boston, MA}, + year = 1989 +} + @incollection{GruFon2009:emaa, editor = { Thomas Bartz-Beielstein and Marco Chiarandini and Lu{\'i}s Paquete and Mike Preuss }, year = 2010, @@ -744,6 +819,13 @@ @techreport{HanJas1998 annote = {Proposed R2 indicator} } +@book{Hol75, + author = { J. Holland }, + title = {Adaptation in Natural and Artificial Systems}, + publisher = {University of Michigan Press}, + year = 1975 +} + @inproceedings{IshHeSha2019regular, month = jun, year = 2019, @@ -786,6 +868,19 @@ @inproceedings{KnoCor2002cec pages = {711--716} } +@inproceedings{KukLam2005:gde3, + address = {Piscataway, NJ}, + publisher = {IEEE Press}, + month = sep, + year = 2005, + booktitle = {Proceedings of the 2005 Congress on Evolutionary Computation (CEC 2005)}, + key = {IEEE CEC}, + author = {Kukkonen, S. and Lampinen, J.}, + title = {{GDE3}: the third evolution step of generalized differential + evolution}, + pages = {443--450} +} + @incollection{Lop2026hvapprox, location = {San Jos{\'e}, Costa Rica}, address = { New York, NY}, @@ -881,6 +976,19 @@ @misc{Rob2018unreasonable annote = {Proposed the Rphi low-discrepancy sequence} } +@incollection{RobFil05demo, + address = {Berlin~/ Heidelberg}, + publisher = {Springer}, + volume = 3410, + series = {Lecture Notes in Computer Science}, + editor = { Carlos A. {Coello Coello} and Hern{\'a}ndez Aguirre, Arturo and Eckart Zitzler }, + year = 2005, + booktitle = { Evolutionary Multi-criterion Optimization, EMO 2005}, + author = {Tea Robi{\v c} and Bogdan Filipi{\v c}}, + title = {{DEMO}: Differential Evolution for Multiobjective Optimization}, + pages = {520--533} +} + @book{RubMel1998simulation, author = {R. Y. Rubinstein and B. Melamed}, title = {Modern simulation and modeling}, @@ -890,6 +998,17 @@ @book{RubMel1998simulation notes = {Uniform sampling from the simplex.} } +@inproceedings{Syswerda89, + publisher = {Morgan Kaufmann Publishers, San Mateo, CA}, + editor = { J. David Schaffer }, + year = 1989, + booktitle = {Proceedings of the Third International Conference on Genetic Algorithms (ICGA'89)}, + author = { Gilbert Syswerda }, + title = {Uniform Crossover in Genetic Algorithms}, + pages = {2--9}, + keywords = {uniform crossover} +} + @inproceedings{VelLam1998gp, year = 1998, publisher = {Stanford University Bookstore},