Mobo enhancement

xopt/generators/bayesian/bayesian_generator.py

@@ -4,13 +4,14 @@
 import warnings
 from abc import ABC, abstractmethod
 from copy import deepcopy
-from typing import Dict, List, Optional
+from typing import Dict, List, Optional, Union
 
 import pandas as pd
 import torch
 from botorch.acquisition import FixedFeatureAcquisitionFunction, qUpperConfidenceBound
 from botorch.models.model import Model
 from botorch.sampling import get_sampler
+from botorch.utils.multi_objective import is_non_dominated
 from botorch.utils.multi_objective.box_decompositions import DominatedPartitioning
 from gpytorch import Module
 from pydantic import Field, field_validator, PositiveInt, SerializeAsAny
@@ -407,8 +408,18 @@ def propose_candidates(self, model, n_candidates=1):
         # get acquisition function
         acq_funct = self.get_acquisition(model)
 
-        # get candidates
-        candidates = self.numerical_optimizer.optimize(acq_funct, bounds, n_candidates)
+        # get initial candidates to start acquisition function optimization
+        initial_points = self._get_initial_conditions(n_candidates)
+
+        # get candidates -- grid optimizer does not support batch_initial_conditions
+        if isinstance(self.numerical_optimizer, GridOptimizer):
+            candidates = self.numerical_optimizer.optimize(
+                acq_funct, bounds, n_candidates
+            )
+        else:
+            candidates = self.numerical_optimizer.optimize(
+                acq_funct, bounds, n_candidates, batch_initial_conditions=initial_points
+            )
         return candidates
 
     def get_training_data(self, data: pd.DataFrame) -> pd.DataFrame:
@@ -536,6 +547,11 @@ def visualize_model(self, **kwargs):
         """displays the GP models"""
         return visualize_generator_model(self, **kwargs)
 
+    def _get_initial_conditions(self, n_candidates=1) -> Union[Tensor, None]:
+        """overwrite if algorithm should specifiy initial candidates for optimizing
+        the acquisition function"""
+        return None
+
     def _process_candidates(self, candidates: Tensor):
         """process pytorch candidates from optimizing the acquisition function"""
         logger.debug(f"Best candidate from optimize {candidates}")
@@ -776,31 +792,52 @@ def torch_reference_point(self):
 
         return torch.tensor(pt, **self._tkwargs)
 
-    def calculate_hypervolume(self):
-        """compute hypervolume given data"""
-        objective_data = torch.tensor(
-            self.vocs.objective_data(self.data, return_raw=True).to_numpy()
+    def _get_scaled_data(self):
+        """get scaled input/objective data for use with botorch logic which assumes
+        maximization for each objective"""
+        var_df, obj_df, _, _ = self.vocs.extract_data(
+            self.data, return_valid=True, return_raw=True
         )
 
-        # hypervolume must only take into account feasible data points
-        if self.vocs.n_constraints > 0:
-            objective_data = objective_data[
-                self.vocs.feasibility_data(self.data)["feasible"].to_list()
-            ]
+        variable_data = torch.tensor(var_df[self.vocs.variable_names].to_numpy())
+        objective_data = torch.tensor(obj_df[self.vocs.objective_names].to_numpy())
+        weights = set_botorch_weights(self.vocs).to(**self._tkwargs)[
+            : self.vocs.n_objectives
+        ]
+        return variable_data, objective_data * weights, weights
 
-        n_objectives = self.vocs.n_objectives
-        weights = torch.zeros(n_objectives)
-        weights = set_botorch_weights(self.vocs).to(**self._tkwargs)
-        objective_data = objective_data * weights
+    def calculate_hypervolume(self):
+        """compute hypervolume given data"""
 
         # compute hypervolume
         bd = DominatedPartitioning(
-            ref_point=self.torch_reference_point, Y=objective_data
+            ref_point=self.torch_reference_point, Y=self._get_scaled_data()[1]
         )
         volume = bd.compute_hypervolume().item()
 
         return volume
 
+    def get_pareto_front(self):
+        """compute the pareto front x/y values given data"""
+        variable_data, objective_data, weights = self._get_scaled_data()
+        obj_data = torch.vstack(
+            (self.torch_reference_point.unsqueeze(0), objective_data)
+        )
+        var_data = torch.vstack(
+            (
+                torch.full_like(variable_data[0], float("Nan")).unsqueeze(0),
+                variable_data,
+            )
+        )
+        non_dominated = is_non_dominated(obj_data)
+
+        # note need to undo weights for real number output
+        # only return values if non nan values exist
+        if torch.all(torch.isnan(var_data[non_dominated])):
+            return None, None
+        else:
+            return var_data[non_dominated], obj_data[non_dominated] / weights
+
 
 def formatted_base_docstring():
     return "\nBase Generator\n---------------\n" + BayesianGenerator.__doc__

xopt/generators/bayesian/mobo.py

@@ -1,16 +1,29 @@
+from typing import Union
+
+import torch
 from botorch.acquisition import FixedFeatureAcquisitionFunction
 from botorch.acquisition.multi_objective import qNoisyExpectedHypervolumeImprovement
 from botorch.acquisition.multi_objective.logei import (
     qLogNoisyExpectedHypervolumeImprovement,
 )
+from botorch.utils import draw_sobol_samples
+from pydantic import Field
+from torch import Tensor
 
 from xopt.generators.bayesian.bayesian_generator import MultiObjectiveBayesianGenerator
 
 from xopt.generators.bayesian.objectives import create_mobo_objective
+from xopt.numerical_optimizer import LBFGSOptimizer
 
 
 class MOBOGenerator(MultiObjectiveBayesianGenerator):
     name = "mobo"
+    supports_batch_generation: bool = True
+    use_pf_as_initial_points: bool = Field(
+        False,
+        description="flag to specify if pareto front points are to be used during "
+        "optimization of the acquisition function",
+    )
     __doc__ = """Implements Multi-Objective Bayesian Optimization using the Expected
             Hypervolume Improvement acquisition function"""
 
@@ -64,3 +77,59 @@ def _get_acquisition(self, model):
         )
 
         return acq
+
+    def _get_initial_conditions(self, n_candidates=1) -> Union[Tensor, None]:
+        """
+        generate initial candidates for optimizing the acquisition function based on
+        the pareto front
+
+        If `use_pf_as_initial_points` flag is set to true then the current
+        Pareto-optimal set is used as initial points for optimizing the acquisition
+        function instead of randomly selected points (random points fill in the set
+        if `num_restarts` is greater than the number of points in the Pareto set).
+
+        Returns:
+            A `num_restarts x q x d` tensor of initial conditions.
+
+        """
+        if self.use_pf_as_initial_points:
+            if isinstance(self.numerical_optimizer, LBFGSOptimizer):
+                bounds = self._get_optimization_bounds()
+                num_restarts = self.numerical_optimizer.n_restarts
+
+                pf_locations, _ = self.get_pareto_front()
+
+                # if there is no pareto front just return None to revert back to
+                # default behavior
+                if pf_locations is None:
+                    return None
+
+                initial_points = torch.clone(pf_locations)
+
+                # add the q dimension
+                initial_points = initial_points.unsqueeze(1)
+                initial_points = initial_points.expand([-1, n_candidates, -1])
+
+                # initial_points must equal the number of restarts
+                if len(initial_points) < num_restarts:
+                    # add random points to the list inside the bounds
+                    sobol_samples = draw_sobol_samples(
+                        bounds, num_restarts - len(initial_points), n_candidates
+                    )
+
+                    initial_points = torch.cat([initial_points, sobol_samples])
+                elif len(initial_points) > num_restarts:
+                    # if there are too many select the first `num_restarts` points at
+                    # random
+                    initial_points = initial_points[
+                        torch.randperm(len(initial_points))
+                    ][:num_restarts]
+
+                return initial_points
+            else:
+                raise RuntimeWarning(
+                    "cannot use PF as initial optimization points "
+                    "for non-LBFGS optimizers, ignoring flag"
+                )
+
+        return None

xopt/numerical_optimizer.py

@@ -15,7 +15,9 @@ class NumericalOptimizer(XoptBaseModel, ABC):
     model_config = ConfigDict(extra="forbid")
 
     @abstractmethod
-    def optimize(self, function: AcquisitionFunction, bounds: Tensor, n_candidates=1):
+    def optimize(
+        self, function: AcquisitionFunction, bounds: Tensor, n_candidates=1, **kwargs
+    ):
         """optimize a function to produce a number of candidate points that
         minimize the function"""
         pass
@@ -33,7 +35,7 @@ class LBFGSOptimizer(NumericalOptimizer):
 
     model_config = ConfigDict(validate_assignment=True)
 
-    def optimize(self, function, bounds, n_candidates=1):
+    def optimize(self, function, bounds, n_candidates=1, **kwargs):
         assert isinstance(bounds, Tensor)
         if len(bounds) != 2:
             raise ValueError("bounds must have the shape [2, ndim]")
@@ -46,6 +48,7 @@ def optimize(self, function, bounds, n_candidates=1):
             num_restarts=self.n_restarts,
             timeout_sec=self.max_time,
             options={"maxiter": self.max_iter},
+            **kwargs,
         )
         return candidates
 

xopt/resources/test_functions/zdt.py

@@ -0,0 +1,51 @@
+import numpy as np
+
+from xopt import VOCS
+
+
+def construct_zdt(n_dims, problem_index=1):
+    """construct Xopt versions of the multiobjective ZDT test problems"""
+    vocs = VOCS(
+        variables={f"x{i + 1}": [0, 1] for i in range(n_dims)},
+        objectives={"f1": "MINIMIZE", "f2": "MINIMIZE"},
+    )
+
+    if problem_index == 1:
+        # ZDT1
+        def evaluate(input_dict):
+            x = np.array([input_dict[f"x{i + 1}"] for i in range(n_dims)])
+
+            f1 = x[0]
+            g = 1 + (9 / (n_dims - 1)) * np.sum(x[1:])
+            h = 1 - np.sqrt(f1 / g)
+            f2 = g * h
+
+            return {"f1": f1, "f2": f2, "g": g}
+    elif problem_index == 2:
+
+        def evaluate(input_dict):
+            x = np.array([input_dict[f"x{i + 1}"] for i in range(n_dims)])
+
+            f1 = x[0]
+            g = 1 + (9 / (n_dims - 1)) * np.sum(x[1:])
+            h = 1 - (f1 / g) ** 2
+            f2 = g * h
+
+            return {"f1": f1, "f2": f2, "g": g}
+    elif problem_index == 3:
+
+        def evaluate(input_dict):
+            x = np.array([input_dict[f"x{i + 1}"] for i in range(n_dims)])
+
+            f1 = x[0]
+            g = 1 + (9 / (n_dims - 1)) * np.sum(x[1:])
+            h = 1 - np.sqrt(f1 / g) - (f1 / g) * np.sin(10 * np.pi * f1)
+            f2 = g * h
+
+            return {"f1": f1, "f2": f2, "g": g}
+    else:
+        raise NotImplementedError()
+
+    reference_point = {"f1": 1.0, "f2": 1.0}
+
+    return vocs, evaluate, reference_point