From 99d6a8727f8c1bc8628c9a5af025d1a35affe1de Mon Sep 17 00:00:00 2001 From: MLopez-Ibanez <2620021+MLopez-Ibanez@users.noreply.github.com> Date: Thu, 17 Sep 2026 14:14:37 +0100 Subject: [PATCH] python/ Implement generate_sequence() * bibkeys.txt (Gla2017fast): New. * python/src/moocore/_generate.py (generate_sequence): New. * python/doc/source/reference/functions.io.rst: Document it. * python/examples/plot_generate_sequence.py: New. * python/doc/source/whatsnew/index.rst: Document. * python/doc/source/REFERENCES.bib: Regenerate. * r/inst/REFERENCES.bib: Regenerate. --- bibkeys.txt | 1 + python/doc/source/REFERENCES.bib | 28 ++ python/doc/source/reference/functions.io.rst | 1 + python/doc/source/whatsnew/index.rst | 1 + python/examples/plot_generate_sequence.py | 432 +++++++++++++++++++ python/src/moocore/__init__.py | 3 +- python/src/moocore/_generate.py | 369 +++++++++++++++- r/inst/REFERENCES.bib | 28 ++ 8 files changed, 861 insertions(+), 2 deletions(-) create mode 100644 python/examples/plot_generate_sequence.py diff --git a/bibkeys.txt b/bibkeys.txt index 472b9efb9..0a0ff5285 100644 --- a/bibkeys.txt +++ b/bibkeys.txt @@ -18,6 +18,7 @@ EmmFon2011emo FanWan1994numtheory FonGueLopPaq2011emo FonPaqLop06:hypervolume +Gla2017fast GruFon2009:emaa Grunert01 GueFon2017hv4d diff --git a/python/doc/source/REFERENCES.bib b/python/doc/source/REFERENCES.bib index 1d482ef7b..b752aab29 100644 --- a/python/doc/source/REFERENCES.bib +++ b/python/doc/source/REFERENCES.bib @@ -653,6 +653,34 @@ @inproceedings{FonPaqLop06:hypervolume exponent even further.} } +@incollection{Gla2017fast, + address = { Cham, Switzerland}, + series = {Lecture Notes in Computer Science}, + volume = 10173, + booktitle = { Evolutionary Multi-criterion Optimization, EMO 2017}, + publisher = {Springer International Publishing}, + year = 2017, + editor = { Heike Trautmann and G{\"u}nther Rudolph and Kathrin Klamroth and Oliver Sch{\"u}tze and Margaret M. Wiecek and Yaochu Jin and Christian Grimme}, + author = { T. Glasmachers }, + title = {A Fast Incremental {BSP} Tree Archive for Non-dominated + Points}, + pages = {252--266}, + doi = {10.1007/978-3-319-54157-0_18}, + abstract = {Maintaining an archive of all non-dominated points is a + standard task in multi-objective optimization. Sometimes it + is sufficient to store all evaluated points and to obtain the + non-dominated subset in a post-processing step. Alternatively + the non-dominated set can be updated on the fly. While + keeping track of many non-dominated points efficiently is + easy for two objectives, we propose an efficient algorithm + based on a binary space partitioning BSP tree for the general + case of three or more objectives. Our analysis and our + empirical results demonstrate the superiority of the method + over the brute-force baseline method, as well as graceful + scaling to large numbers of objectives.}, + keywords = {archiving} +} + @incollection{GruFon2009:emaa, editor = { Thomas Bartz-Beielstein and Marco Chiarandini and Lu{\'i}s Paquete and Mike Preuss }, year = 2010, diff --git a/python/doc/source/reference/functions.io.rst b/python/doc/source/reference/functions.io.rst index e9f5ae6de..8158c444a 100644 --- a/python/doc/source/reference/functions.io.rst +++ b/python/doc/source/reference/functions.io.rst @@ -30,3 +30,4 @@ Generate data :toctree: generated/ generate_ndset + generate_sequence diff --git a/python/doc/source/whatsnew/index.rst b/python/doc/source/whatsnew/index.rst index ea4fcea5f..87bb9d453 100644 --- a/python/doc/source/whatsnew/index.rst +++ b/python/doc/source/whatsnew/index.rst @@ -11,6 +11,7 @@ Version 0.4.0 - :func:`~moocore.vorob_t` returns a :class:`~typing.NamedTuple` instead of a dictionary. - :func:`~moocore.is_nondominated` is up to 10x faster in some inputs thanks to a customized radixsort implementation. - New shapes ``"cliff-concave"`` and ``"cliff-convex"`` added to :func:`~moocore.generate_ndset`. +- New function :func:`~moocore.generate_sequence` to generate sequences of dominated and nondominated points. Version 0.3.2 (11/07/2026) diff --git a/python/examples/plot_generate_sequence.py b/python/examples/plot_generate_sequence.py new file mode 100644 index 000000000..5248135d3 --- /dev/null +++ b/python/examples/plot_generate_sequence.py @@ -0,0 +1,432 @@ +r""" +Sampling Sequences of Dominated and Nondominated Points +======================================================= + +This example illustrates how to sample sequences of multi-dimensional points with various dominance properties using :func:`~moocore.generate_sequence`. + +First we define a few functions useful for plotting. +""" + +# sphinx_gallery_multi_image = "single" +import moocore +import numpy as np +import matplotlib.pyplot as plt +import plotly.graph_objects as go +from plotly.subplots import make_subplots +from matplotlib.patches import Arc +from matplotlib.colors import to_hex + + +def plot_3d(what, x, title, plotly=False, show_index=False): + """Scatter plot of 3D points.""" + if not plotly: + fig = plt.figure() + ax = fig.add_subplot(projection="3d") + + ranks = moocore.pareto_rank(x) + n_ranks = ranks.max() + 1 + + # One distinct color for each rank + cmap = plt.get_cmap("viridis", n_ranks) + rank_colors = [to_hex(cmap(i)) for i in range(n_ranks)] + + match what: + case "simplex": + # Standard 2-simplex vertices in 3D + x_s, y_s, z_s = np.eye(3, dtype=int) + if plotly: + surface = go.Mesh3d( + x=x_s, + y=y_s, + z=z_s, + i=[0], + j=[1], + k=[2], + color="cyan", + opacity=0.2, + flatshading=True, + name="Simplex", + hoverinfo="none", + ) + else: + ax.plot_trisurf( + x_s, + y_s, + z_s, + triangles=[[0, 1, 2]], + color="cyan", + alpha=0.2, + edgecolor="gray", + ) + + case "concave" | "convex": + # Generate points on the positive orthant of the sphere. + phi = np.linspace(0, np.pi / 2, 50) + theta = np.linspace(0, np.pi / 2, 50) + phi, theta = np.meshgrid(phi, theta) + # Convert spherical to Cartesian coordinates (unit sphere) + x_s = np.sin(phi) * np.cos(theta) + y_s = np.sin(phi) * np.sin(theta) + z_s = np.cos(phi) + + if what == "convex": + x_s = 1 - x_s + y_s = 1 - y_s + z_s = 1 - z_s + + if plotly: + surface = go.Surface( + x=x_s, + y=y_s, + z=z_s, + colorscale=[[0, "cyan"], [1, "cyan"]], + opacity=0.2, + showscale=False, + name="Surface", + ) + else: + ax.plot_surface( + x_s, y_s, z_s, color="cyan", alpha=0.2, edgecolor="gray" + ) + + case _: + raise ValueError(f"Unknown plot type {what}") + + if plotly: + scatter = go.Scatter3d( + x=x[:, 0], + y=x[:, 1], + z=x[:, 2], + mode="markers+text" if show_index else "markers", + text=[str(i + 1) for i in range(len(x))], + textposition="top center", + marker=dict(size=2, color=[rank_colors[r] for r in ranks]), + ) + + if x.max() <= 1: + limits = [0, 1] + else: + limits = [x.min(), x.max()] + layout = go.Layout( + title=title, + scene=dict( + xaxis=dict(title="X", range=limits), + yaxis=dict(title="Y", range=limits), + zaxis=dict(title="Z", range=limits), + # Approx. elev=30, azim=25 + camera=dict(eye=dict(x=1.2, y=1.2, z=0.8)), + ), + margin=dict(l=0, r=0, b=0, t=40), + showlegend=False, + ) + fig = go.Figure(data=[surface, scatter], layout=layout) + else: + ax.scatter( + x[:, 0], + x[:, 1], + x[:, 2], + color="blue", + s=20, + marker="o", + depthshade=False, + ) + ax.set( + xlabel="X", + ylabel="Y", + zlabel="Z", + xlim=(0, 1), + ylim=(0, 1), + zlim=(0, 1), + title=title, + ) + ax.view_init(elev=30, azim=25) + + return fig + + +def plotly_3d(what, x, title, show_index=False): + """Scatter plot of 3D points using plotly.""" + return plot_3d( + what=what, x=x, title=title, plotly=True, show_index=show_index + ) + + +def plotly_3d_side_by_side(fig1, fig2): + """Show two plotly 3D figures side-by-side.""" + fig = make_subplots( + rows=1, + cols=2, + specs=[[{"type": "scene"}, {"type": "scene"}]], + subplot_titles=(fig1.layout.title.text, fig2.layout.title.text), + ) + len1 = len(fig1.data) + fig.add_traces(fig1.data, rows=[1] * len1, cols=[1] * len1) + len2 = len(fig2.data) + fig.add_traces(fig2.data, rows=[1] * len2, cols=[2] * len2) + fig.update_layout( + height=400, title="", margin=dict(l=0, r=0, b=0, t=40), showlegend=False + ) + fig.update_scenes(fig1.layout.scene.to_plotly_json()) + return fig + + +def plot_sequence_2d(points, title, ax, show_index=True, show_circle=False): + """Plot a sequence in 2D""" + ranks = moocore.pareto_rank(points) + n_ranks = ranks.max() + 1 + + ax.scatter( + points[:, 0], + points[:, 1], + c=ranks, + cmap=plt.get_cmap("viridis", n_ranks), + vmin=-0.5, + vmax=n_ranks - 0.5, + ) + + if show_index: + for i, (x, y) in enumerate(points): + ax.text(x, y, str(i + 1), fontsize=10, ha="left", va="bottom") + + if show_circle: + radius = 1 if points.max() <= 1 else 2**31 + ax.add_patch( + Arc( + (0, 0), + 2 * radius, + 2 * radius, + theta1=0, + theta2=90, + fill=False, + linewidth=1.5, + linestyle="--", + ) + ) + + if points.min() >= 0 and points.max() <= 1: + ax.set_xlim(0, 1.05) + ax.set_ylim(0, 1.05) + + ax.set_title(title) + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.grid(True) + + +# %% +# +# Sequences in 2D +# --------------- +# +# First, we plot completely ordered sequences. We use different seeds to get +# different sequences. Points are colored according to their Pareto rank, with +# lower ranks having darker colors. In these two sequences, each point has a +# different color (rank) from the rest. +# +n = 10 +fig, axes = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True) +for ax, method, seed in zip( + axes, ["each_dominates_previous", "each_dominates_next"], [42, 43] +): + points = moocore.generate_sequence(n, 2, method=method, seed=seed) + plot_sequence_2d(points, title=f'method="{method}"', ax=ax) + +fig.tight_layout() +plt.show() + +# %% +# Random sequences sampled in the unit hypercube or the unit hypersphere. +# +n = 10 +for method in ["cube", "sphere"]: + fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=True) + show_circle = method == "sphere" + for ax, ndsort in zip(axes, [1, 0, -1]): + points = moocore.generate_sequence( + n, 2, method=method, ndsort=ndsort, seed=42 + ) + plot_sequence_2d( + points, + title=f'method="{method}", ndsort={ndsort}', + ax=ax, + show_circle=show_circle, + ) + + fig.tight_layout() + plt.show() + +# %% +# Naive sampling within the hypersphere produces a bias towards the origin. +# +n = 200 +method = "sphere" +fig, axes = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True) +points = moocore.generate_sequence(n, 2, method=method, seed=42) +plot_sequence_2d( + points, + title=f'method="{method}"', + ax=axes[0], + show_circle=True, + show_index=False, +) + +# Naive +rng = np.random.default_rng(42) +points = np.abs(rng.normal(size=(n, 2))) +points /= np.linalg.norm(points, axis=1, keepdims=True) +points *= rng.uniform(0, 1, size=(n, 1)) +plot_sequence_2d( + points, title="Naive", ax=axes[1], show_circle=True, show_index=False +) +fig.tight_layout() +plt.show() + + +# %% +# Generate the analytic sequence proposed by :cite:t:`Gla2017fast`. +# + +n = 15 +method = "glas2017" +fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=True) +for ax, c_value in zip(axes, [0.5, 1, 2]): + points = moocore.generate_sequence( + n, 2, method=method, c_value=c_value, seed=42 + ) + plot_sequence_2d( + points, title=f'method="{method}", c_value={c_value}', ax=ax + ) + +fig.tight_layout() +plt.show() + + +# %% +# +# Sequences in integer space +# ----------------------------------------- +# +# We can also generate points in integer space. + +n = 10 +fig, axes = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True) +for ax, method, seed in zip( + axes, ["each_dominates_previous", "each_dominates_next"], [42, 43] +): + points = moocore.generate_sequence( + n, 2, method=method, seed=seed, integer=True + ) + plot_sequence_2d(points, title=f'method="{method}"', ax=ax) + +fig.tight_layout() +plt.show() + +for method in ["cube", "sphere"]: + fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=True) + show_circle = method == "sphere" + for ax, ndsort in zip(axes, [1, 0, -1]): + points = moocore.generate_sequence( + n, 2, method=method, ndsort=ndsort, seed=42, integer=True + ) + plot_sequence_2d( + points, + title=f'method="{method}", ndsort={ndsort}', + ax=ax, + show_circle=show_circle, + ) + + fig.tight_layout() + plt.show() + +n = 15 +method = "glas2017" +fig, axes = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=True) +for ax, c_value in zip(axes, [0.5, 1, 2]): + points = moocore.generate_sequence( + n, 2, method=method, c_value=c_value, seed=42, integer=True + ) + plot_sequence_2d( + points, title=f'method="{method}", c_value={c_value}', ax=ax + ) + +fig.tight_layout() +plt.show() + + +# %% +# +# Sequences in 3D +# --------------- +# +# Completely ordered sequences. We use different seeds to get different +# sequences. Points are colored according to their Pareto rank, with lower +# ranks having darker colors. In these two sequences, each point has a +# different color (rank) from the rest. +# +n = 10 +method = "each_dominates_previous" +fig1 = plotly_3d( + "simplex", + moocore.generate_sequence(n, 3, method, seed=42), + title=f'method="{method}"', + show_index=True, +) +method = "each_dominates_next" +fig2 = plotly_3d( + "simplex", + moocore.generate_sequence(n, 3, method, seed=43), + title=f'method="{method}"', + show_index=True, +) + +plotly_3d_side_by_side(fig1, fig2) + + +# %% +# +# Random sequences sampled in the unit hypercube or the unit hypersphere. If +# you hover over a point, a tooltip shows the coordinates and the index of the +# point in the sequence. +# +n = 100 +method = "cube" +fig1 = plotly_3d( + "concave", + moocore.generate_sequence(n, 3, method, seed=42), + title=f'method="{method}"', +) +method = "sphere" +fig2 = plotly_3d( + "concave", + moocore.generate_sequence(n, 3, method, seed=42), + title=f'method="{method}"', +) + +plotly_3d_side_by_side(fig1, fig2) + + +# %% +# +# Generate the analytic sequence proposed by :cite:t:`Gla2017fast`. +# + +n = 100 +method = "glas2017" +c_value = 0.9 +fig1 = plotly_3d( + "simplex", + moocore.generate_sequence(n, 3, method, seed=42, c_value=c_value), + title=f'method="{method}", c_value={c_value}', +) +c_value = 1.1 +fig2 = plotly_3d( + "simplex", + moocore.generate_sequence(n, 3, method, seed=42, c_value=c_value), + title=f'method="{method}", c_value={c_value}', +) + +plotly_3d_side_by_side(fig1, fig2) + +# %% +# .. rubric:: Related examples +# .. minigallery:: ../../examples/plot_generate.py diff --git a/python/src/moocore/__init__.py b/python/src/moocore/__init__.py index 9662f0d06..d07f559f0 100644 --- a/python/src/moocore/__init__.py +++ b/python/src/moocore/__init__.py @@ -36,7 +36,7 @@ get_dataset_path, ) -from ._generate import generate_ndset +from ._generate import generate_ndset, generate_sequence from importlib.metadata import version as _metadata_version @@ -59,6 +59,7 @@ "filter_dominated", "filter_dominated_within_sets", "generate_ndset", + "generate_sequence", "get_dataset", "get_dataset_path", "hv_approx", diff --git a/python/src/moocore/_generate.py b/python/src/moocore/_generate.py index 29a06c07b..1dbf01dc5 100644 --- a/python/src/moocore/_generate.py +++ b/python/src/moocore/_generate.py @@ -6,7 +6,7 @@ # https://github.com/renaudlr/moo-nondominated-sets/blob/master/wfgHardGenerator.R) from ._docsubstitute import DocSubstitute -from ._moocore import any_dominated +from ._moocore import any_dominated, pareto_rank from ._utils import is_integer_value @@ -44,6 +44,12 @@ def generate_ndset( ------- A numeric matrix of size :math:`n \times d` containing nondominated points. + + See Also + -------- + generate_sequence : Generate a sequence of dominated and nondominated points. + + Notes ----- The available methods are: @@ -239,3 +245,364 @@ def _sample_concave_simplex(n: int, d: int) -> np.ndarray: if not any_dominated(y): return y x *= 2 + + +def sort_by_pareto_rank(points: np.ndarray, *, reverse: bool) -> np.ndarray: + """Return points ordered by final front assignment.""" + ranks = pareto_rank(points) + if reverse: + ranks = -ranks + return points[np.argsort(ranks), :] + + +def _generate_glas2017_sequence( + n_points: int, + dim: int, + rng: np.random.Generator, + *, + c_value: float, + dominated_fraction: float = 0.5, + dominated_offset: float = 1.0, +) -> np.ndarray: + """Generate objective vectors following the analytic sequence model of section 6.1 in :footcite:t:`Gla2017fast`.""" + if not 0.0 <= dominated_fraction <= 1.0: + raise ValueError("'dominated_fraction' must be in [0, 1]") + dominated_remaining = int(n_points * dominated_fraction) + is_dom = np.zeros(n_points, dtype=bool) + for index in range(n_points): + if dominated_remaining == 0: + break + + total_remaining = n_points - index + if total_remaining == dominated_remaining: + # nondominated remaining == 0 + is_dom[index:] = True + break + + probability = c_value * dominated_remaining / total_remaining + # If probability >= 1, then dominated=True + if probability >= 1 or rng.random() < probability: + dominated_remaining -= 1 + is_dom[index] = True + + # Sample N(0, I - 1/dim * 11^T) by sampling standard normals and ... + points = rng.normal(size=(n_points, dim)) + # ... projecting each row onto the subspace orthogonal to the all-ones vector: + points -= points.mean(axis=1, keepdims=True) + # Now shift those points that should be dominated. + indices = np.arange(1, n_points + 1, dtype=float) + offsets = np.zeros(n_points, dtype=float) + offsets[is_dom] = (dominated_offset * n_points) / indices[is_dom] + points += offsets[:, None] + return points + + +def _generate_within_hypercube(n, d, rng): + return rng.uniform(size=(n, d)) + + +def _generate_within_hypersphere(n, d, rng, r_min): + x = np.abs(rng.normal(size=(n, d))) + x /= np.linalg.norm(x, axis=1, keepdims=True) + # These are on the surface, so now we add noise to move them between the + # surface and the origin. We need to transform according to the d-root, + # otherwise the distribution will be biased towards the origin. + u = rng.uniform(r_min**d, 1, size=(n, 1)) + u = u ** (1 / d) + x *= u + return x + + +def _generate_total_dominance_chain( + n: int, d: int, *, rng, dominates_previous: bool +) -> np.ndarray: + """Generate a total dominance chain.""" + x = _generate_within_hypercube(n, d, rng) + x.sort(axis=0) + if dominates_previous: + return x[::-1] + return x + + +@DocSubstitute() +def generate_sequence( + n: int, + d: int, + /, + method: str, + *, + seed: int | np.random.Generator | None = None, + integer: bool = False, + n_rep: int = 0, + ndsort: int = 0, + c_value: float = 0.9, + r_min: float = 0.0, +) -> np.ndarray: + r"""Generate a sequence of ``n`` points of dimension ``d`` with the properties defined by ``method``. + + When ``ndsort=1``, the points are sorted according to their Pareto rank + using :func:`pareto_rank`. With ``ndsort=-1``, the order is reversed. The + ranking assumes that all dimensions are minimised. + + When ``integer=False`` (the default), the points are generated within the + hypercube :math:`(0,1)^d`, except for ``method="glas2017"``, which may + generate points outside this range. These points can be scaled to + another range using :func:`normalise`. + + When ``integer=True``, points are scaled to the non-negative integers in + the range :math:`[0,2^{31}]`, except for ``method="glas2017"``, which may + generate points outside this range. + + + Parameters + ---------- + n : + Number of rows in the output. + d : + Number of columns in the output. + method : + Method used to generate the point sequence. See the Notes below for more details. + seed : + ${random_seed} + integer: + If ``True``, return integer-valued points. + n_rep: + If non-zero, then repeat the generated sequence so that the output will have ``n_rep`` rows. If ``0 < n_rep <= n``, raise ``ValueError``. + ndsort : + Whether points are sorted according to Pareto rank (``1``), sorted in reverse (``-1``) or not sorted at all (``0``). + c_value: + Parameter of ``method="glas2017``. It controls the probability of a + dominated point appearing earlier in the sequence than a nondominated + one, with ``c = 1`` giving equal probability, ``c > 1`` increasing the + probability of dominated points and ``c < 1`` increasing the probability of nondominated ones. + r_min: + Minimum distance to the origin (``method="sphere`` only). + + + Returns + ------- + A numeric matrix of size :math:`n \times d` containing a sequence of points. + + + See Also + -------- + generate_ndset : Generate a nondominated set. + + + Notes + ----- + The available methods are: + + ``'cube'`` + Uniformly samples points within the unit hypercube :math:`(0,1)^d`. + + ``'sphere'`` + Uniformly samples points within the positive orthant of the unit hypersphere. + + Each point :math:`\vec{z} \in (0,1)^d \subset \mathbb{R}^d` is generated + by sampling :math:`d` independent and identically distributed values + :math:`\vec{x}=(x_1,x_2, \dots, x_d)` from the standard normal + distribution, then dividing each value by the l2-norm of the vector, + :math:`z_i = \frac{|x_i|}{\|\vec{x}\|_2}` + :footcite:p:`Muller1959sphere`. The absolute value in the numerator + ensures that points are sampled on the surface of the positive orthant of + the hypersphere. Then each point is moved into the interior of the + hypersphere by sampling :math:`\vec{u} \in \mathbb{R}^d`, with each + component uniformly sampled within the interval :math:`(r_\min^d, 1)`, + and returning :math:`\vec{z}\cdot \sqrt[d]{u}`. The :math:`d`-root + transformation avoids biasing the distribution towards the + origin. Parameter :math:`r_\min` (``r_min``) restricts the minimum + distance to the origin. + + ``'each_dominates_previous'|'each_dominates_next'`` + Each point in the sequence dominates the previous or next one. + + A matrix :math:`n\times d` is sampled uniformly within :math:`(0,1)`. + Then, each column is sorted independently in increasing + (``each_dominates_next'``) or decreasing (``each_dominates_previous'``) + order. Argument ``ndsort`` has no effect for these sequences + because they already sorted. + + ``'glas2017'`` + Generate objective vectors following the analytic sequence model of Section 6.1 in :footcite:t:`Gla2017fast`. + + This model constructs a sequence :math:`\vec{z}^(k) \in \mathbb{R}^d` of + length :math:`N`, where :math:`D=\lfloor fN \rfloor` points are dominated + by another point in the sequence and :math:`N - D` are nondominated, and + :math:`f=0.5` in our implementation. First, points :math:`\vec{x}^{(k)}`, + :math:`\forall k=1,\dots,N`, are generated by sampling :math:`N \times d` + independent values from the standard normal distribution. Second, these + points are projected onto the :math:`(d-1)`-dimensional hyperplane that + satisfies :math:`\{\vec{x}\in\mathbb{R}^d\mid \sum_{i=1}^d x_i = 0\}` by + calculating :math:`\vec{y}^{(k)} = \vec{x}^{(k)} - \bar{x}^{(k)}`, where + :math:`\bar{x} = \frac{1}{d}\sum_{i=1}^d x_i`. Finally, the projected + points are shifted by a constant amount in all dimensions + :math:`\vec{z}^{(k)} = \vec{y}^{(k)} + a^{(k)}`, so that exactly + :math:`D` points become dominated, with + + .. math:: + a^{(k)}=\begin{cases}\frac{\delta N}{k} & \text{if }D_k = 1,\\ + 0 & \text{otherwise}.\end{cases} + + where :math:`D_k \in \{0,1\}` determines whether point :math:`k` is + marked dominated and :math:`\delta = 1` in our implementation. + + Point :math:`k` is marked dominated :math:`(D_k=1)` with probability + :math:`c\frac{n^\text{dom}_k}{n_k}`, where :math:`n^\text{dom}_k = D - + \sum_{i=1}^{k-1}D_k`, that is, the remaining points needed to reach + :math:`D` dominated points in the sequence, and :math:`n_k = N - k + 1`, + that is, the remaining points in the sequence. The parameter :math:`c` + (``c_value``) controls the probability of having dominated points early + in the sequence, with :math:`c > 1` increasing this probability and + :math:`c = 1` giving equal probability to dominated and nondominated + points. + + + The argument ``ndsort=1`` sorts the sequence according to Pareto rank using + :func:`~moocore.pareto_rank`, or in reverse order with ``ndsort=-1``. That + is, earlier points in the sequence will have a lower (or higher in reverse + order) or equal rank than later points. Algorithms that expect points to + have increasing quality should perform worse with ``ndsort=1``, whereas + algorithms that expect new points to be often dominated by previous ones + should perform worse with ``ndsort=-1``. + + + References + ---------- + .. footbibliography:: + + + Examples + -------- + Points within the unit 2D-sphere, i.e., within the circle: + + >>> generate_sequence(5, 2, "sphere", ndsort=0, seed=42) + array([[0.17121976, 0.58436446], + [0.60040891, 0.75251188], + [0.66741195, 0.44545079], + [0.33995168, 0.84094855], + [0.01311259, 0.66576442]]) + + Points within the unit 2D-cube, i.e., unit square, sorted by increasing Pareto rank: + + >>> generate_sequence(5, 2, "cube", ndsort=1, seed=42) + array([[0.77395605, 0.43887844], + [0.09417735, 0.97562235], + [0.12811363, 0.45038594], + [0.85859792, 0.69736803], + [0.7611397 , 0.78606431]]) + + Same points but in different order: + + >>> generate_sequence(5, 2, "cube", ndsort=-1, seed=42) + array([[0.85859792, 0.69736803], + [0.7611397 , 0.78606431], + [0.77395605, 0.43887844], + [0.09417735, 0.97562235], + [0.12811363, 0.45038594]]) + + We can add duplicated points to the sequence by repeating it: + + >>> generate_sequence(5, 2, "cube", ndsort=-1, seed=42, n_rep=11) + array([[0.85859792, 0.69736803], + [0.7611397 , 0.78606431], + [0.77395605, 0.43887844], + [0.09417735, 0.97562235], + [0.12811363, 0.45038594], + [0.85859792, 0.69736803], + [0.7611397 , 0.78606431], + [0.77395605, 0.43887844], + [0.09417735, 0.97562235], + [0.12811363, 0.45038594], + [0.85859792, 0.69736803]]) + + These two sequences are already sorted, so ``ndsort`` is not needed: + + >>> generate_sequence(5, 2, "each_dominates_previous", seed=42) + array([[0.85859792, 0.97562235], + [0.77395605, 0.78606431], + [0.7611397 , 0.69736803], + [0.12811363, 0.45038594], + [0.09417735, 0.43887844]]) + >>> generate_sequence(5, 2, "each_dominates_next", seed=42) + array([[0.09417735, 0.43887844], + [0.12811363, 0.45038594], + [0.7611397 , 0.69736803], + [0.77395605, 0.78606431], + [0.85859792, 0.97562235]]) + + A more complicated sequence, not in the unit cube: + + >>> generate_sequence(5, 2, "glas2017", seed=42) + array([[-0.32442784, 0.32442784], + [ 2.7220415 , 2.2779585 ], + [ 0.41812139, -0.41812139], + [ 0.05080302, -0.05080302], + [ 0.46939475, 1.53060525]]) + + Instead of floating-point values, we can generate an integer matrix: + + >>> generate_sequence(5, 2, "cube", ndsort=1, seed=42, integer=True) + array([[1662057958, 942484272], + [ 202244314, 2095133046], + [ 275121931, 967196436], + [1843824993, 1497586439], + [1634535063, 1688060241]]) + >>> generate_sequence(5, 2, "cube", ndsort=-1, seed=42, integer=True) + array([[1843824993, 1497586439], + [1634535063, 1688060241], + [1662057958, 942484272], + [ 202244314, 2095133046], + [ 275121931, 967196436]]) + >>> generate_sequence(5, 2, "glas2017", seed=42, integer=True) + array([[-696703483, 696703483], + [5845539605, 4891878634], + [ 897908837, -897908837], + [ 109098654, -109098654], + [1008017539, 3286949756]]) + + """ + if seed is None or is_integer_value(seed): + seed = np.random.default_rng(seed) + + match method: + case "cube": + sample = _generate_within_hypercube(n, d, seed) + case "sphere": + sample = _generate_within_hypersphere(n, d, seed, r_min=r_min) + case "each_dominates_previous": + sample = _generate_total_dominance_chain( + n, d, rng=seed, dominates_previous=True + ) + case "each_dominates_next": + sample = _generate_total_dominance_chain( + n, d, rng=seed, dominates_previous=False + ) + case "glas2017": + sample = _generate_glas2017_sequence(n, d, seed, c_value=c_value) + case _: + raise ValueError(f"unknown method={method}") + + if ndsort != 0: + if method in [ + "each_dominates_previous", + "each_dominates_next", + "glas2017", + ]: + raise ValueError( + f'`ndsort != {ndsort}` does not make sense with `method == "{method}"' + ) + sample = sort_by_pareto_rank(sample, reverse=(ndsort < 0)) + + if integer: + # FIXME: Glas2017 is problematic because it is not in [0,1], but + # forcing it may change its properties. + sample *= 2**31 + sample = sample.astype(int) + + if n_rep > 0: + if n_rep <= n: + raise ValueError(f"'n_rep' ({n_rep}) must be larger than 'n' ({n})") + return np.resize(sample, (n_rep, d)) + + return sample diff --git a/r/inst/REFERENCES.bib b/r/inst/REFERENCES.bib index 1d482ef7b..b752aab29 100644 --- a/r/inst/REFERENCES.bib +++ b/r/inst/REFERENCES.bib @@ -653,6 +653,34 @@ @inproceedings{FonPaqLop06:hypervolume exponent even further.} } +@incollection{Gla2017fast, + address = { Cham, Switzerland}, + series = {Lecture Notes in Computer Science}, + volume = 10173, + booktitle = { Evolutionary Multi-criterion Optimization, EMO 2017}, + publisher = {Springer International Publishing}, + year = 2017, + editor = { Heike Trautmann and G{\"u}nther Rudolph and Kathrin Klamroth and Oliver Sch{\"u}tze and Margaret M. Wiecek and Yaochu Jin and Christian Grimme}, + author = { T. Glasmachers }, + title = {A Fast Incremental {BSP} Tree Archive for Non-dominated + Points}, + pages = {252--266}, + doi = {10.1007/978-3-319-54157-0_18}, + abstract = {Maintaining an archive of all non-dominated points is a + standard task in multi-objective optimization. Sometimes it + is sufficient to store all evaluated points and to obtain the + non-dominated subset in a post-processing step. Alternatively + the non-dominated set can be updated on the fly. While + keeping track of many non-dominated points efficiently is + easy for two objectives, we propose an efficient algorithm + based on a binary space partitioning BSP tree for the general + case of three or more objectives. Our analysis and our + empirical results demonstrate the superiority of the method + over the brute-force baseline method, as well as graceful + scaling to large numbers of objectives.}, + keywords = {archiving} +} + @incollection{GruFon2009:emaa, editor = { Thomas Bartz-Beielstein and Marco Chiarandini and Lu{\'i}s Paquete and Mike Preuss }, year = 2010,