Pareto optimization

CHANGELOG.md

@@ -9,6 +9,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `BCUT2D` encoding for `SubstanceParameter`
 - Stored benchmarking results now include the Python environment and version
 - `qPSTD` acquisition function
+- `ParetoObjective` class for Pareto optimization of multiple targets and corresponding
+  `qLogNoisyExpectedHypervolumeImprovement` acquisition function
+- `BroadcastingSurrogate` class and corresponding `Surrogate.broadcast` method for
+  making single-target surrogate models multi-target compatible
+- `CompositeSurrogate` class for composing multi-target surrogates from single-target
+  surrogates
+- `supports_multi_output` attribute/property to `Surrogate`/`AcquisitionFunction`
 
 ### Changed
 - Acquisition function indicator `is_mc` has been removed in favor of new indicators 

README.md

@@ -37,7 +37,8 @@ The following provides a non-comprehensive overview:
 
 - 🛠️ Custom parameter encodings: Improve your campaign with domain knowledge
 - 🧪 Built-in chemical encodings: Improve your campaign with chemical knowledge
-- 🎯 Single and multiple targets with min, max and match objectives
+- 🎯 Numerical and binary targets with min, max and match objectives
+- ⚖️  Multi-target support via Pareto optimization and desirability scalarization
 - 🔍 Insights: Easily analyze feature importance and model behavior
 - 🎭 Hybrid (mixed continuous and discrete) spaces
 - 🚀 Transfer learning: Mix data from multiple campaigns and accelerate optimization
@@ -78,8 +79,8 @@ target = NumericalTarget(
 objective = SingleTargetObjective(target=target)
 ```
 In cases where we are confronted with multiple (potentially conflicting) targets,
-the `DesirabilityObjective` can be used instead. It allows to define additional
-settings, such as how these targets should be balanced.
+the `ParetoObjective` or `DesirabilityObjective` can be used instead.
+These allow to define additional settings, such as how the targets should be balanced.
 For more details, see the
 [objectives section](https://emdgroup.github.io/baybe/stable/userguide/objectives.html)
 of the user guide.

baybe/acquisition/__init__.py

@@ -10,6 +10,7 @@
     qExpectedImprovement,
     qKnowledgeGradient,
     qLogExpectedImprovement,
+    qLogNoisyExpectedHypervolumeImprovement,
     qLogNoisyExpectedImprovement,
     qNegIntegratedPosteriorVariance,
     qNoisyExpectedImprovement,
@@ -37,6 +38,7 @@
 UCB = UpperConfidenceBound
 qUCB = qUpperConfidenceBound
 qTS = qThompsonSampling
+qLogNEHVI = qLogNoisyExpectedHypervolumeImprovement
 
 __all__ = [
     ######################### Acquisition functions
@@ -64,6 +66,8 @@
     "qUpperConfidenceBound",
     # Thompson Sampling
     "qThompsonSampling",
+    # Hypervolume Improvement
+    "qLogNoisyExpectedHypervolumeImprovement",
     ######################### Abbreviations
     # Knowledge Gradient
     "qKG",
@@ -89,4 +93,6 @@
     "qUCB",
     # Thompson Sampling
     "qTS",
+    # Hypervolume Improvement
+    "qLogNEHVI",
 ]

baybe/acquisition/acqfs.py

@@ -4,6 +4,8 @@
 import math
 from typing import ClassVar
 
+import numpy as np
+import numpy.typing as npt
 import pandas as pd
 from attr.converters import optional as optional_c
 from attr.validators import optional as optional_v
@@ -13,7 +15,7 @@
 
 from baybe.acquisition.base import AcquisitionFunction
 from baybe.searchspace import SearchSpace
-from baybe.utils.basic import classproperty
+from baybe.utils.basic import classproperty, convert_to_float
 from baybe.utils.sampling_algorithms import (
     DiscreteSamplingMethod,
     sample_numerical_df,
@@ -320,5 +322,82 @@ def supports_batching(cls) -> bool:
         return False
 
 
+########################################################################################
+### Hypervolume Improvement
+@define(frozen=True)
+class qLogNoisyExpectedHypervolumeImprovement(AcquisitionFunction):
+    """Logarithmic Monte Carlo based noisy expected hypervolume improvement."""
+
+    abbreviation: ClassVar[str] = "qLogNEHVI"
+
+    ref_point: float | tuple[float, ...] | None = field(
+        default=None, converter=optional_c(convert_to_float)
+    )
+    """The reference point for computing the hypervolume improvement.
+
+    * When omitted, a default reference point is computed based on the provided data.
+    * When specified as a float, the value is interpreted as a multiplicative factor
+      determining the reference point location based on the difference between the best
+      and worst target configuration in the provided data.
+    * When specified as a vector, the input is taken as is.
+    """
+
+    prune_baseline: bool = field(default=True, validator=instance_of(bool))
+    """Auto-prune candidates that are unlikely to be the best."""
+
+    @staticmethod
+    def compute_ref_point(
+        array: npt.ArrayLike, maximize: npt.ArrayLike, factor: float = 0.1
+    ) -> np.ndarray:
+        """Compute a reference point for a given set of of target configurations.
+
+        The reference point is positioned in relation to the worst target configuration
+        within the provided array. The distance in each target dimension is adjusted by
+        a specified multiplication factor, which scales the reference point away from
+        the worst target configuration based on the maximum observed differences in
+        target values.
+
+        Example:
+            >>> from baybe.acquisition import qLogNEHVI
+
+            >>> qLogNEHVI.compute_ref_point([[0, 10], [2, 20]], [True, True], 0.1)
+            array([-0.2,  9. ])
+
+            >>> qLogNEHVI.compute_ref_point([[0, 10], [2, 20]], [True, False], 0.2)
+            array([ -0.4, 22. ])
+
+        Args:
+            array: A 2-D array-like where each row represents a target configuration.
+            maximize: A 1-D boolean array indicating which targets are to be maximized.
+            factor: A numeric value controlling the location of the reference point.
+
+        Raises:
+            ValueError: If the given target configuration array is not two-dimensional.
+            ValueError: If the given Boolean array is not one-dimensional.
+
+        Returns:
+            The computed reference point.
+        """
+        if np.ndim(array) != 2:
+            raise ValueError(
+                "The specified data array must have exactly two dimensions."
+            )
+        if np.ndim(maximize) != 1:
+            raise ValueError(
+                "The specified Boolean array must have exactly one dimension."
+            )
+
+        # Convert arrays
+        array = np.asarray(array)
+        maximize = np.where(maximize, 1.0, -1.0)
+
+        # Compute bounds
+        array = array * maximize[None, :]
+        min = np.min(array, axis=0)
+        max = np.max(array, axis=0)
+
+        return (min - factor * (max - min)) * maximize
+
+
 # Collect leftover original slotted classes processed by `attrs.define`
 gc.collect()

baybe/acquisition/base.py

@@ -5,6 +5,7 @@
 import gc
 import warnings
 from abc import ABC
+from collections.abc import Iterable
 from inspect import signature
 from typing import TYPE_CHECKING, ClassVar
 
@@ -17,6 +18,7 @@
 )
 from baybe.objectives.base import Objective
 from baybe.objectives.desirability import DesirabilityObjective
+from baybe.objectives.pareto import ParetoObjective
 from baybe.objectives.single import SingleTargetObjective
 from baybe.searchspace.core import SearchSpace
 from baybe.serialization.core import (
@@ -56,6 +58,11 @@ def supports_pending_experiments(cls) -> bool:
         """
         return cls.supports_batching
 
+    @classproperty
+    def supports_multi_output(cls) -> bool:
+        """Flag indicating whether multiple outputs are supported."""
+        return "Hypervolume" in cls.__name__  # type: ignore[attr-defined]
+
     @classproperty
     def _non_botorch_attrs(cls) -> tuple[str, ...]:
         """Names of attributes that are not passed to the BoTorch constructor."""
@@ -76,9 +83,13 @@ def to_botorch(
         """
         import botorch.acquisition as bo_acqf
         import torch
+        from botorch.acquisition.multi_objective import WeightedMCMultiOutputObjective
         from botorch.acquisition.objective import LinearMCObjective
 
-        from baybe.acquisition.acqfs import qThompsonSampling
+        from baybe.acquisition.acqfs import (
+            qLogNoisyExpectedHypervolumeImprovement,
+            qThompsonSampling,
+        )
 
         # Retrieve botorch acquisition function class and match attributes
         acqf_cls = _get_botorch_acqf_class(type(self))
@@ -151,6 +162,39 @@ def to_botorch(
                     additional_params["best_f"] = (
                         bo_surrogate.posterior(train_x).mean.max().item()
                     )
+            case ParetoObjective():
+                if not isinstance(self, qLogNoisyExpectedHypervolumeImprovement):
+                    raise IncompatibleAcquisitionFunctionError(
+                        f"Pareto optimization currently supports the "
+                        f"'{qLogNoisyExpectedHypervolumeImprovement.__name__}' "
+                        f"acquisition function only."
+                    )
+                if not all(
+                    isinstance(t, NumericalTarget)
+                    and t.mode in (TargetMode.MAX, TargetMode.MIN)
+                    for t in objective.targets
+                ):
+                    raise NotImplementedError(
+                        "Pareto optimization currently supports "
+                        "maximization/minimization targets only."
+                    )
+                maximize = [t.mode is TargetMode.MAX for t in objective.targets]  # type: ignore[attr-defined]
+                multiplier = [1.0 if m else -1.0 for m in maximize]
+                additional_params["objective"] = WeightedMCMultiOutputObjective(
+                    torch.tensor(multiplier)
+                )
+                train_y = measurements[[t.name for t in objective.targets]].to_numpy()
+                if isinstance(ref_point := params_dict["ref_point"], Iterable):
+                    ref_point = [
+                        p * m for p, m in zip(ref_point, multiplier, strict=True)
+                    ]
+                else:
+                    kwargs = {"factor": ref_point} if ref_point is not None else {}
+                    ref_point = (
+                        self.compute_ref_point(train_y, maximize, **kwargs) * multiplier
+                    )
+                params_dict["ref_point"] = ref_point
+
             case _:
                 raise ValueError(f"Unsupported objective type: {objective}")
 
@@ -172,7 +216,9 @@ def _get_botorch_acqf_class(
     import botorch
 
     for cls in baybe_acqf_cls.mro():
-        if acqf_cls := getattr(botorch.acquisition, cls.__name__, False):
+        if acqf_cls := getattr(botorch.acquisition, cls.__name__, False) or getattr(
+            botorch.acquisition.multi_objective, cls.__name__, False
+        ):
             if is_abstract(acqf_cls):
                 continue
             return acqf_cls  # type: ignore

baybe/exceptions.py

@@ -25,6 +25,10 @@ class IncompatibleSearchSpaceError(IncompatibilityError):
     """
 
 
+class IncompatibleSurrogateError(IncompatibilityError):
+    """An incompatible surrogate was selected."""
+
+
 class IncompatibleAcquisitionFunctionError(IncompatibilityError):
     """An incompatible acquisition function was selected."""
 

baybe/objectives/__init__.py

@@ -1,9 +1,11 @@
 """BayBE objectives."""
 
 from baybe.objectives.desirability import DesirabilityObjective
+from baybe.objectives.pareto import ParetoObjective
 from baybe.objectives.single import SingleTargetObjective
 
 __all__ = [
     "SingleTargetObjective",
     "DesirabilityObjective",
+    "ParetoObjective",
 ]

baybe/objectives/desirability.py

@@ -18,7 +18,7 @@
 from baybe.targets.base import Target
 from baybe.targets.numerical import NumericalTarget
 from baybe.utils.basic import is_all_instance, to_tuple
-from baybe.utils.dataframe import get_transform_objects, pretty_print_df
+from baybe.utils.dataframe import pretty_print_df, transform_target_columns
 from baybe.utils.numerical import geom_mean
 from baybe.utils.plotting import to_string
 from baybe.utils.validation import finite_float
@@ -145,7 +145,7 @@ def transform(
         # >>>>>>>>>> Deprecation
         if not ((df is None) ^ (data is None)):
             raise ValueError(
-                "Provide the dataframe to be transformed as argument to `df`."
+                "Provide the dataframe to be transformed as first positional argument."
             )
 
         if data is not None:
@@ -172,15 +172,10 @@ def transform(
                 )
         # <<<<<<<<<< Deprecation
 
-        # Extract the relevant part of the dataframe
-        targets = get_transform_objects(
+        # Transform all targets individually
+        transformed = transform_target_columns(
             df, self.targets, allow_missing=allow_missing, allow_extra=allow_extra
         )
-        transformed = df[[t.name for t in targets]].copy()
-
-        # Transform all targets individually
-        for target in self.targets:
-            transformed[target.name] = target.transform(df[target.name])
 
         # Scalarize the transformed targets into desirability values
         vals = scalarize(transformed.values, self.scalarizer, self._normalized_weights)

baybe/objectives/pareto.py

@@ -0,0 +1,74 @@
+"""Functionality for multi-target objectives."""
+
+import warnings
+
+import pandas as pd
+from attrs import define, field
+from attrs.validators import deep_iterable, instance_of, min_len
+from typing_extensions import override
+
+from baybe.objectives.base import Objective
+from baybe.targets.base import Target
+from baybe.utils.basic import to_tuple
+from baybe.utils.dataframe import transform_target_columns
+
+
+@define(frozen=True, slots=False)
+class ParetoObjective(Objective):
+    """An objective handling multiple targets in a Pareto sense."""
+
+    _targets: tuple[Target, ...] = field(
+        converter=to_tuple,
+        validator=[min_len(2), deep_iterable(member_validator=instance_of(Target))],
+        alias="targets",
+    )
+    "The targets considered by the objective."
+
+    @override
+    @property
+    def targets(self) -> tuple[Target, ...]:
+        return self._targets
+
+    @override
+    def transform(
+        self,
+        df: pd.DataFrame | None = None,
+        /,
+        *,
+        allow_missing: bool = False,
+        allow_extra: bool | None = None,
+        data: pd.DataFrame | None = None,
+    ) -> pd.DataFrame:
+        # >>>>>>>>>> Deprecation
+        if not ((df is None) ^ (data is None)):
+            raise ValueError(
+                "Provide the dataframe to be transformed as first positional argument."
+            )
+
+        if data is not None:
+            df = data
+            warnings.warn(
+                "Providing the dataframe via the `data` argument is deprecated and "
+                "will be removed in a future version. Please pass your dataframe "
+                "as positional argument instead.",
+                DeprecationWarning,
+            )
+
+        # Mypy does not infer from the above that `df` must be a dataframe here
+        assert isinstance(df, pd.DataFrame)
+
+        if allow_extra is None:
+            allow_extra = True
+            if set(df.columns) - {p.name for p in self.targets}:
+                warnings.warn(
+                    "For backward compatibility, the new `allow_extra` flag is set "
+                    "to `True` when left unspecified. However, this behavior will be "
+                    "changed in a future version. If you want to invoke the old "
+                    "behavior, please explicitly set `allow_extra=True`.",
+                    DeprecationWarning,
+                )
+        # <<<<<<<<<< Deprecation
+
+        return transform_target_columns(
+            df, self.targets, allow_missing=allow_missing, allow_extra=allow_extra
+        )