from __future__ import annotations
# NOTE: if we ever start using SciPy, we can use
# from scipy.special import gamma_function
import lzma
import numbers
import os
import shutil
import tempfile
from collections.abc import Callable, Sequence
from io import StringIO
from typing import Any, Literal, NamedTuple
import numpy as np
from numpy.typing import ArrayLike # For type hints
from ._docsubstitute import DocSubstitute
## The CFFI library is used to create C bindings.
from ._libmoocore import ffi, lib
from ._utils import (
_get_seed_for_c,
array_1d_of_length_n,
asarray_maybe_copy,
is_integer_value,
np1d_to_double_array,
np1d_to_int_array,
np2d_to_double_array,
unique_nosort,
)
# Maximum number of objectives (columns) supported by the C library. These
# mirror the limits defined in c/config.h: most functions support up to 255
# (the range of dimension_t == uint_fast8_t), while the hypervolume and
# hypervolume-approximation functions use a smaller limit for efficiency.
DIMENSION_MAX = lib.MOOCORE_DIMENSION_MAX
HV_DIMENSION_MAX = lib.MOOCORE_HV_DIMENSION_MAX
HVAPPROX_DIMENSION_MAX = lib.MOOCORE_HVAPPROX_DIMENSION_MAX
[docs]
class ReadDatasetsError(Exception):
"""Custom exception class for an error returned by the :func:`read_datasets()` function.
Parameters
----------
error_code:
Error code returned by ``read_datasets()`` C function.
"""
_error_strings = (
"NO_ERROR",
"READ_INPUT_FILE_EMPTY",
"READ_INPUT_WRONG_INITIAL_DIM",
"ERROR_FOPEN",
"ERROR_CONVERSION",
"ERROR_COLUMNS",
)
def __init__(self, error_code: int):
self.error = error_code
self.message = self._error_strings[abs(error_code)]
super().__init__(self.message)
def _check_dimension_max(dim: int, max_dim: int) -> None:
if dim > max_dim:
raise ValueError(
f"This function supports at most {max_dim} columns, but input has {dim}"
)
[docs]
def read_datasets(filename: str | os.PathLike[str] | StringIO) -> np.ndarray:
"""Read an input dataset file, parsing the file and returning a numpy array.
Parameters
----------
filename:
Filename of the dataset file or :class:`io.StringIO` directly containing the file contents.
If it does not contain an absolute path, the filename is relative to the current working directory.
If the filename has extension ``.xz``, it is decompressed to a temporary file before reading it.
Each line of the file corresponds to one point of one dataset. Different datasets are separated by an empty line.
Returns
-------
An array containing a representation of the data in the file.
The first :math:`n-1` columns contain the numerical data for each of the objectives.
The last column contains an identifier for which set the data is relevant to.
Examples
--------
>>> filename = moocore.get_dataset_path("input1.dat")
>>> moocore.read_datasets(filename) # doctest: +ELLIPSIS
array([[ 8.07559653, 2.40702554, 1. ],
[ 8.66094446, 3.64050144, 1. ],
[ 0.20816431, 4.62275469, 1. ],
...
[ 4.92599726, 2.70492519, 10. ],
[ 1.22234394, 5.68950311, 10. ],
[ 7.99466959, 2.81122537, 10. ],
[ 2.12700289, 2.43114174, 10. ]])
The numpy array represents this data:
+-------------+-------------+------------+
| Objective 1 | Objective 2 | Set Number |
+-------------+-------------+------------+
| 8.07559653 | 2.40702554 | 1. |
+-------------+-------------+------------+
| 8.66094446 | 3.64050144 | 1. |
+-------------+-------------+------------+
| ... | ... | ... |
+-------------+-------------+------------+
| 7.99466959 | 2.81122537 | 10. |
+-------------+-------------+------------+
| 2.12700289 | 2.43114174 | 10. |
+-------------+-------------+------------+
It is also possible to read datasets from a string:
>>> from io import StringIO
>>> fh = StringIO("0.5 0.5\\n\\n1 0\\n0 1\\n\\n0.5 0.5")
>>> moocore.read_datasets(fh)
array([[0.5, 0.5, 1. ],
[1. , 0. , 2. ],
[0. , 1. , 2. ],
[0.5, 0.5, 3. ]])
""" # noqa: D301
if isinstance(filename, os.PathLike):
filename = os.fspath(filename)
elif isinstance(filename, StringIO):
with tempfile.NamedTemporaryFile(mode="wt", delete=False) as fdst:
shutil.copyfileobj(filename, fdst)
# FIXME: Avoid recursion
return read_datasets(fdst.name)
else:
filename = os.path.expanduser(filename)
if not os.path.isfile(filename):
raise FileNotFoundError(f"file '{filename}' not found")
if filename.endswith(".xz"):
with lzma.open(filename, "rb") as fsrc:
with tempfile.NamedTemporaryFile(delete=False) as fdst:
shutil.copyfileobj(fsrc, fdst)
filename = fdst.name
else:
fdst = None
# Encode filename to a binary string
filename = filename.encode("utf-8")
# Create return pointers for function
data_p = ffi.new("double **")
ncols_p = ffi.new("int *")
datasize_p = ffi.new("int *")
err_code = lib.read_datasets(filename, data_p, ncols_p, datasize_p)
if fdst:
os.remove(fdst.name)
if err_code != 0:
raise ReadDatasetsError(err_code)
datasize = datasize_p[0]
# Ensure that Python can free the data allocated by the C library.
data_p = ffi.gc(data_p[0], lib.free, size=datasize)
# Convert 1D numpy array to 2D array with (n obj... , sets) columns
return np.frombuffer(ffi.buffer(data_p, datasize)).reshape(-1, ncols_p[0])
def _parse_maximise(
maximise: bool | Sequence[bool | int], nobj: int
) -> np.ndarray:
"""Convert maximise array or single bool to ndarray format."""
return array_1d_of_length_n(maximise, nobj, name="maximise").astype(bool)
def _parse_maximise_to_bool_array(
maximise: bool | Sequence[bool | int], nobj: int
) -> ffi.CData:
"""Convert maximise array or single bool to C array.
In the C library, boolean arrays are of type boolvec, which is uint8_t.
"""
return ffi.from_buffer(
"uint8_t []", _parse_maximise(maximise, nobj).view(np.uint8)
)
def _all_positive(x: np.ndarray) -> bool:
return x.min() > 0
def _igd_python(
points: ArrayLike,
ref: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
p: int = 1,
) -> float:
points = np.asarray(points, dtype=float)
points_sq = (points * points).sum(axis=1)
ref = np.atleast_2d(np.asarray(ref, dtype=float))
ref_sq = (ref * ref).sum(axis=1)
# Exploit that ||r - a||^2 = ||r||^2 + ||a||^2 - 2 r^T \cdot a
sq_dists = ref_sq[:, None] + points_sq[None, :] - 2 * (ref @ points.T)
# (avoid tiny negative values from numerical cancellation)
np.maximum(sq_dists, 0, out=sq_dists)
res = (sq_dists.min(axis=1) ** (p / 2.0)).mean() ** (1.0 / p)
return float(res)
def _igd_plus_python(
points: ArrayLike,
ref: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
) -> float:
points, points_copied = asarray_maybe_copy(points)
ref, ref_copied = asarray_maybe_copy(ref)
ref = np.atleast_2d(ref)
nobj = points.shape[1]
maximise = _parse_maximise(maximise, nobj)
# FIXME: Do this in C.
if maximise.any():
if not points_copied:
points = points.copy()
points[:, maximise] = -points[:, maximise]
if not ref_copied:
ref = ref.copy()
ref[:, maximise] = -ref[:, maximise]
res = np.sqrt(
[(np.maximum(points - r, 0) ** 2).sum(axis=1).min() for r in ref]
).mean()
return float(res)
def _avg_hausdorff_dist_python(
points: ArrayLike,
ref: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
p: int = 1,
) -> float:
return max(
_igd_python(points, ref, maximise, p),
_igd_python(ref, points, maximise, p),
)
def _epsilon_addi_python(
points: ArrayLike,
ref: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
) -> float:
points, points_copied = asarray_maybe_copy(points)
ref, ref_copied = asarray_maybe_copy(ref)
ref = np.atleast_2d(ref)
nobj = points.shape[1]
maximise = _parse_maximise(maximise, nobj)
# FIXME: Do this in C.
if maximise.any():
if not points_copied:
points = points.copy()
points[:, maximise] = -points[:, maximise]
if not ref_copied:
ref = ref.copy()
ref[:, maximise] = -ref[:, maximise]
res = np.array([np.max(points - r, axis=1).min() for r in ref]).max()
return float(res)
def _epsilon_mult_python(
points: ArrayLike,
ref: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
) -> float:
points, points_copied = asarray_maybe_copy(points)
ref, ref_copied = asarray_maybe_copy(ref)
ref = np.atleast_2d(ref)
nobj = points.shape[1]
maximise = _parse_maximise(maximise, nobj)
# FIXME: Do this in C.
if maximise.any():
if not points_copied:
points = points.copy()
points[:, maximise] = 1.0 / points[:, maximise]
if not ref_copied:
ref = ref.copy()
ref[:, maximise] = 1.0 / ref[:, maximise]
res = np.array([np.max(points / r, axis=1).min() for r in ref]).max()
return float(res)
def _unary_refset_common(
points: ArrayLike,
ref: ArrayLike,
maximise: bool | Sequence[bool | int],
check_all_positive: bool = False,
) -> tuple[
int, ffi.CData, ffi.CData, ffi.CData, ffi.CData, ffi.CData, np.ndarray
]:
# Convert to numpy.array in case the user provides a list. We use
# np.asarray(dtype=float) to convert it to floating-point, otherwise if a
# user inputs something like ref = np.array([10, 10]) then numpy would
# interpret it as an int array.
points = np.asarray(points, dtype=float)
ref = np.atleast_2d(np.asarray(ref, dtype=float))
if check_all_positive:
if not _all_positive(points):
raise ValueError(
"All values must be larger than 0 in the input points"
)
if not _all_positive(ref):
raise ValueError(
"All values must be larger than 0 in the reference set"
)
nobj = points.shape[1]
if nobj != ref.shape[1]:
raise ValueError(
f"points and ref need to have the same number of columns ({nobj} != {ref.shape[1]})"
)
if nobj == 0:
raise ValueError("The number of columns cannot be 0")
points, npoints, nobj_c = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
ref, ref_size, _ = np2d_to_double_array(
ref, ctype_shape=("size_t", "uint_fast8_t")
)
maximise = _parse_maximise_to_bool_array(maximise, nobj)
return nobj, points, npoints, nobj_c, ref, ref_size, maximise
[docs]
@DocSubstitute()
def igd(
points: ArrayLike,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
) -> float:
"""Inverted Generational Distance (IGD).
.. seealso:: For examples, see :func:`igd_plus`.
For details of the calculation, see :ref:`igd_hausdorf`.
Parameters
----------
points :
${points}
ref :
${ref_set}
maximise :
${maximise}
Returns
-------
A single numerical value.
"""
nobj, points_p, n, d, ref_p, ref_size, maximise_p = _unary_refset_common(
points, ref, maximise
)
if nobj == 1 or nobj > DIMENSION_MAX:
return _igd_python(points, ref, maximise)
return lib.IGD(points_p, n, d, ref_p, ref_size, maximise_p)
[docs]
@DocSubstitute()
def igd_plus(
points: ArrayLike,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
) -> float:
r"""Modified IGD (IGD+).
IGD+ is a Pareto-compliant version of :func:`igd` proposed by :cite:t:`IshMasTanNoj2015igd`.
.. seealso:: For details of the calculation, see :ref:`igd_hausdorf`.
Parameters
----------
points :
${points}
ref :
${ref_set}
maximise :
${maximise}
Returns
-------
A single numerical value.
Examples
--------
>>> dat = np.array([[3.5, 5.5], [3.6, 4.1], [4.1, 3.2], [5.5, 1.5]])
>>> ref = np.array([[1, 6], [2, 5], [3, 4], [4, 3], [5, 2], [6, 1]])
>>> moocore.igd(dat, ref=ref)
1.0627908666722465
>>> moocore.igd_plus(dat, ref=ref)
0.9855036468106652
>>> moocore.avg_hausdorff_dist(dat, ref)
1.0627908666722465
Example 4 by :cite:t:`IshMasTanNoj2015igd` shows a case where IGD gives the
wrong answer. In this case A is better than B in terms of Pareto optimality.
>>> ref = np.array([10, 0, 6, 1, 2, 2, 1, 6, 0, 10]).reshape(-1, 2)
>>> A = np.array([4, 2, 3, 3, 2, 4]).reshape(-1, 2)
>>> B = np.array([8, 2, 4, 4, 2, 8]).reshape(-1, 2)
>>> print(
... f"IGD(A)={moocore.igd(A, ref)} > IGD(B)={moocore.igd(B, ref)}\n"
... + f"and AvgHausdorff(A)={moocore.avg_hausdorff_dist(A, ref)} > "
... + f"AvgHausdorff(B)={moocore.avg_hausdorff_dist(B, ref)},\n"
... + f"which both contradict Pareto optimality. By contrast,\n"
... + f"IGD+(A)={moocore.igd_plus(A, ref)} < IGD+(B)={moocore.igd_plus(B, ref)},"
... + " which is correct."
... )
IGD(A)=3.707092031609239 > IGD(B)=2.59148346584763
and AvgHausdorff(A)=3.707092031609239 > AvgHausdorff(B)=2.59148346584763,
which both contradict Pareto optimality. By contrast,
IGD+(A)=1.482842712474619 < IGD+(B)=2.260112615949154, which is correct.
""" # noqa: D401
nobj, points_p, n, d, ref_p, ref_size, maximise_p = _unary_refset_common(
points, ref, maximise
)
if nobj == 1 or nobj > DIMENSION_MAX:
return _igd_plus_python(points, ref, maximise)
return lib.IGD_plus(points_p, n, d, ref_p, ref_size, maximise_p)
[docs]
@DocSubstitute()
def avg_hausdorff_dist(
points: ArrayLike,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
p: int = 1,
) -> float:
"""Average Hausdorff distance.
.. seealso:: For examples, see :func:`igd_plus`. For details of the calculation, see :ref:`igd_hausdorf`.
Parameters
----------
points :
${points}
ref :
${ref_set}
maximise :
${maximise}
p :
Hausdorff distance parameter. Must be larger than 0.
Returns
-------
A single numerical value.
"""
nobj, points_p, n, d, ref_p, ref_size, maximise_p = _unary_refset_common(
points, ref, maximise
)
if not is_integer_value(p):
raise ValueError("'p' must be an integer")
if p <= 0:
raise ValueError("'p' must be larger than zero")
if nobj == 1 or nobj > DIMENSION_MAX:
return _avg_hausdorff_dist_python(points, ref, maximise, p)
return lib.avg_Hausdorff_dist(
points_p, n, d, ref_p, ref_size, maximise_p, p
)
[docs]
@DocSubstitute()
def epsilon_additive(
points: ArrayLike,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
) -> float:
r"""Additive epsilon metric.
The current implementation uses the naive algorithm that requires
:math:`O(m \cdot |A| \cdot |R|)`, where :math:`m` is the number of
objectives (dimension of vectors), :math:`A` is the input set and :math:`R`
is the reference set.
.. seealso:: For details of the calculation, see :ref:`epsilon_metric`.
Parameters
----------
points :
${points}
ref :
${ref_set}
maximise :
${maximise}
Returns
-------
A single numerical value.
Examples
--------
>>> import numpy as np
>>> dat = np.array([[3.5, 5.5], [3.6, 4.1], [4.1, 3.2], [5.5, 1.5]])
>>> ref = np.array([[1, 6], [2, 5], [3, 4], [4, 3], [5, 2], [6, 1]])
>>> moocore.epsilon_additive(dat, ref=ref)
2.5
>>> moocore.epsilon_mult(dat, ref=ref)
3.5
>>> float(np.exp(moocore.epsilon_additive(np.log(dat), ref=np.log(ref))))
3.5
"""
nobj, points_p, n, d, ref_p, ref_size, maximise_p = _unary_refset_common(
points, ref, maximise
)
if nobj == 1 or nobj > DIMENSION_MAX:
return _epsilon_addi_python(points, ref, maximise)
return lib.epsilon_additive(points_p, n, d, ref_p, ref_size, maximise_p)
[docs]
@DocSubstitute()
def epsilon_mult(
points: ArrayLike,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
) -> float:
"""Multiplicative epsilon metric.
.. warning::
All values in ``points`` and ``ref`` must be larger than 0.
.. seealso:: For examples, see :func:`epsilon_additive`. For details of the calculation, see :ref:`epsilon_metric`.
Parameters
----------
points :
${points}
ref :
${ref_set}
maximise :
${maximise}
Returns
-------
A single numerical value.
"""
nobj, points_p, n, d, ref_p, ref_size, maximise_p = _unary_refset_common(
points, ref, maximise, check_all_positive=True
)
if nobj == 1 or nobj > DIMENSION_MAX:
return _epsilon_mult_python(points, ref, maximise)
return lib.epsilon_mult(points_p, n, d, ref_p, ref_size, maximise_p)
def _hypervolume(points: ArrayLike, ref: ArrayLike) -> float:
_check_dimension_max(points.shape[1], HV_DIMENSION_MAX)
points_p, npoints, nobj = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
ref_buf = ffi.from_buffer("double []", ref)
hv = lib.fpli_hv(points_p, npoints, nobj, ref_buf)
if hv < 0:
raise MemoryError("memory allocation failed")
return hv
[docs]
@DocSubstitute()
def hypervolume(
points: ArrayLike,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
) -> float:
r"""Hypervolume indicator.
Compute the hypervolume metric with respect to a given reference point
assuming minimization of all objectives.
.. seealso:: For details about the hypervolume, see :ref:`hypervolume_metric`.
See :ref:`Benchmarking the exact computation of hypervolume <bench-hv>`.
Parameters
----------
points :
${points}
ref :
${ref_point}
maximise :
${maximise}
Returns
-------
A single numerical value, the hypervolume indicator.
See Also
--------
Hypervolume : Object-oriented interface.
RelativeHypervolume : Compute hypervolume relative to a reference set.
hv_contributions: Hypervolume contribution of individual points.
Notes
-----
For 2D and 3D, the algorithms used
:footcite:p:`FonPaqLop06:hypervolume,BeuFonLopPaqVah09:tec` have :math:`O(n
\log n)` complexity, where :math:`n` is the number of input points. The 3D
case uses the HV3D\ :sup:`+` algorithm :footcite:p:`GueFon2017hv4d`, which
has the same complexity as the HV3D algorithm
:footcite:p:`FonPaqLop06:hypervolume,BeuFonLopPaqVah09:tec`, but it is
faster in practice.
For 4D, the algorithm used is HV4D\ :sup:`+` :footcite:p:`GueFon2017hv4d`,
which has :math:`O(n^2)` complexity. Compared to the `original
implementation <https://github.com/apguerreiro/HVC/>`_, this implementation
correctly handles weakly dominated points and has been further optimized
for speed.
For 5D or higher and up to |HV_INEX_MAX_ROWS| points, the implementation
uses the inclusion-exclusion algorithm :footcite:p:`WuAza2001metrics`,
which has :math:`O(m 2^{n})` time and :math:`O(n\cdot m)` space complexity,
where :math:`m` is the dimension of the points, but it is very fast for
such small sets. For larger number of points, it uses a recursive
algorithm :footcite:p:`FonPaqLop06:hypervolume` that computes 4D
contributions :footcite:p:`GueFon2017hv4d` as the base case, resulting in a
:math:`O(n^{m-2})` time complexity and :math:`O(n)` space complexity in the
worst-case. Experimental results show that the pruning techniques used may
reduce the time complexity even further. The original proposal
:footcite:p:`FonPaqLop06:hypervolume` had the HV3D algorithm as the base
case, giving a time complexity of :math:`O(n^{m-2} \log n)`. Andreia
P. Guerreiro enhanced the numerical stability of the recursive algorithm by
avoiding floating-point comparisons of partial hypervolumes.
The hypervolume of 1D inputs is defined as :code:`max(0, ref - points.min())`.
References
----------
.. footbibliography::
Examples
--------
>>> dat = np.array([[5, 5], [4, 6], [2, 7], [7, 4]])
>>> moocore.hypervolume(dat, ref=[10, 10])
38.0
This function assumes that objectives must be minimized by default. We can
easily specify maximization:
>>> moocore.hypervolume(dat, ref=0, maximise=True)
39.0
Merge all the sets of a dataset by removing the set number column:
>>> dat = moocore.get_dataset("input1.dat")[:, :-1]
>>> len(dat)
100
Dominated points are ignored, so this:
>>> moocore.hypervolume(dat, ref=[10, 10])
93.55331425585321
gives the same hypervolume as this:
>>> dat = moocore.filter_dominated(dat)
>>> len(dat)
6
>>> moocore.hypervolume(dat, ref=[10, 10])
93.55331425585321
"""
# Convert to numpy.array in case the user provides a list. We use
# np.asarray to convert it to floating-point, otherwise if a user inputs
# something like ref = np.array([10, 10]) then numpy would interpret it as
# an int array.
points, points_copied = asarray_maybe_copy(points)
nobj = points.shape[1]
if nobj == 0:
raise ValueError("input points must have at least 1 column")
# Make sure it is a 1D array of length nobj.
ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref")
maximise = _parse_maximise(maximise, nobj)
# FIXME: Do this in C.
if maximise.any():
if not points_copied:
points = points.copy()
points[:, maximise] = -points[:, maximise]
ref = ref.copy()
ref[maximise] = -ref[maximise]
return _hypervolume(points, ref)
[docs]
@DocSubstitute()
class Hypervolume:
"""Object-oriented interface for the hypervolume indicator.
.. seealso:: For examples, see :func:`hypervolume`.
Parameters
----------
ref :
${ref_point}
maximise :
${maximise}
Examples
--------
This function assumes that objectives must be minimized by default. We can
easily specify maximization:
>>> hv_ind = moocore.Hypervolume(ref=0, maximise=True)
>>> hv_ind([[5, 5], [4, 6], [2, 7], [7, 4]])
39.0
>>> hv_ind([[5, 5], [4, 6], [7, 4]])
37.0
See Also
--------
RelativeHypervolume : Compute hypervolume relative to a reference set.
"""
def __init__(
self, ref: ArrayLike, maximise: bool | Sequence[bool | int] = False
) -> None:
self._ref = np.array(ref, dtype=float, ndmin=1)
self._maximise = np.array(maximise, dtype=bool, ndmin=1)
n = len(self._maximise)
if n == 1:
n = len(self._ref)
if n > 1:
self._maximise = np.full((n), self._maximise[0])
elif len(self._ref) == 1:
self._ref = np.full((n), self._ref[0])
elif n != len(self._ref):
raise ValueError(
f"ref and maximise need to have the same length ({len(self._ref)} != {n})"
)
if self._maximise.any():
self._ref[self._maximise] = -self._ref[self._maximise]
self._nobj = n
[docs]
@DocSubstitute()
def __call__(self, points: ArrayLike) -> float:
r"""Compute hypervolume indicator.
Parameters
----------
points :
${points}
Returns
-------
A single numerical value, the hypervolume indicator.
"""
points, points_copied = asarray_maybe_copy(points)
nobj = self._nobj
ref = self._ref
maximise = self._maximise
if nobj == 1:
nobj = points.shape[1]
ref = np.full((nobj), ref[0])
maximise = np.full((nobj), maximise[0])
elif nobj != points.shape[1]:
raise ValueError(
f"points and ref need to have the same number of objectives ({points.shape[1]} != {nobj})"
)
# FIXME: Do this in C.
if maximise.any():
if not points_copied:
points = points.copy()
points[:, maximise] = -points[:, maximise]
return _hypervolume(points, ref)
[docs]
@DocSubstitute()
class RelativeHypervolume(Hypervolume):
r"""Computes the hypervolume value of fronts relative to the hypervolume of a reference front.
The value is computed as :math:`1 - \frac{\text{hyp}(X)}{\text{hyp}(R)}`,
where :math:`X` is an input set (front), :math:`R` is the reference front
and :math:`\text{hyp}()` is calculated by :func:`hypervolume`. Thus, lower
values are better, in contrast to the usual hypervolume. The metric has
also been called *hypervolume relative deviation*
:footcite:p:`BezLopStu2017assessment`.
Parameters
----------
ref :
${ref_point}
ref_set :
Reference set (front) as 2D array. The reference front does not need to
weakly dominate the input sets. If this is required, the user must ensure it.
maximise :
${maximise}
References
----------
.. footbibliography::
Examples
--------
This function assumes that objectives must be minimized by default. We can
easily specify maximization:
>>> ref_set = [[6, 6], [2, 7], [7, 2]]
>>> hv_ind = moocore.RelativeHypervolume(ref=0, ref_set=ref_set, maximise=True)
>>> hv_ind([[5, 5], [4, 6], [2, 7], [7, 4]])
0.0250000
>>> hv_ind([[5, 5], [4, 6], [7, 4]])
0.0749999
>>> hv_ind([[0, 0]])
1.0
>>> hv_ind(ref_set)
0.0
"""
def __init__(
self,
ref: ArrayLike,
ref_set: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
) -> None:
super().__init__(ref=ref, maximise=maximise)
self._ref_set_hv = super().__call__(ref_set)
if self._ref_set_hv == 0:
raise ValueError("hypervolume of 'ref_set' is zero")
[docs]
@DocSubstitute()
def __call__(self, points: ArrayLike) -> float:
r"""Compute relative hypervolume indicator as ``1 - (hypervolume(points, ref=ref) / hypervolume(ref_set, ref=ref))``.
Parameters
----------
points :
${points}
Returns
-------
A single numerical value, the relative hypervolume indicator, which must be minimized.
"""
points_hv = super().__call__(points)
return 1.0 - points_hv / self._ref_set_hv
[docs]
@DocSubstitute()
def hv_contributions(
points: ArrayLike,
/,
ref: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
ignore_dominated: bool = True,
) -> np.ndarray:
r"""Hypervolume contributions of a set of points.
Computes the hypervolume contribution of each point of a set of points with
respect to a given reference point. Duplicated and dominated points have
zero contribution. By default, dominated points are ignored, that is, they
do not affect the contribution of other points. See the Notes below for
more details.
.. seealso:: For details about the hypervolume, see :ref:`hypervolume_metric`.
Parameters
----------
points :
${points}
ref :
${ref_point}
maximise :
${maximise}
ignore_dominated :
Whether dominated points are ignored when computing the contribution of nondominated points.
The value of this parameter has an effect on the return values only if the input contains dominated points.
Setting this to ``False`` slows down the computation significantly.
See the Notes below for a detailed explanation.
Returns
-------
An array of floating-point values as long as the number of rows in ``points``.
Each value is the contribution of the corresponding point in ``points``.
Notes
-----
The hypervolume contribution of point :math:`\vec{p} \in X` is defined as
:math:`\text{hvc}(\vec{p}) = \text{hyp}(X) - \text{hyp}(X \setminus
\{\vec{p}\})`. This definition implies that duplicated points have zero
contribution even if not dominated, because removing one of the duplicates
does not change the hypervolume of the remaining set. Moreover, dominated
points also have zero contribution. However, a point that is dominated by
a single (dominating) nondominated point reduces the contribution of the
latter, because removing the dominating point makes the dominated one
become nondominated.
Handling this special case is non-trivial and makes the computation more
expensive, thus the default (``ignore_dominated=True``) ignores all
dominated points in the input, that is, their contribution is set to zero
and their presence does not affect the contribution of any other point. Setting
``ignore_dominated=False`` will consider dominated points according to the
mathematical definition given above, but the computation will be slower.
When the input only consists of mutually nondominated points, the value of
``ignore_dominated`` does not change the result, but the default value is
significantly faster.
The current implementation uses a :math:`O(n \log n)` dimension-sweep
algorithm for 2D. With ``ignore_dominated=True``, the 3D case uses the
HVC3D algorithm :footcite:p:`GueFon2017hv4d`, which has :math:`O(n \log n)`
complexity. Otherwise, the implementation uses the naive algorithm that
requires calculating the hypervolume :math:`|X|+1` times.
References
----------
.. footbibliography::
Examples
--------
>>> x = np.array([[5, 1], [1, 5], [4, 2], [4, 4], [5, 1]])
>>> moocore.hv_contributions(x, ref=(6, 6), ignore_dominated=True)
... ## Explanation:
... # hvc[(5,1)] = 0 = duplicated
... # hvc[(1,5)] = 3 = (4 - 1) * (6 - 5)
... # hvc[(4,2)] = 3 = (5 - 4) * (5 - 2)
... # hvc[(4,4)] = 0 = dominated
... # hvc[(5,1)] = 0 = duplicated
array([0., 3., 3., 0., 0.])
>>> moocore.hv_contributions(x, ref=(6, 6), ignore_dominated=False)
... ## Explanation:
... # hvc[(5,1)] = 0 = duplicated
... # hvc[(1,5)] = 3 = (4 - 1) * (6 - 5)
... # hvc[(4,2)] = 2 = (5 - 4) * (4 - 2)
... # hvc[(4,4)] = 0 = dominated
... # hvc[(5,1)] = 0 = duplicated
array([0., 3., 2., 0., 0.])
>>> x = np.array([[5, 5], [4, 6], [2, 7], [7, 4]])
>>> moocore.hv_contributions(x, ref=[10, 10])
array([2., 1., 6., 3.])
"""
# Convert to numpy.array in case the user provides a list. We use
# np.asarray to convert it to floating-point, otherwise if a user inputs
# something like ref = np.array([10, 10]) then numpy would interpret it as
# an int array.
points, points_copied = asarray_maybe_copy(points)
nobj = points.shape[1]
if nobj < 2:
raise ValueError("input points must have at least 2 columns")
_check_dimension_max(nobj, HV_DIMENSION_MAX)
ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref")
maximise = _parse_maximise(maximise, nobj)
# FIXME: Do this in C.
if maximise.any():
if not points_copied:
points = points.copy()
points[:, maximise] = -points[:, maximise]
ref = ref.copy()
ref[maximise] = -ref[maximise]
hvc = np.empty(len(points), dtype=float)
hvc_p, _ = np1d_to_double_array(hvc)
points_p, npoints, nobj = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
ref_buf = ffi.from_buffer("double []", ref)
ignore_dominated = bool(ignore_dominated)
lib.hv_contributions(
hvc_p, points_p, npoints, nobj, ref_buf, ignore_dominated
)
return hvc
[docs]
@DocSubstitute()
def hv_approx(
points: ArrayLike,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
nsamples: int = 262_144,
seed: int | np.random.Generator | None = None,
method: Literal["DZ2019-HW", "DZ2019-MC", "Rphi-FWE+"] = "Rphi-FWE+",
) -> float:
r"""Approximate the hypervolume indicator.
Approximate the value of the hypervolume metric with respect to a given
reference point assuming minimization of all objectives. Methods
``"Rphi-FWE+"`` and ``"DZ2019-HW"`` are deterministic and ignore the
parameter ``seed``, while ``method="DZ2019-MC"`` relies on Monte-Carlo
sampling :footcite:p:`DenZha2019approxhv`. All methods tend to get more
accurate with higher values of ``nsamples``, but the increase in accuracy
is not monotonic as shown in the example
:ref:`sphx_glr_auto_examples_plot_hv_approx.py`.
.. seealso:: For details of the calculation, see the Notes section below.
See :ref:`Benchmarks: Approximation of the hypervolume <bench-hvapprox>`.
Parameters
----------
points :
${points}
ref :
${ref_point}
maximise :
${maximise}
nsamples :
Number of samples for Monte-Carlo sampling. Higher values typically produce more
accurate approximations of the true hypervolume, but require more time.
seed :
${random_seed}
method :
Method to approximate the hypervolume.
Returns
-------
A single numerical value, the approximate hypervolume indicator.
See Also
--------
hypervolume, whv_hype
Notes
-----
All available methods approximate the hypervolume as a :math:`(m-1)`-dimensional
integral over the surface of hypersphere :footcite:p:`DenZha2019approxhv`:
.. math::
\widehat{HV}_r(A) = \frac{2\pi^\frac{m}{2}}{\Gamma(\frac{m}{2})}\frac{1}{m 2^m}\frac{1}{n}\sum_{i=1}^n \max_{y \in A} s(w^{(i)}, y)^m
where :math:`m` is the number of objectives, :math:`w^{(i)}` are weights
uniformly distributed on :math:`S_{+}`, i.e., the positive orthant of the
:math:`(m-1)`-D unit hypersphere, :math:`n` is the number of weights
sampled, :math:`\Gamma()` is the gamma function :func:`math.gamma()`, i.e.,
the analytical continuation of the factorial function, and :math:`s(w, y) =
\min_{k=1}^m (r_k - y_k)/w_k`.
In the default ``method="Rphi-FWE+"`` :footcite:p:`Lop2026hvapprox`, the
weights :math:`w^{(i)}, i=1\ldots n` are defined using the deterministic
low-discrepancy sequence :math:`R_\phi` :footcite:p:`Rob2018unreasonable`
mapped to the positive orthant of the hypersphere using a modified version
of Fang and Wang efficient mapping :footcite:p:`FanWan1994numtheory`.
In ``method="DZ2019-HW"`` :footcite:p:`DenZha2019approxhv`, the weights
:math:`w^{(i)}, i=1\ldots n` are defined using a deterministic
low-discrepancy sequence. The weight values depend on their number
(``nsamples``), thus increasing the number of weights may not necessarily
increase accuracy because the set of weights would be different.
In ``method="DZ2019-MC"`` :footcite:p:`DenZha2019approxhv`, the weights
:math:`w^{(i)}, i=1\ldots n` are sampled from the unit normal vector such
that each weight :math:`w = \frac{|x|}{\|x\|_2}` where each component of
:math:`x` is independently sampled from the standard normal distribution
:footcite:p:`Muller1959sphere`.
The original source code in C++/MATLAB for both ``"DZ2019-HW"`` and
``"DZ2019-MC"`` methods can be found `here
<https://github.com/Ksrma/Hypervolume-Approximation-using-polar-coordinate>`_.
:footcite:t:`Lop2026hvapprox` empirically shows that ``"Rphi-FWE+"``
typically produces an approximation error as low as the other methods, with
a computational cost similar to ``"DZ2019-MC"`` and significantly faster
than ``"DZ2019-HW"``.
References
----------
.. footbibliography::
Examples
--------
>>> x = np.array([[5, 5], [4, 6], [2, 7], [7, 4]])
>>> moocore.hypervolume(x, ref=[10, 10])
38.0
>>> moocore.hv_approx(x, ref=[10, 10], method="Rphi-FWE+")
37.99998
>>> moocore.hv_approx(x, ref=[10, 10], method="DZ2019-HW")
37.99996
>>> moocore.hv_approx(x, ref=[10, 10], seed=42, method="DZ2019-MC")
38.00081
Merge all the sets of a dataset by removing the set number column:
>>> x = moocore.get_dataset("input1.dat")[:, :-1]
Dominated points are ignored, so this:
>>> moocore.hv_approx(x, ref=10)
93.5533
>>> moocore.hv_approx(x, ref=10, method="DZ2019-HW")
93.5533
gives the same hypervolume approximation as this:
>>> x = moocore.filter_dominated(x)
>>> moocore.hv_approx(x, ref=10)
93.5533
>>> moocore.hv_approx(x, ref=10, method="DZ2019-HW")
93.5533
The approximation is far from perfect for large number of dimensions:
>>> x = moocore.get_dataset("ran.10pts.9d.10")[:, :-1]
>>> x = moocore.filter_dominated(-x, maximise=True)
>>> x = moocore.normalise(x, to_range=[1, 2])
>>> reference = 0.9
>>> moocore.hypervolume(x, ref=reference, maximise=True)
0.483633123747
>>> moocore.hv_approx(x, ref=reference, maximise=True)
0.48397532301
>>> moocore.hv_approx(x, ref=reference, maximise=True, method="DZ2019-HW")
0.483083547763
"""
# Convert to numpy.array in case the user provides a list. We use
# np.asarray to convert it to floating-point, otherwise if a user inputs
# something like ref = np.array([10, 10]) then numpy would interpret it as
# an int array.
points = np.asarray(points, dtype=float)
nobj = points.shape[1]
if nobj < 2:
if nobj == 1: # Just return the exact hypervolume.
return hypervolume(points, ref, maximise=maximise)
else: # nobj == 0
raise ValueError("input points must have at least 1 column")
_check_dimension_max(nobj, HVAPPROX_DIMENSION_MAX)
ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref")
if not is_integer_value(nsamples) or nsamples <= 0 or nsamples > 2147483648:
raise ValueError(
f"nsamples ({nsamples}) must be a positive integer value smaller than 2147483648"
)
maximise_p = _parse_maximise_to_bool_array(maximise, nobj)
points_p, npoints, nobj = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
ref = ffi.from_buffer("double []", ref)
match method:
case "DZ2019-MC":
seed = _get_seed_for_c(seed)
hv = lib.hv_approx_normal(
points_p, npoints, nobj, ref, maximise_p, nsamples, seed
)
case "DZ2019-HW":
hv = lib.hv_approx_hua_wang(
points_p, npoints, nobj, ref, maximise_p, nsamples
)
case "Rphi-FWE+":
hv = lib.hv_approx_rphi_fang_wang_plus(
points_p, npoints, nobj, ref, maximise_p, nsamples
)
case _:
raise ValueError("Unknown method = {method}")
return hv
[docs]
@DocSubstitute()
def is_nondominated(
points: ArrayLike,
/,
*,
maximise: bool | Sequence[bool | int] = False,
keep_weakly: bool = False,
) -> np.ndarray:
r"""Identify dominated points according to Pareto optimality.
.. seealso:: :ref:`Benchmarks of identifying nondominated points <bench-ndom>`.
Parameters
----------
points :
${points}
maximise :
${maximise}
keep_weakly:
If ``False``, return ``False`` for any duplicates of nondominated points except the first one.
Returns
-------
Returns a boolean array of the same length as the number of rows of points,
where ``True`` means that the point is not dominated by any other point.
See Also
--------
filter_dominated : to filter out dominated points.
any_dominated : to quickly test if any point is (weakly)-dominated.
pareto_rank : to rank points according to Pareto dominance (nondominated sorting).
Notes
-----
Given :math:`n` points of dimension :math:`m`, the current implementation
of :func:`is_nondominated`, :func:`filter_dominated` and
:func:`any_dominated` always uses the best-known :math:`O(n \log n)`
dimension-sweep algorithm :footcite:p:`KunLucPre1975jacm` for :math:`m \leq 3`.
For :math:`m \geq 4`, functions :func:`is_nondominated` and
:func:`filter_dominated` use the best-known :math:`O(n \log^{m-2} n)`
algorithm :footcite:p:`KunLucPre1975jacm` when :math:`n > \KUNGSMALLTHRESHOLD`,
and the naive :math:`O(m n^2)` algorithm otherwise.
Function :func:`any_dominated` always uses the naive algorithm for :math:`m
\geq 4`.
References
----------
.. footbibliography::
Examples
--------
>>> S = np.array([[1, 1], [0, 1], [1, 0], [1, 0]])
>>> moocore.is_nondominated(S)
array([False, True, True, False])
>>> moocore.filter_dominated(S)
array([[0, 1],
[1, 0]])
>>> moocore.is_nondominated(S, keep_weakly=True)
array([False, True, True, True])
>>> moocore.filter_dominated(S, keep_weakly=True)
array([[0, 1],
[1, 0],
[1, 0]])
"""
points = np.asarray(points, dtype=float)
if len(points.shape) != 2:
raise ValueError("'points' must be a matrix")
nrows, nobj = points.shape
if nrows == 0:
return np.array([], dtype=bool)
if nobj < 2:
if nobj == 1: # Handle single-objective inputs
maximise = array_1d_of_length_n(maximise, 1).astype(bool)
if keep_weakly:
best = points.max() if maximise else points.min()
return (points == best).ravel()
else:
nondom = np.zeros(len(points), dtype=bool)
nondom[points.argmax() if maximise else points.argmin()] = True
return nondom
else: # nobj == 0
raise ValueError("input points must have at least 1 column")
_check_dimension_max(nobj, DIMENSION_MAX)
keep_weakly = bool(keep_weakly)
maximise_p = _parse_maximise_to_bool_array(maximise, nobj)
points_p, npoints, nobj = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
nondom = np.empty(nrows, dtype=bool)
# The C library uses uint8_t for boolean vectors.
lib.is_nondominated(
ffi.from_buffer("uint8_t []", nondom.view(np.uint8)),
points_p,
npoints,
nobj,
keep_weakly,
maximise_p,
)
return nondom
[docs]
@DocSubstitute()
def any_dominated(
points: ArrayLike,
/,
*,
maximise: bool | Sequence[bool | int] = False,
keep_weakly: bool = False,
) -> bool:
r"""Test whether the input points contains any (weakly-)dominated point.
This function is equivalent to ``np.logical_not(moocore.is_nondominated(x)).any()``, but significantly faster.
.. seealso:: For details about the calculation, see :func:`is_nondominated`.
Parameters
----------
points :
${points}
maximise :
${maximise}
keep_weakly:
If ``False``, return ``True`` if any point is duplicated.
Returns
-------
Returns ``True`` if the input array contains at least one point that is weakly-dominated by another point in the array. With ``keep_weakly=True``, the point must be dominated not just a duplicate.
Examples
--------
>>> S = np.array([[1, 1], [1, 0], [0, 1], [1, 0]])
>>> moocore.is_nondominated(S)
array([False, True, True, False])
>>> moocore.any_dominated(S)
True
>>> moocore.any_dominated(moocore.filter_dominated(S))
False
>>> S = np.array([[0, 1], [1, 0], [1, 0]])
>>> moocore.is_nondominated(S, keep_weakly=True)
array([ True, True, True])
>>> moocore.any_dominated(S, keep_weakly=True)
False
"""
if keep_weakly:
points = np.unique(points, axis=0)
points = np.asarray(points, dtype=float)
if len(points.shape) != 2:
raise ValueError("'points' must be a matrix")
nrows, nobj = points.shape
if nrows < 2:
if nrows == 0:
raise ValueError("no points in the input points")
else:
return False
if nobj < 2:
if nobj == 1: # Handle single-objective inputs
# If there are more than one row, then something is dominated.
return True
else: # nobj == 0
raise ValueError("input points must have at least 1 column")
_check_dimension_max(nobj, DIMENSION_MAX)
maximise_p = _parse_maximise_to_bool_array(maximise, nobj)
points_p, npoints, nobj = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
res = lib.find_weakly_dominated_point(points_p, npoints, nobj, maximise_p)
return res < nrows
[docs]
@DocSubstitute()
def is_nondominated_within_sets(
points: ArrayLike,
/,
sets: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
keep_weakly: bool = False,
) -> np.ndarray:
r"""Identify dominated points according to Pareto optimality within each set.
Executes the :func:`is_nondominated` function within each set in a dataset
and returns back a 1D array of booleans. This is equivalent to
``apply_within_sets(points, sets, is_nondominated, ...)`` but slightly
faster.
Parameters
----------
points :
${points}
sets :
${sets_of_points}
maximise :
${maximise}
keep_weakly :
If ``False``, return ``False`` for any duplicates of nondominated points.
Returns
-------
Returns a boolean array of the same length as the number of rows of points,
where ``True`` means that the point is not dominated by any other point.
See Also
--------
filter_dominated_within_sets : filter out dominated points.
apply_within_sets : a more general way to apply any function to each set.
Examples
--------
>>> x = np.array([[1, 2, 1], [1, 3, 1], [2, 1, 1], [2, 2, 2]])
>>> moocore.is_nondominated_within_sets(x[:, :-1], x[:, -1])
array([ True, False, True, True])
>>> x = moocore.get_dataset("input1.dat")
>>> nondom_per_set = moocore.is_nondominated_within_sets(x[:, :-1], x[:, -1])
>>> len(nondom_per_set)
100
>>> nondom_per_set # doctest: +ELLIPSIS
array([False, False, True, False, True, False, False, False, False,
True, False, True, True, True, False, True, True, True,
False, True, False, False, False, False, True, False, True,
...
True, True, True, False, True, False, True, True, False,
True, False, False, True, True, False, False, False, False,
False])
>>> x[nondom_per_set, :] # doctest: +ELLIPSIS
array([[ 0.20816431, 4.62275469, 1. ],
[ 0.22997367, 1.11772205, 1. ],
[ 0.58799475, 0.73891181, 1. ],
[ 1.5964888 , 5.98825094, 2. ],
[ 5.2812367 , 3.47800969, 2. ],
[ 2.16315952, 4.7394435 , 2. ],
...
[ 0.6510164 , 9.42381213, 9. ],
[ 1.30291449, 4.50417698, 9. ],
[ 0.62230271, 3.56945324, 10. ],
[ 0.86723965, 1.58599089, 10. ],
[ 6.43135537, 1.00153569, 10. ]])
"""
points = np.asarray(points, dtype=float)
ncols = points.shape[1]
if ncols < 2:
raise ValueError("'points' must have at least 2 columns (2 objectives)")
# FIXME: How can we make this faster?
_, idx, inv = np.unique(sets, return_index=True, return_inverse=True)
# Remember the original position of each element of each set.
idx = [np.flatnonzero(inv == i) for i in idx.argsort()]
points = np.concatenate(
[
is_nondominated(
points.take(g_idx, axis=0),
maximise=maximise,
keep_weakly=keep_weakly,
)
for g_idx in idx
]
)
idx = np.concatenate(idx).argsort()
return points.take(idx, axis=0)
[docs]
@DocSubstitute()
def filter_dominated(
points: ArrayLike,
/,
*,
maximise: bool | Sequence[bool | int] = False,
keep_weakly: bool = False,
) -> np.ndarray:
"""Remove dominated points according to Pareto optimality.
.. seealso:: For examples, see :func:`is_nondominated`.
Parameters
----------
points :
${points}
maximise :
${maximise}
keep_weakly :
If ``False``, delete duplicates of nondominated points.
Returns
-------
Returns the rows of ``points`` where :func:`is_nondominated` is ``True``.
"""
return points[
is_nondominated(points, maximise=maximise, keep_weakly=keep_weakly)
]
[docs]
@DocSubstitute()
def filter_dominated_within_sets(
data: ArrayLike,
/,
*,
maximise: bool | Sequence[bool | int] = False,
keep_weakly: bool = False,
) -> np.ndarray:
"""Given a matrix that represents multiple sets of points (last column gives the set index), filter dominated points within each set.
Executes the :func:`filter_dominated` function within each set in a dataset
and returns back a dataset. This is roughly equivalent to partitioning ``dataset`` according to the last column,
filtering dominated solutions within each partition, and joining back the result.
Parameters
----------
data :
Numpy array of numerical values and set numbers, containing multiple datasets. For example the output of the :func:`read_datasets` function.
maximise :
${maximise}
keep_weakly :
If ``False``, delete duplicates of nondominated points.
Returns
-------
A numpy array where each set only contains nondominated points with respect to the set (last column is the set index).
Points from one set can still dominate points from another set.
Examples
--------
>>> x = moocore.get_dataset("input1.dat")
>>> pf_per_set = moocore.filter_dominated_within_sets(x)
>>> len(pf_per_set)
42
>>> pf_per_set # doctest: +ELLIPSIS
array([[ 0.20816431, 4.62275469, 1. ],
[ 0.22997367, 1.11772205, 1. ],
[ 0.58799475, 0.73891181, 1. ],
[ 1.5964888 , 5.98825094, 2. ],
[ 5.2812367 , 3.47800969, 2. ],
[ 2.16315952, 4.7394435 , 2. ],
...
[ 0.6510164 , 9.42381213, 9. ],
[ 1.30291449, 4.50417698, 9. ],
[ 0.62230271, 3.56945324, 10. ],
[ 0.86723965, 1.58599089, 10. ],
[ 6.43135537, 1.00153569, 10. ]])
>>> pf = moocore.filter_dominated(x[:, :-1])
>>> len(pf)
6
>>> pf
array([[0.20816431, 4.62275469],
[0.22997367, 1.11772205],
[0.58799475, 0.73891181],
[1.54506255, 0.38303122],
[0.17470556, 8.89066343],
[8.57911868, 0.35169752]])
>>> moocore.filter_dominated_within_sets(x[(x[:, 2] >= 3) & (x[:, 2] <= 5), :])
array([[2.60764118, 6.31309852, 3. ],
[3.22509709, 6.1522834 , 3. ],
[0.37731545, 9.02211752, 3. ],
[4.61023932, 2.29231998, 3. ],
[0.2901393 , 8.32259412, 4. ],
[1.54506255, 0.38303122, 4. ],
[4.43498452, 4.13150648, 5. ],
[9.78758589, 1.41238277, 5. ],
[7.85344142, 3.02219054, 5. ],
[0.9017068 , 7.49376946, 5. ],
[0.17470556, 8.89066343, 5. ]])
The above returns sets 3,4,5 with dominated points within each set removed.
See Also
--------
read_datasets : read datasets from a file.
filter_dominated : to be used with a single dataset.
"""
data = np.asarray(data, dtype=float)
if data.shape[1] < 3:
raise ValueError(
"'data' must have at least 3 columns (2 objectives + set column)"
)
is_nondom = is_nondominated_within_sets(
data[:, :-1],
sets=data[:, -1],
maximise=maximise,
keep_weakly=keep_weakly,
)
return data[is_nondom, :]
[docs]
@DocSubstitute()
def pareto_rank(
points: ArrayLike, /, *, maximise: bool | Sequence[bool | int] = False
) -> np.ndarray:
r"""Rank points according to Pareto-optimality (nondominated sorting).
The function :func:`pareto_rank` is meant to be used like
:func:`numpy.argsort`, but it assigns indexes according to Pareto
dominance, where rank 0 indicates those solutions not dominated by any
other solution in the input set. Duplicated points are assigned the same
rank. The resulting ranking can be used to partition points into different
lists or arrays, each of them being mutually nondominated
:footcite:p:`Deb02nsga2` (see examples below).
.. seealso:: :ref:`Benchmarks of nondominated Sorting (Pareto ranking) <bench-ndsort>`.
Parameters
----------
points :
${points}
maximise :
${maximise}
Returns
-------
An integer vector of the same length as the number of rows of ``points`` with values within ``[0, len(points) - 1]``, where each value gives the rank of each point (lower is better).
Notes
-----
Given :math:`x,y \in\mathbb{R}^m`, let :math:`x \prec y` denote that
:math:`x` dominates :math:`y` according to Pareto optimality.
Nondominated sorting, also called the *layers-of-maxima problem* :footcite:p:`BucGoo2004maxima`,
partitions a finite set of points :math:`X \subset \mathbb{R}^m` into an
ordered sequence of `fronts`, :math:`F_0, F_1, \dots, F_k`, where
:math:`0 \leq k < |X|`, satisfying the following conditions:
#. All points are allocated to a front: :math:`\bigcup_{i=0}^{k} F_i = X`
#. Each point is allocated to only one front: :math:`F_i \cap F_j = \emptyset,\; \forall i,j \in\{0,1,\dots,k\},\,i\neq j`
#. Points within a front are mutually nondominated: :math:`\nexists x,y \in F_i, \,x \prec y,\; \forall i\in\{0,1,\dots,k\}`
#. The first front contains points not dominated by any other point: :math:`\forall y \in F_0,\,\nexists x \in X,\, x\prec y`
#. Every point in front :math:`i` is dominated by at least one point in front :math:`i-1`: :math:`\forall y \in F_i,\, \exists x \in F_{i-1},\, x\prec y,\; \forall i\in\{1,2,\dots,k\}`
The rank returned by :func:`pareto_rank` is the index :math:`i \in
\{0,1,\dots,k\}` of the front :math:`F_i` allocated to each point. If all
points are mutually nondominated, there is only one front :math:`(k=0)`. If
each front contains one point, then :math:`k=|X| - 1`.
With :math:`m=2`, the code uses the best-known :math:`O(n \log n)`
algorithm by :footcite:t:`Jen03`, where :math:`n=|X|`. When :math:`m \geq
3`, it uses the naive algorithm that identifies one `front` at a time,
which requires :math:`O(n^2\log n)` for :math:`m=3`, and :math:`O(n^2
\log^{m-2} n)` for :math:`m \geq 4`.
References
----------
.. footbibliography::
Examples
--------
>>> x = moocore.get_dataset("input1.dat")[:, :2]
>>> ranks = moocore.pareto_rank(x)
>>> ranks
array([ 4, 8, 0, 11, 0, 3, 7, 1, 3, 0, 8, 4, 5, 4, 11, 4, 4,
5, 7, 3, 8, 12, 8, 9, 5, 10, 6, 2, 7, 3, 10, 7, 2, 5,
2, 7, 1, 2, 9, 0, 11, 6, 7, 10, 13, 3, 6, 3, 0, 9, 9,
0, 2, 13, 1, 6, 7, 6, 6, 10, 4, 13, 6, 8, 12, 13, 4, 8,
5, 1, 12, 10, 3, 8, 9, 6, 7, 6, 6, 9, 5, 2, 3, 4, 7,
3, 8, 3, 2, 10, 1, 2, 12, 1, 2, 9, 4, 2, 5, 2],
dtype=int32)
We can now split the original set into a list of nondominated sets ordered by Pareto rank:
>>> fronts = [x.compress((g == ranks), axis=0) for g in np.unique(ranks)]
>>> len(fronts)
14
The first element is the set of points not dominated by anything else:
>>> np.array_equal(fronts[0], moocore.filter_dominated(x))
True
"""
points, points_copied = asarray_maybe_copy(points)
nrows, nobj = points.shape
if nobj < 2:
if nobj == 1:
points = points.ravel()
if maximise:
points = -points
# FIXME: Can we do the same faster?
return np.unique(points, return_inverse=True)[1]
else: # nobj == 0
return np.zeros(shape=nrows, dtype=int)
_check_dimension_max(nobj, DIMENSION_MAX)
maximise = _parse_maximise(maximise, nobj)
if maximise.any():
# FIXME: Do this in C.
if not points_copied:
points = points.copy()
points[:, maximise] = -points[:, maximise]
points_p, npoints, nobj = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
ranks = np.empty(nrows, dtype=np.intc())
lib.pareto_rank(ffi.from_buffer("int []", ranks), points_p, npoints, nobj)
return ranks
[docs]
@DocSubstitute()
def normalise(
points: ArrayLike,
/,
to_range: ArrayLike = (0.0, 1.0),
*,
lower: ArrayLike = np.nan,
upper: ArrayLike = np.nan,
maximise: bool | Sequence[bool | int] = False,
) -> np.ndarray:
"""Normalise points per coordinate to a range, e.g., ``to_range = [1,2]``, where the minimum value will correspond to 1 and the maximum to 2.
Parameters
----------
points :
${points}
to_range :
Range composed of two numerical values. Normalise values to this range. If the objective is maximised, it is normalised to ``(to_range[1], to_range[0])`` instead.
upper, lower:
Bounds on the values. If :data:`numpy.nan`, the maximum and minimum values of each coordinate are used.
maximise :
${maximise}
Returns
-------
Returns the points normalised as requested.
Examples
--------
>>> dat = np.array([[3.5, 5.5], [3.6, 4.1], [4.1, 3.2], [5.5, 1.5]])
>>> moocore.normalise(dat)
array([[0. , 1. ],
[0.05 , 0.65 ],
[0.3 , 0.425],
[1. , 0. ]])
>>> moocore.normalise(dat, to_range=[1, 2], lower=[3.5, 3.5], upper=5.5)
array([[1. , 2. ],
[1.05, 1.3 ],
[1.3 , 0.85],
[2. , 0. ]])
"""
# normalise() modifies the points, so we need to create a copy.
# order='C' is needed for np2d_to_double_array()
points = np.array(points, dtype=float, order="C")
npoints, nobj = points.shape
if nobj < 2:
raise ValueError("'points' must have at least two columns")
_check_dimension_max(nobj, DIMENSION_MAX)
to_range = np.asarray(to_range, dtype=float)
if to_range.shape[0] != 2:
raise ValueError("'to_range' must have length 2")
lower = array_1d_of_length_n(
np.asarray(lower, dtype=float), nobj, name="lower"
)
upper = array_1d_of_length_n(
np.asarray(upper, dtype=float), nobj, name="upper"
)
if np.any(np.isnan(lower)):
lower = np.where(np.isnan(lower), points.min(axis=0), lower)
if np.any(np.isnan(upper)):
upper = np.where(np.isnan(upper), points.max(axis=0), upper)
maximise_p = _parse_maximise_to_bool_array(maximise, nobj)
points_p, npoints, nobj = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
lbound_p = ffi.from_buffer("double []", lower)
ubound_p = ffi.from_buffer("double []", upper)
lib.agree_normalise(
points_p,
npoints,
nobj,
maximise_p,
to_range[0],
to_range[1],
lbound_p,
ubound_p,
)
# We can return points directly because we only changed the values, not the shape.
return points
[docs]
@DocSubstitute()
def eaf(
points: ArrayLike, /, sets: ArrayLike, *, percentiles: Sequence[float] = ()
) -> np.ndarray:
r"""Exact computation of the Empirical Attainment Function (EAF).
.. seealso:: For the definition of EAF, see :ref:`eaf_computation`.
Parameters
----------
points :
${points}
sets :
${sets_of_points}
percentiles :
List indicating which percentiles are computed. By default, all possible percentiles are calculated.
Returns
-------
EAF data points, with the same number of columns as the input argument, but a different number of rows. The last column represents the EAF percentile for that data point.
See Also
--------
mooplot.plot_eaf: Plotting the EAF.
Notes
-----
In the current implementation, the EAF is computed using the algorithms
proposed by :footcite:t:`FonGueLopPaq2011emo`, which have complexity
:math:`\Theta(m\log m + nm)` in 2D and :math:`O(n^2 m \log m)` in 3D, where
:math:`n` is the number of input sets and :math:`m` is the total number of
input points.
References
----------
.. footbibliography::
Examples
--------
>>> x = moocore.get_dataset("input1.dat")
>>> moocore.eaf(x[:, :-1], x[:, -1]) # doctest: +ELLIPSIS
array([[ 0.17470556, 8.89066343, 10. ],
[ 0.20816431, 4.62275469, 10. ],
[ 0.22997367, 1.11772205, 10. ],
[ 0.58799475, 0.73891181, 10. ],
[ 1.54506255, 0.38303122, 10. ],
[ 8.57911868, 0.35169752, 10. ],
[ 0.20816431, 8.89066343, 20. ],
[ 0.2901393 , 8.32259412, 20. ],
...
[ 9.78758589, 2.8124162 , 90. ],
[ 1.13096306, 9.72645436, 100. ],
[ 2.71891214, 8.84691923, 100. ],
[ 3.34035397, 7.49376946, 100. ],
[ 4.43498452, 6.94327481, 100. ],
[ 4.96525837, 6.20957074, 100. ],
[ 7.92511295, 3.92669598, 100. ]])
>>> moocore.eaf(x[:, :-1], x[:, -1], percentiles=[0, 50, 100]) # doctest: +ELLIPSIS
array([[ 0.17470556, 8.89066343, 0. ],
[ 0.20816431, 4.62275469, 0. ],
[ 0.22997367, 1.11772205, 0. ],
[ 0.58799475, 0.73891181, 0. ],
[ 1.54506255, 0.38303122, 0. ],
[ 8.57911868, 0.35169752, 0. ],
[ 0.53173087, 9.73244829, 50. ],
[ 0.62230271, 9.02211752, 50. ],
[ 0.79293574, 8.89066343, 50. ],
[ 0.9017068 , 8.32259412, 50. ],
[ 0.97468676, 7.65893644, 50. ],
[ 1.06855707, 7.49376946, 50. ],
[ 1.54506255, 6.7102429 , 50. ],
[ 1.5964888 , 5.98825094, 50. ],
[ 2.16315952, 4.7394435 , 50. ],
[ 2.85891341, 4.49240941, 50. ],
[ 3.34035397, 2.89377444, 50. ],
[ 4.61023932, 2.87955367, 50. ],
[ 4.96525837, 2.29231998, 50. ],
[ 7.04694467, 1.83484358, 50. ],
[ 9.7398055 , 1.00153569, 50. ],
[ 1.13096306, 9.72645436, 100. ],
[ 2.71891214, 8.84691923, 100. ],
[ 3.34035397, 7.49376946, 100. ],
[ 4.43498452, 6.94327481, 100. ],
[ 4.96525837, 6.20957074, 100. ],
[ 7.92511295, 3.92669598, 100. ]])
"""
points = np.asarray(points, dtype=float)
ncols = points.shape[1]
if ncols < 2:
raise ValueError("'points' must have at least 2 columns")
if ncols > 3:
raise NotImplementedError(
"Only 2D or 3D datasets are currently supported for computing the EAF"
)
if len(sets) != points.shape[0]:
raise ValueError(
"'sets' must have the same length as the number of rows of 'points'"
)
_, cumsizes = np.unique(sets, return_counts=True)
nsets = len(cumsizes)
cumsizes = np.cumsum(cumsizes)
cumsizes_p, ncumsizes = np1d_to_int_array(cumsizes)
percentiles = np.ravel(percentiles)
if len(percentiles) == 0:
percentiles = np.arange(1.0, nsets + 1) * (100.0 / nsets)
else:
percentiles = np.unique(np.asarray(percentiles, dtype=float))
percentile_p, npercentiles = np1d_to_double_array(percentiles)
# Get C pointers + matrix size for calling CFFI generated extension module
points_p, _, nobj = np2d_to_double_array(points)
eaf_npoints = ffi.new("int *")
eaf_data_p = lib.eaf_compute_matrix(
eaf_npoints,
points_p,
nobj,
cumsizes_p,
ncumsizes,
percentile_p,
npercentiles,
)
eaf_npoints = eaf_npoints[0]
eaf_buf = ffi.buffer(
eaf_data_p, ffi.sizeof("double") * eaf_npoints * (ncols + 1)
)
return np.frombuffer(eaf_buf).reshape((eaf_npoints, -1))
class VorobtResult(NamedTuple):
threshold: float
ve: np.ndarray
avg_hyp: float
[docs]
@DocSubstitute()
def vorob_t(
points: ArrayLike, /, sets: ArrayLike, *, ref: ArrayLike
) -> VorobtResult:
r"""Compute Vorob'ev threshold and expectation.
.. seealso:: For the definition of these concepts, see :ref:`vorobev_expectation_and_deviation`.
Parameters
----------
points :
${points}
sets :
${sets_of_points}
ref :
${ref_point}
Returns
-------
The return value is a :class:`~typing.NamedTuple` with the following attributes:
threshold : float
The threshold is returned as a percentile :math:`\beta^{*} \in [0,100]`.
avg_hyp : float
The hypervolume of the expectation (average hypervolume).
ve : np.ndarray
The expectation :math:`\mathcal{Q}_{\beta^{*}}` as a 2D array.
See Also
--------
vorob_dev : Compute Vorob'ev deviation.
Examples
--------
>>> CPFs = moocore.get_dataset("CPFs.txt.xz")
>>> res = moocore.vorob_t(CPFs[:, :-1], sets=CPFs[:, -1], ref=(2, 200))
>>> res.threshold
44.140625
>>> res.avg_hyp
8943.333191728081
>>> res.ve.shape
(213, 2)
"""
points = np.asarray(points, dtype=float)
nobj = points.shape[1]
if nobj < 2:
raise ValueError("'points' must have at least 2 columns")
hv_ind = Hypervolume(ref=ref)
avg_hyp = np.mean(apply_within_sets(points, sets, hv_ind))
prev_hyp = diff = np.inf # hypervolume of quantile at previous step
a = 0.0
b = 100.0
while diff != 0:
c = (a + b) / 2.0
eaf_res = eaf(points, sets=sets, percentiles=c)[:, :nobj]
tmp = hv_ind(eaf_res)
if tmp > avg_hyp:
a = c
else:
b = c
diff = prev_hyp - tmp
prev_hyp = tmp
return VorobtResult(threshold=c, ve=eaf_res, avg_hyp=float(avg_hyp))
[docs]
@DocSubstitute()
def vorob_dev(
points: ArrayLike,
/,
sets: ArrayLike,
*,
ref: ArrayLike,
ve: ArrayLike | None = None,
) -> float:
r"""Compute Vorob'ev deviation.
Parameters
----------
points :
${points}
sets :
${sets_of_points}
ref :
${ref_point}
ve :
Vorob'ev expectation, e.g., as returned by :func:`vorob_t`. If not
provided, it is calculated by calling ``vorob_t(data, sets, ref)``.
Returns
-------
Vorob'ev deviation.
See Also
--------
vorob_t : Compute Vorob'ev threshold and expectation.
Notes
-----
For more background, see :footcite:t:`BinGinRou2015gaupar,Molchanov2005theory,CheGinBecMol2013moda`.
References
----------
.. footbibliography::
Examples
--------
>>> CPFs = moocore.get_dataset("CPFs.txt.xz")
>>> vd = moocore.vorob_dev(CPFs[:, :-1], sets=CPFs[:, -1], ref=(2, 200))
>>> vd
3017.12989402326
"""
points = np.asarray(points, dtype=float)
nobj = points.shape[1]
if nobj < 2:
raise ValueError("'points' must have at least 2 columns")
# Hypervolume of the symmetric difference between A and B:
# 2 * H(AUB) - H(A) - H(B)
hv_ind = Hypervolume(ref=ref)
if ve is None:
res = vorob_t(points, sets=sets, ref=ref)
ve = res.ve
h1 = res.avg_hyp
else:
h1 = np.mean(apply_within_sets(points, sets, hv_ind))
h2 = hv_ind(ve)
vd = 2 * np.mean(
apply_within_sets(points, sets, lambda g: hv_ind(np.vstack((g, ve))))
)
return float(vd - h1 - h2)
[docs]
@DocSubstitute()
def eafdiff(
x: ArrayLike,
y: ArrayLike,
/,
*,
intervals: int | None = None,
maximise: bool | Sequence[bool | int] = False,
rectangles: bool = False,
) -> np.ndarray:
"""Compute empirical attainment function (EAF) differences.
Calculate the differences between the empirical attainment functions of two
data sets.
.. warning::
The current implementation only supports 2 objectives.
Parameters
----------
x, y : ArrayLike
Numpy matrices corresponding to the input data of left and right sides,
respectively. Each data frame has at least three columns, the third one
being the set of each point. See also :func:`read_datasets`.
intervals :
The absolute range of the differences :math:`[0, 1]` is partitioned into the number of intervals provided.
maximise :
${maximise}
rectangles :
If ``True``, the output is in the form of rectangles of the same color.
Returns
-------
With ``rectangles=False``, a matrix with three columns, The first two
columns describe the points where there is a transition in the value of
the EAF differences.
With ``rectangles=True``, a matrix with five columns, where the first 4
columns give the coordinates of two corners of each rectangle. In both
cases, the last column gives the difference in terms of sets in ``x``
minus sets in ``y`` that attain each point (i.e., negative values are
differences in favour ``y``).
See Also
--------
mooplot.eafdiffplot: Plotting EAF differences.
Examples
--------
>>> from io import StringIO
>>> A1 = moocore.read_datasets(
... StringIO('''
... 3 2
... 2 3
...
... 2.5 1
... 1 2
...
... 1 2''')
... )
>>> A2 = moocore.read_datasets(
... StringIO('''
... 4 2.5
... 3 3
... 2.5 3.5
...
... 3 3
... 2.5 3.5
...
... 2 1''')
... )
>>> moocore.eafdiff(A1, A2)
array([[ 1. , 2. , 2. ],
[ 2. , 1. , -1. ],
[ 2.5, 1. , 0. ],
[ 2. , 2. , 1. ],
[ 2. , 3. , 2. ],
[ 3. , 2. , 2. ],
[ 2.5, 3.5, 0. ],
[ 3. , 3. , 0. ],
[ 4. , 2.5, 1. ]])
>>> moocore.eafdiff(A1, A2, rectangles=True)
array([[ 2. , 1. , 2.5, 2. , -1. ],
[ 1. , 2. , 2. , inf, 2. ],
[ 2.5, 1. , inf, 2. , 0. ],
[ 2. , 2. , 3. , 3. , 1. ],
[ 2. , 3.5, 2.5, inf, 2. ],
[ 2. , 3. , 3. , 3.5, 2. ],
[ 3. , 2.5, 4. , 3. , 2. ],
[ 3. , 2. , inf, 2.5, 2. ],
[ 4. , 2.5, inf, 3. , 1. ]])
"""
x = np.asarray(x, dtype=float)
y = np.asarray(y, dtype=float)
assert x.shape[1] == y.shape[1], (
"'x' and 'y' must have the same number of columns"
)
nobj = x.shape[1] - 1
assert nobj == 2
maximise = _parse_maximise(maximise, nobj=nobj)
# The C code expects points within a set to be contiguous.
x = x[x[:, -1].argsort(), :]
y = y[y[:, -1].argsort(), :]
_, cumsizes_x = np.unique(x[:, -1], return_counts=True)
_, cumsizes_y = np.unique(y[:, -1], return_counts=True)
cumsizes_x = np.cumsum(cumsizes_x)
cumsizes_y = np.cumsum(cumsizes_y)
cumsizes = np.concatenate((cumsizes_x, cumsizes_x[-1] + cumsizes_y))
nsets = len(cumsizes)
data = np.vstack((x[:, :-1], y[:, :-1]))
if maximise.any():
data[:, maximise] = -data[:, maximise]
if intervals is None:
intervals = int(nsets / 2.0)
else:
assert is_integer_value(intervals)
intervals = min(intervals, int(nsets / 2.0))
data_p, _, nobj_int = np2d_to_double_array(data)
cumsizes_p, nsets = np1d_to_int_array(cumsizes)
eaf_npoints = ffi.new("int *")
if rectangles:
eaf_data = lib.eafdiff_compute_rectangles(
eaf_npoints, data_p, nobj_int, cumsizes_p, nsets, intervals
)
ncols = 2 * nobj + 1 # 2x2D points + color
else:
eaf_data = lib.eafdiff_compute_matrix(
eaf_npoints, data_p, nobj_int, cumsizes_p, nsets, intervals
)
ncols = nobj + 1 # 2D points + color
eaf_npoints = eaf_npoints[0]
eaf_data = ffi.buffer(eaf_data, ffi.sizeof("double") * eaf_npoints * ncols)
eaf_data = np.frombuffer(eaf_data).reshape((eaf_npoints, -1))
# FIXME: We should remove duplicated rows in C code.
eaf_data = unique_nosort(eaf_data, axis=0)
# FIXME: Check that we do not generate duplicated nor overlapping
# rectangles with different colors. That would be a bug.
if maximise.any():
if rectangles:
maximise = np.concatenate((maximise, maximise))
maximise = np.flatnonzero(
maximise
) # Using bool directly misses the color column.
eaf_data[:, maximise] = -eaf_data[:, maximise]
return eaf_data
# def get_diff_eaf(x, y, intervals=None, debug=False):
# if np.min(x[:, -1]) != 1 or np.min(y[:, -1]) != 1:
# raise ValueError("x and y should contain set numbers starting from 1")
# ycopy = np.copy(
# y
# ) # Do hard copy so that the matrix is not corrupted. This could be optimised
# ycopy[:, -1] = ycopy[:, -1] + np.max(
# x[:, -1]
# ) # Make Y sets start from end of X sets
# data = np.vstack((x, ycopy)) # Combine X and Y datasets to one matrix
# nsets = len(np.unique(data[:, -1]))
# if intervals is None:
# intervals = nsets / 2.0
# else:
# intervals = min(intervals, nsets / 2.0)
# intervals = int(intervals)
# data = np.ascontiguousarray(
# np.asarray(data, dtype=float)
# ) # C function requires contiguous data
# num_data_columns = data.shape[1]
# data_p, npoints, ncols = np2d_to_double_array(data)
# eaf_npoints = ffi.new("int *", 0)
# sizeof_eaf = ffi.new("int *", 0)
# nsets = ffi.cast("int", nsets) # Get num of sets from data
# intervals = ffi.cast("int", intervals)
# debug = ffi.cast("bool", debug)
# eaf_diff_data = lib.compute_eafdiff_(
# data_p, ncols, npoints, nsets, intervals, eaf_npoints, sizeof_eaf, debug
# )
# eaf_buf = ffi.buffer(eaf_diff_data, sizeof_eaf[0])
# eaf_arr = np.frombuffer(eaf_buf)
# # The C code gets diff EAF in Column Major order so I return it in column major order than transpose to fix into row major order
# return np.reshape(eaf_arr, (num_data_columns, -1)).T
[docs]
@DocSubstitute()
def whv_rect(
points: ArrayLike,
/,
rectangles: ArrayLike,
*,
ref: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
) -> float:
"""Compute weighted hypervolume given a set of rectangles.
.. seealso:: For examples, see :func:`total_whv_rect`.
.. warning::
The current implementation only supports 2 objectives.
Parameters
----------
points :
${points}
rectangles :
${whv_rectangles}
ref :
${ref_point}
maximise :
${maximise}
Returns
-------
A single numerical value, the weighted hypervolume.
See Also
--------
total_whv_rect, whv_hype
"""
points, _ = asarray_maybe_copy(points)
nobj = points.shape[1]
if nobj != 2:
raise NotImplementedError("Only 2D points are currently supported")
if rectangles.shape[1] != 5:
raise ValueError(
"Invalid number of columns in 'rectangles' (should be 5)"
)
maximise = _parse_maximise(maximise, nobj)
if maximise.any():
raise NotImplementedError("Only minimization is currently supported")
# FIXME: Move to C code.
# maximise = _parse_maximise(maximise, nobj)
# if maximise.any():
# x[:, maximise] = -x[:, maximise]
# ref[maximise] = -ref[maximise]
# if maximise.all():
# rectangles[:, :4] = -rectangles[:, :4]
# else:
# pos = np.flatnonzero(maximise) + [0,2]
# rectangles[:, pos] = -rectangles[:, pos]
ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref")
ref = ffi.from_buffer("double []", ref)
points, npoints, _ = np2d_to_double_array(points)
rectangles, rectangles_nrow, _ = np2d_to_double_array(rectangles)
hv = lib.rect_weighted_hv2d(
points, npoints, rectangles, rectangles_nrow, ref
)
return hv
def get_ideal(x: np.ndarray, maximise: np.ndarray) -> np.ndarray:
# FIXME: Is there a better way to do this?
lower = x.min(axis=0)
upper = x.max(axis=0)
return np.where(maximise, upper, lower)
[docs]
@DocSubstitute()
def total_whv_rect(
points: ArrayLike,
/,
rectangles: ArrayLike,
*,
ref: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
ideal: ArrayLike | None = None,
scalefactor: float = 0.1,
) -> float:
r"""Compute total weighted hypervolume given a set of rectangles.
Calculates the hypervolume weighted by a set of rectangles (with zero
weight outside the rectangles). The function :func:`total_whv_rect()`
calculates the total weighted hypervolume as :func:`hypervolume()` ``+
scalefactor * abs(prod(ref - ideal)) *`` :func:`whv_rect()`. The
details of the computation are given by :footcite:t:`DiaLop2020ejor`.
.. warning::
The current implementation only supports 2 objectives.
Parameters
----------
points :
${points}
rectangles :
${whv_rectangles}
ref :
${ref_point}
maximise :
${maximise}
ideal :
Ideal point as a vector of numerical values. If ``None``, it is calculated as minimum (or maximum if maximising that objective) of each objective in the input points.
scalefactor :
Real value within :math:`(0,1]` that scales the overall weight of the differences. This is parameter psi (:math:`\psi`) in :footcite:t:`DiaLop2020ejor`.
Returns
-------
A single numerical value, the weighted hypervolume.
See Also
--------
hypervolume, whv_rect, whv_hype
References
----------
.. footbibliography::
Examples
--------
>>> rectangles = np.array(
... [
... [1.0, 3.0, 2.0, np.inf, 1],
... [2.0, 3.5, 2.5, np.inf, 2],
... [2.0, 3.0, 3.0, 3.5, 3],
... ]
... )
>>> whv_rect([[2, 2]], rectangles, ref=6)
4.0
>>> whv_rect([[2, 1]], rectangles, ref=6)
4.0
>>> whv_rect([[1, 2]], rectangles, ref=6)
7.0
>>> total_whv_rect([[2, 2]], rectangles, ref=6, ideal=1)
26.0
>>> total_whv_rect([[2, 1]], rectangles, ref=6, ideal=1)
30.0
>>> total_whv_rect([[1, 2]], rectangles, ref=6, ideal=1)
37.5
"""
points = np.asarray(points, dtype=float)
nobj = points.shape[1]
if nobj != 2:
raise NotImplementedError("Only 2D datasets are currently supported")
if rectangles.shape[1] != 5:
raise ValueError(
"Invalid number of columns in 'rectangles' (should be 5)"
)
if scalefactor <= 0 or scalefactor > 1:
raise ValueError("'scalefactor' must be within (0,1]")
ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref")
maximise = _parse_maximise(maximise, nobj)
whv = whv_rect(points, rectangles, ref=ref, maximise=maximise)
hv = hypervolume(points, ref=ref) # FIXME: maximise = maximise)
if ideal is None:
# FIXME: Should we include the range of the rectangles here?
ideal = get_ideal(points, maximise=maximise)
else:
ideal = array_1d_of_length_n(
np.asarray(ideal, dtype=float), nobj, name="ideal"
)
beta = scalefactor * abs((ref - ideal).prod())
return float(hv + beta * whv)
[docs]
@DocSubstitute()
def largest_eafdiff(
x: list,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
intervals: int = 5,
ideal: ArrayLike | None = None,
) -> tuple[tuple[int, int], float]:
"""Identify largest EAF differences.
Given a list of datasets, return the indexes of the pair with the largest
EAF differences according to the method proposed by :footcite:t:`DiaLop2020ejor`.
.. warning::
The current implementation only supports 2 objectives.
Parameters
----------
x :
A list of matrices with at least 3 columns (last column indicates the set).
ref :
Reference point as a 1D vector. Must be same length as a single point in the input data.
maximise :
${maximise}
intervals :
The absolute range of the differences :math:`[0, 1]` is partitioned into the number of intervals provided.
ideal :
Ideal point as a vector of numerical values. If ``None``, it is calculated as minimum (or maximum if maximising that objective) of each objective in the input data.
Returns
-------
pair : tuple[int,int]
Pair of indexes into the list that give the largest EAF difference.
value : float
Value of the largest difference.
See Also
--------
eafdiff, whv_rect
References
----------
.. footbibliography::
Examples
--------
>>> import pandas as pd
>>> df = pd.read_csv(moocore.get_dataset_path("tpls50x20_1_MWT.csv"))
>>> df
algorithm Makespan WeightedTardiness run
0 1to2 4280.0 10231.0 1.0
1 1to2 4238.0 10999.0 1.0
2 1to2 4137.0 11737.0 1.0
3 1to2 4024.0 14871.0 1.0
4 1to2 4014.0 17825.0 1.0
... ... ... ... ...
1506 double 4048.0 14755.0 15.0
1507 double 3923.0 25507.0 15.0
1508 double 3890.0 29567.0 15.0
1509 double 3862.0 31148.0 15.0
1510 double 4413.0 9894.0 15.0
<BLANKLINE>
[1511 rows x 4 columns]
>>> nadir = df.iloc[:, 1:3].max().to_numpy()
>>> # Split by column 'algorithm'
>>> x = [
... g.drop("algorithm", axis=1) for i, g in df.groupby("algorithm", sort=False)
... ]
>>> moocore.largest_eafdiff(x, ref=nadir)
((2, 5), 777017.0)
"""
n = len(x)
if n == 0:
raise ValueError("Empty list")
x = [np.asarray(z, dtype=float) for z in x]
nobj = x[0].shape[1] - 1
if nobj != 2:
raise NotImplementedError("Only 2D datasets are currently supported")
maximise = _parse_maximise(maximise, nobj)
if ideal is None:
ideal = get_ideal(
np.concatenate(
[[z[:, :-1].min(axis=0), z[:, :-1].max(axis=0)] for z in x],
axis=0,
),
maximise=maximise,
)
ideal = ideal.reshape(1, -1)
best_value = 0
for a in range(n - 1):
for b in range(a + 1, n):
diff = eafdiff(
x[a],
x[b],
intervals=intervals,
maximise=maximise,
rectangles=True,
)
# Set color to 1
a_rectangles = diff[diff[:, -1] >= 1, :]
a_rectangles[:, -1] = 1
a_value = whv_rect(
ideal, rectangles=a_rectangles, ref=ref, maximise=maximise
)
b_rectangles = diff[diff[:, -1] <= -1, :]
b_rectangles[:, -1] = 1
b_value = whv_rect(
ideal, rectangles=b_rectangles, ref=ref, maximise=maximise
)
value = min(a_value, b_value)
if value > best_value:
best_value = value
best_pair = (a, b)
return best_pair, best_value
[docs]
@DocSubstitute()
def whv_hype(
points: ArrayLike,
/,
*,
ref: ArrayLike,
ideal: ArrayLike,
maximise: bool | Sequence[bool | int] = False,
nsamples: int = 100000,
seed: int | np.random.Generator | None = None,
dist: Literal["uniform", "point", "exponential"] = "uniform",
mu: numbers.Real | ArrayLike | None = None,
) -> float:
r"""Approximation of the (weighted) hypervolume by Monte-Carlo sampling (2D only).
Return an estimation of the hypervolume of the space dominated by the input
points following the procedure described by :footcite:t:`AugBadBroZit2009gecco`.
A weight distribution describing user preferences may be specified.
.. warning::
The current implementation only supports 2 objectives.
Parameters
----------
points :
${points}
ref :
${ref_point}
ideal :
Ideal point as a numpy array or list. Must have same length as the number of columns of the dataset.
maximise :
${maximise}
nsamples :
Number of samples for Monte-Carlo sampling. Higher values give more
accurate approximation of the true hypervolume but require more time.
seed :
${random_seed}
dist :
Weight distribution :footcite:p:`AugBadBroZit2009gecco`. The ones currently supported are:
``'uniform'``
corresponds to the default hypervolume (unweighted). ``mu`` should be ``None``.
``'point'``
describes a goal in the objective space, where ``mu`` gives the
coordinates of the goal (1D Numpy array). The resulting weight distribution is a
multivariate normal distribution centred at the goal.
``'exponential'``
describes an exponential distribution with rate parameter ``1/mu``, i.e., :math:`\lambda = \frac{1}{\mu}`.
mu :
Parameter of ``dist``. See above for details.
Returns
-------
A single numerical value, the approximated (weighted) hypervolume.
See Also
--------
hypervolume, hv_approx
References
----------
.. footbibliography::
Examples
--------
>>> moocore.hypervolume([[2, 2]], ref=4)
4.0
>>> moocore.whv_hype([[2, 2]], ref=4, ideal=1, seed=42)
3.99807
>>> moocore.hypervolume([[3, 1]], ref=4)
3.0
>>> moocore.whv_hype([[3, 1]], ref=4, ideal=1, seed=42)
3.00555
>>> moocore.whv_hype([[2, 2]], ref=4, ideal=1, dist="exponential", mu=0.2, seed=42)
1.14624
>>> moocore.whv_hype([[3, 1]], ref=4, ideal=1, dist="exponential", mu=0.2, seed=42)
1.66815
>>> moocore.whv_hype(
... [[2, 2]],
... ref=4,
... ideal=1,
... dist="point",
... mu=[2.9, 0.9],
... seed=42,
... )
0.64485
>>> moocore.whv_hype(
... [[3, 1]],
... ref=4,
... ideal=1,
... dist="point",
... mu=[2.9, 0.9],
... seed=42,
... )
4.03632
"""
# Convert to numpy.array in case the user provides a list. We use
# np.asarray to convert it to floating-point, otherwise if a user inputs
# something like [10, 10] then numpy would interpret it as an int array.
points, points_copied = asarray_maybe_copy(points)
nobj = points.shape[1]
if nobj != 2:
raise NotImplementedError("Only 2D points are currently supported")
ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref")
ideal = array_1d_of_length_n(
np.asarray(ideal, dtype=float), nobj, name="ideal"
)
maximise = _parse_maximise(maximise, nobj)
# FIXME: Do this in C.
if maximise.any():
if not points_copied:
points = points.copy()
points[:, maximise] = -points[:, maximise]
# These are so small that is ok to just copy them.
ref = ref.copy()
ref[maximise] = -ref[maximise]
ideal = ideal.copy()
ideal[maximise] = -ideal[maximise]
points_p, npoints, nobj = np2d_to_double_array(points)
ref = ffi.from_buffer("double []", ref)
ideal = ffi.from_buffer("double []", ideal)
# FIXME: Check ranges.
seed = _get_seed_for_c(seed)
match dist:
case "uniform":
hv = lib.whv_hype_unif(
points_p, npoints, ideal, ref, nsamples, seed
)
case "exponential":
if not isinstance(mu, numbers.Real):
raise ValueError("'mu' must be a single number")
hv = lib.whv_hype_expo(
points_p, npoints, ideal, ref, nsamples, seed, float(mu)
)
case "point":
mu = array_1d_of_length_n(
np.asarray(mu, dtype=float), nobj, name="mu"
)
mu, _ = np1d_to_double_array(mu)
hv = lib.whv_hype_gaus(
points_p, npoints, ideal, ref, nsamples, seed, mu
)
case _:
raise ValueError("Unknown value of dist = '{dist}'")
return hv
[docs]
def apply_within_sets(
x: ArrayLike, sets: ArrayLike, func: Callable[..., Any], **kwargs: Any
) -> np.ndarray:
"""Split ``x`` by row according to ``sets`` and apply ``func`` to each sub-array.
Parameters
----------
x :
2D array to be divided into sub-arrays.
sets :
A list or 1D array of length equal to the number of rows of ``x``. The values are used as-is to determine the groups and do not need to be sorted.
func :
A function that can take a 2D array as input. This function may return (1) a 2D array with the same number of rows as the input,
(2) a 1D array as long as the number of input rows,
(3) a scalar value, or
(4) a 2D array with a single row.
kwargs :
Additional keyword arguments to ``func``.
Returns
-------
An array whose shape depends on the output of ``func``. See Examples below.
See Also
--------
is_nondominated_within_sets, filter_dominated_within_sets
Examples
--------
>>> sets = np.array([3, 1, 2, 4, 2, 3, 1])
>>> x = np.arange(len(sets) * 2).reshape(-1, 2)
>>> x = np.hstack((x, sets.reshape(-1, 1)))
If ``func`` returns an array with the same number of rows as the input (case 1),
then the output is ordered in exactly the same way as the input.
>>> moocore.apply_within_sets(x, sets, lambda x: x)
array([[ 0, 1, 3],
[ 2, 3, 1],
[ 4, 5, 2],
[ 6, 7, 4],
[ 8, 9, 2],
[10, 11, 3],
[12, 13, 1]])
This is also the behavior if ``func`` returns a 1D array with one value per input row (case 2).
>>> moocore.apply_within_sets(x, sets, lambda x: x.sum(axis=1))
array([ 4, 6, 11, 17, 19, 24, 26])
If ``func`` returns a single scalar (case 3) or a 2D array with a single row (case 4),
then the order of the output is the order of the unique values as found in
``sets``, without sorting the unique values, which is what
:meth:`pandas.Series.unique` returns and NOT what :func:`numpy.unique`
returns.
>>> moocore.apply_within_sets(x, sets, lambda x: x.max())
array([11, 13, 9, 7])
>>> moocore.apply_within_sets(x, sets, lambda x: [x.max(axis=0)])
array([[10, 11, 3],
[12, 13, 1],
[ 8, 9, 2],
[ 6, 7, 4]])
In the previous example, ``func`` returns a 2D array with a single row. The
following will produce an error because it returns a 1D array, which is
interpreted as case 2, but the number of values does not match the number
of input rows.
>>> moocore.apply_within_sets(
... x, sets, lambda x: x.max(axis=0)
... ) # doctest: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
ValueError: `func` returned an array of length 3 but the input has length 2 for rows [0 5]
"""
x = np.asarray(x)
_, idx, inv = np.unique(sets, return_index=True, return_inverse=True)
# Remember the original position of each element of each set.
idx = [np.flatnonzero(inv == i) for i in idx.argsort()]
res = []
shorter = False
for g_idx in idx:
z = func(x.take(g_idx, axis=0), **kwargs)
z = np.atleast_1d(z)
if len(z) != len(g_idx):
if len(z) != 1:
raise ValueError(
f"`func` returned an array of length {len(z)} but the input has length {len(g_idx)} for rows {g_idx}"
)
shorter = True
res.append(z)
res = np.concatenate(res)
if not shorter:
res = res.take(np.concatenate(idx).argsort(), axis=0)
return res
[docs]
@DocSubstitute()
def r2_exact(
points: ArrayLike,
/,
ref: ArrayLike,
*,
maximise: bool | Sequence[bool | int] = False,
) -> float:
r"""Exact R2 indicator.
Computes the exact R2 indicator with respect to a given ideal/utopian reference point
assuming minimization of all objectives.
.. seealso:: For details of the R2 indicator, see :ref:`r2_indicator`.
.. warning::
The current implementation only supports 2 objectives.
Parameters
----------
points :
${points}
ref :
${ref_point}
maximise :
${maximise}
Returns
-------
A single numerical value, the exact R2 indicator.
Notes
-----
The current implementation exclusively supports bi-objective solution sets and runs in :math:`O(n \log n)`.
For more details on the computation of the exact R2 indicator refer to :footcite:t:`SchKer2025r2v2`.
References
----------
.. footbibliography::
Examples
--------
>>> dat = np.array([[5, 5], [4, 6], [2, 7], [7, 4]])
>>> moocore.r2_exact(dat, ref=[0, 0])
2.594191919191919
This function assumes that objectives must be minimized by default. We can
easily specify maximization:
>>> dat = np.array([[5, 5], [4, 6], [2, 7], [7, 4]])
>>> moocore.r2_exact(dat, ref=[10, 10], maximise=True)
2.519696969696969
Merge all the sets of a dataset by removing the set number column:
>>> dat = moocore.get_dataset("input1.dat")[:, :-1]
>>> len(dat)
100
Dominated points are ignored, so this:
>>> moocore.r2_exact(dat, ref=0)
0.3336076878950565
gives the same exact R2 value as this:
>>> dat = moocore.filter_dominated(dat)
>>> len(dat)
6
>>> moocore.r2_exact(dat, ref=0)
0.3336076878950565
"""
# Convert to numpy.array in case the user provides a list. We use
# np.asarray to convert it to floating-point, otherwise if a user inputs
# something like ref = np.array([10, 10]) then numpy would interpret it as
# an int array.
points, points_copied = asarray_maybe_copy(points)
nobj = points.shape[1]
if nobj != 2:
raise NotImplementedError("Only 2D points are currently supported")
# Make sure it is a 1D array of length nobj.
ref = array_1d_of_length_n(np.asarray(ref, dtype=float), nobj, name="ref")
maximise = _parse_maximise(maximise, nobj)
# FIXME: Do this in C.
if maximise.any():
if not points_copied:
points = points.copy()
points[:, maximise] = -points[:, maximise]
ref = ref.copy()
ref[maximise] = -ref[maximise]
points_p, npoints, nobj = np2d_to_double_array(
points, ctype_shape=("size_t", "uint_fast8_t")
)
ref = ffi.from_buffer("double []", ref)
r2 = lib.r2_exact(points_p, npoints, nobj, ref)
# if r2 < 0:
# raise MemoryError("memory allocation failed")
return r2