Merge pull request #344 from till-m/constrained-optimization

Advanced constrained optimization

Author: bwheelz36

bayes_opt/__init__.py

@@ -2,9 +2,11 @@
 from .domain_reduction import SequentialDomainReductionTransformer
 from .util import UtilityFunction
 from .logger import ScreenLogger, JSONLogger
+from .constraint import ConstraintModel
 
 __all__ = [
     "BayesianOptimization",
+    "ConstraintModel"
     "UtilityFunction",
     "Events",
     "ScreenLogger",

bayes_opt/bayesian_optimization.py

@@ -1,7 +1,9 @@
 import warnings
 from queue import Queue, Empty
 
-from .target_space import TargetSpace
+from bayes_opt.constraint import ConstraintModel
+
+from .target_space import TargetSpace, ConstrainedTargetSpace
 from .event import Events, DEFAULT_EVENTS
 from .logger import _get_default_logger
 from .util import UtilityFunction, acq_max, ensure_rng
@@ -16,6 +18,7 @@ class Observable(object):
     Inspired/Taken from
         https://www.protechtraining.com/blog/post/879#simple-observer
     """
+
     def __init__(self, events):
         # maps event names to subscribers
         # str -> dict
@@ -52,9 +55,12 @@ class BayesianOptimization(Observable):
         Dictionary with parameters names as keys and a tuple with minimum
         and maximum values.
 
+    constraint: A ConstraintModel. Note that the names of arguments of the
+        constraint function and of f need to be the same.
+
     random_state: int or numpy.random.RandomState, optional(default=None)
         If the value is an integer, it is used as the seed for creating a
-        numpy.random.RandomState. Otherwise the random state provided it is used.
+        numpy.random.RandomState. Otherwise the random state provided is used.
         When set to None, an unseeded random state is generated.
 
     verbose: int, optional(default=2)
@@ -76,14 +82,16 @@ class BayesianOptimization(Observable):
     set_bounds()
         Allows changing the lower and upper searching bounds
     """
-    def __init__(self, f, pbounds, random_state=None, verbose=2,
+
+    def __init__(self,
+                 f,
+                 pbounds,
+                 constraint=None,
+                 random_state=None,
+                 verbose=2,
                  bounds_transformer=None):
         self._random_state = ensure_rng(random_state)
 
-        # Data structure containing the function to be optimized, the bounds of
-        # its domain, and a record of the evaluations we have done so far
-        self._space = TargetSpace(f, pbounds, random_state)
-
         self._queue = Queue()
 
         # Internal GP regressor
@@ -95,6 +103,27 @@ def __init__(self, f, pbounds, random_state=None, verbose=2,
             random_state=self._random_state,
         )
 
+        if constraint is None:
+            # Data structure containing the function to be optimized, the
+            # bounds of its domain, and a record of the evaluations we have
+            # done so far
+            self._space = TargetSpace(f, pbounds, random_state)
+            self.is_constrained = False
+        else:
+            constraint_ = ConstraintModel(
+                constraint.fun,
+                constraint.lb,
+                constraint.ub,
+                random_state=random_state
+            )
+            self._space = ConstrainedTargetSpace(
+                f,
+                constraint_,
+                pbounds,
+                random_state
+            )
+            self.is_constrained = True
+
         self._verbose = verbose
         self._bounds_transformer = bounds_transformer
         if self._bounds_transformer:
@@ -110,6 +139,12 @@ def __init__(self, f, pbounds, random_state=None, verbose=2,
     def space(self):
         return self._space
 
+    @property
+    def constraint(self):
+        if self.is_constrained:
+            return self._space.constraint
+        return None
+
     @property
     def max(self):
         return self._space.max()
@@ -136,6 +171,7 @@ def probe(self, params, lazy=True):
             If True, the optimizer will evaluate the points when calling
             maximize(). Otherwise it will evaluate it at the moment.
         """
+
         if lazy:
             self._queue.put(params)
         else:
@@ -152,15 +188,17 @@ def suggest(self, utility_function):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             self._gp.fit(self._space.params, self._space.target)
+            if self.is_constrained:
+                self.constraint.fit(self._space.params,
+                                    self._space._constraint_values)
 
         # Finding argmax of the acquisition function.
-        suggestion = acq_max(
-            ac=utility_function.utility,
-            gp=self._gp,
-            y_max=self._space.target.max(),
-            bounds=self._space.bounds,
-            random_state=self._random_state
-        )
+        suggestion = acq_max(ac=utility_function.utility,
+                             gp=self._gp,
+                             constraint=self.constraint,
+                             y_max=self._space.target.max(),
+                             bounds=self._space.bounds,
+                             random_state=self._random_state)
 
         return self._space.array_to_params(suggestion)
 
@@ -211,15 +249,15 @@ def maximize(self,
         kappa: float, optional(default=2.576)
             Parameter to indicate how closed are the next parameters sampled.
                 Higher value = favors spaces that are least explored.
-                Lower value = favors spaces where the regression function is the
-                highest.
+                Lower value = favors spaces where the regression function is
+                the highest.
 
         kappa_decay: float, optional(default=1)
             `kappa` is multiplied by this factor every iteration.
 
         kappa_decay_delay: int, optional(default=0)
-            Number of iterations that must have passed before applying the decay
-            to `kappa`.
+            Number of iterations that must have passed before applying the
+            decay to `kappa`.
 
         xi: float, optional(default=0.0)
             [unused]
@@ -242,12 +280,11 @@ def maximize(self,
                 util.update_params()
                 x_probe = self.suggest(util)
                 iteration += 1
-
             self.probe(x_probe, lazy=False)
 
             if self._bounds_transformer and iteration > 0:
-                # The bounds transformer should only modify the bounds after the init_points points (only for the true
-                # iterations)
+                # The bounds transformer should only modify the bounds after
+                # the init_points points (only for the true iterations)
                 self.set_bounds(
                     self._bounds_transformer.transform(self._space))
 

bayes_opt/constraint.py

@@ -0,0 +1,154 @@
+import numpy as np
+from sklearn.gaussian_process.kernels import Matern
+from sklearn.gaussian_process import GaussianProcessRegressor
+from scipy.stats import norm
+
+
+class ConstraintModel():
+    """
+    This class takes the function to optimize as well as the parameters bounds
+    in order to find which values for the parameters yield the maximum value
+    using bayesian optimization.
+
+    Parameters
+    ----------
+    fun: function
+        Constraint function. If multiple constraints are handled, this should
+        return a numpy.ndarray of appropriate size.
+
+    lb: numeric or numpy.ndarray
+        Upper limit(s) for the constraints. The return value of `fun` should
+        have exactly this shape.
+
+    ub: numeric or numpy.ndarray
+        Upper limit(s) for the constraints. The return value of `fun` should
+        have exactly this shape.
+
+    random_state: int or numpy.random.RandomState, optional(default=None)
+        If the value is an integer, it is used as the seed for creating a
+        numpy.random.RandomState. Otherwise the random state provided is used.
+        When set to None, an unseeded random state is generated.
+
+    Note
+    ----
+    In case of multiple constraints, this model assumes conditional
+    independence. This means that for each constraint, the probability of
+    fulfillment is the cdf of a univariate Gaussian. The overall probability
+    is a simply the product of the individual probabilities.
+    """
+
+    def __init__(self, fun, lb, ub, random_state=None):
+        self.fun = fun
+
+        if isinstance(lb, float):
+            self._lb = np.array([lb])
+        else:
+            self._lb = lb
+        
+        if isinstance(ub, float):
+            self._ub = np.array([ub])
+        else:
+            self._ub = ub
+        
+        
+        basis = lambda: GaussianProcessRegressor(
+            kernel=Matern(nu=2.5),
+            alpha=1e-6,
+            normalize_y=True,
+            n_restarts_optimizer=5,
+            random_state=random_state,
+        )
+        self._model = [basis() for _ in range(len(self._lb))]
+
+    @property
+    def lb(self):
+        return self._lb
+
+    @property
+    def ub(self):
+        return self._ub
+    
+    @property
+    def model(self):
+        return self._model
+
+    def eval(self, **kwargs):
+        """
+        Evaluates the constraint function.
+        """
+        try:
+            return self.fun(**kwargs)
+        except TypeError as e:
+            msg = (
+                "Encountered TypeError when evaluating constraint " +
+                "function. This could be because your constraint function " +
+                "doesn't use the same keyword arguments as the target " +
+                f"function. Original error message:\n\n{e}"
+                )
+            e.args = (msg,)
+            raise
+
+    def fit(self, X, Y):
+        """
+        Fits internal GaussianProcessRegressor's to the data.
+        """
+        if len(self._model) == 1:
+            self._model[0].fit(X, Y)
+        else:
+            for i, gp in enumerate(self._model):
+                gp.fit(X, Y[:, i])
+
+    def predict(self, X):
+        """
+        Returns the probability that the constraint is fulfilled at `X` based
+        on the internal Gaussian Process Regressors.
+
+        Note that this does not try to approximate the values of the constraint
+        function, but probability that the constraint function is fulfilled.
+        For the former, see `ConstraintModel.approx()`.
+        """
+        X_shape = X.shape
+        X = X.reshape((-1, self._model[0].n_features_in_))
+        if len(self._model) == 1:
+            y_mean, y_std = self._model[0].predict(X, return_std=True)
+
+            p_lower = (norm(loc=y_mean, scale=y_std).cdf(self._lb[0])
+                            if self._lb[0] != -np.inf else np.array([0]))
+            p_upper = (norm(loc=y_mean, scale=y_std).cdf(self._ub[0])
+                            if self._lb[0] != np.inf else np.array([1]))
+            result = p_upper - p_lower
+            return result.reshape(X_shape[:-1])
+        else:
+            result = np.ones(X.shape[0])
+            for j, gp in enumerate(self._model):
+                y_mean, y_std = gp.predict(X, return_std=True)
+                p_lower = (norm(loc=y_mean, scale=y_std).cdf(self._lb[j])
+                           if self._lb[j] != -np.inf else np.array([0]))
+                p_upper = (norm(loc=y_mean, scale=y_std).cdf(self._ub[j])
+                           if self._lb[j] != np.inf else np.array([1]))
+                result = result * (p_upper - p_lower)
+            return result.reshape(X_shape[:-1])
+
+    def approx(self, X):
+        """
+        Returns the approximation of the constraint function using the internal
+        Gaussian Process Regressors.
+        """
+        X_shape = X.shape
+        X = X.reshape((-1, self._model[0].n_features_in_))
+        if len(self._model) == 1:
+            return self._model[0].predict(X).reshape(X_shape[:-1])
+        else:
+            result = np.column_stack([gp.predict(X) for gp in self._model])
+            return result.reshape(X_shape[:-1] + (len(self._lb), ))
+
+    def allowed(self, constraint_values):
+        """
+        Checks whether `constraint_values` are below the specified limits.
+        """
+        if self._lb.size == 1:
+            return (np.less_equal(self._lb, constraint_values)
+                    & np.less_equal(constraint_values, self._ub))
+
+        return (np.all(constraint_values <= self._ub, axis=-1)
+                    & np.all(constraint_values >= self._lb, axis=-1))

bayes_opt/logger.py

@@ -6,7 +6,6 @@
 from .event import Events
 from .util import Colours
 
-
 def _get_default_logger(verbose):
     return ScreenLogger(verbose=verbose)
 
@@ -81,6 +80,11 @@ def _header(self, instance):
         return line + "\n" + ("-" * self._header_length)
 
     def _is_new_max(self, instance):
+        if instance.max["target"] is None:
+            # During constrained optimization, there might not be a maximum
+            # value since the optimizer might've not encountered any points
+            # that fulfill the constraints.
+            return False
         if self._previous_max is None:
             self._previous_max = instance.max["target"]
         return instance.max["target"] > self._previous_max

bayes_opt/observer.py

@@ -25,8 +25,8 @@ def _update_tracker(self, event, instance):
             self._iterations += 1
 
             current_max = instance.max
-            if (self._previous_max is None or
-                current_max["target"] > self._previous_max):
+            if (self._previous_max is None
+                    or current_max["target"] > self._previous_max):
                 self._previous_max = current_max["target"]
                 self._previous_max_params = current_max["params"]