Use scipy NonLinearConstraint, change the way ConstraintModel is constructed and stored

Author: till-m

bayes_opt/bayesian_optimization.py

@@ -1,6 +1,8 @@
 import warnings
 from queue import Queue, Empty
 
+from bayes_opt.constraint import ConstraintModel
+
 from .target_space import TargetSpace, ConstrainedTargetSpace
 from .event import Events, DEFAULT_EVENTS
 from .logger import _get_default_logger
@@ -106,10 +108,21 @@ def __init__(self,
             # 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:
-            self._space = ConstrainedTargetSpace(f, constraint, pbounds,
-                                                 random_state)
-        self.constraint = constraint
+            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
@@ -126,6 +139,12 @@ def __init__(self,
     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()
@@ -169,7 +188,7 @@ def suggest(self, utility_function):
         with warnings.catch_warnings():
             warnings.simplefilter("ignore")
             self._gp.fit(self._space.params, self._space.target)
-            if self.constraint is not None:
+            if self.is_constrained:
                 self.constraint.fit(self._space.params,
                                     self._space._constraint_values)
 

bayes_opt/constraint.py

@@ -12,12 +12,16 @@ class ConstraintModel():
 
     Parameters
     ----------
-    func: function
+    fun: function
         Constraint function. If multiple constraints are handled, this should
         return a numpy.ndarray of appropriate size.
 
-    limits: numeric or numpy.ndarray
-        Upper limit(s) for the constraints. The return value of `func` should
+    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)
@@ -33,41 +37,56 @@ class ConstraintModel():
     is a simply the product of the individual probabilities.
     """
 
-    def __init__(self, func, limits, random_state=None):
-        self.func = func
+    def __init__(self, fun, lb, ub, random_state=None):
+        self.fun = fun
 
-        if isinstance(limits, float):
-            self._limits = np.array([limits])
+        if isinstance(lb, float):
+            self._lb = np.array([lb])
         else:
-            self._limits = limits
-
+            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._limits))]
+        self._model = [basis() for _ in range(len(self._lb))]
 
     @property
-    def limits(self):
-        return self._limits
+    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.func(**kwargs)
+            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}"
                 )
-            raise TypeError(msg)
+            e.args = (msg,)
+            raise
 
     def fit(self, X, Y):
         """
@@ -92,14 +111,22 @@ def predict(self, X):
         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)
-            result = norm(loc=y_mean, scale=y_std).cdf(self._limits[0])
+
+            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)
-                result = result * norm(loc=y_mean, scale=y_std).cdf(
-                    self._limits[j])
+                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):
@@ -113,13 +140,15 @@ def approx(self, X):
             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._limits), ))
+            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._limits.size == 1:
-            return np.less_equal(constraint_values, self._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._limits, axis=-1)
+        return (np.all(constraint_values <= self._ub, axis=-1)
+                    & np.all(constraint_values >= self._lb, axis=-1))

bayes_opt/target_space.py

@@ -267,10 +267,14 @@ def __init__(self,
         self._constraint = constraint
 
         # preallocated memory for constraint fulfillement
-        if constraint.limits.size == 1:
+        if constraint.lb.size == 1:
             self._constraint_values = np.empty(shape=(0), dtype=float)
         else:
-            self._constraint_values = np.empty(shape=(0, constraint.limits.size), dtype=float)
+            self._constraint_values = np.empty(shape=(0, constraint.lb.size), dtype=float)
+
+    @property
+    def constraint(self):
+        return self._constraint
 
     @property
     def constraint_values(self):

examples/constraints.ipynb

@@ -1,26 +1,27 @@
 import numpy as np
 from bayes_opt import BayesianOptimization, ConstraintModel
 from pytest import approx, raises
+from scipy.optimize import NonlinearConstraint
 
 np.random.seed(42)
 
 
-def test_single_constraint():
+def test_single_constraint_upper():
 
     def target_function(x, y):
         return np.cos(2 * x) * np.cos(y) + np.sin(x)
 
     def constraint_function(x, y):
         return np.cos(x) * np.cos(y) - np.sin(x) * np.sin(y)
 
-    constraint_limit = 0.5
+    constraint_limit_upper = 0.5
 
-    conmod = ConstraintModel(constraint_function, constraint_limit)
+    constraint = NonlinearConstraint(constraint_function, -np.inf, constraint_limit_upper)
     pbounds = {'x': (0, 6), 'y': (0, 6)}
 
     optimizer = BayesianOptimization(
         f=target_function,
-        constraint=conmod,
+        constraint=constraint,
         pbounds=pbounds,
         verbose=0,
         random_state=1,
@@ -31,23 +32,25 @@ def constraint_function(x, y):
         n_iter=10,
     )
 
+    assert constraint_function(**optimizer.max["params"]) <= constraint_limit_upper
+
 
-def test_single_constraint_max_is_allowed():
+def test_single_constraint_lower():
 
     def target_function(x, y):
         return np.cos(2 * x) * np.cos(y) + np.sin(x)
 
     def constraint_function(x, y):
         return np.cos(x) * np.cos(y) - np.sin(x) * np.sin(y)
 
-    constraint_limit = 0.5
+    constraint_limit_lower = -0.5
 
-    conmod = ConstraintModel(constraint_function, constraint_limit)
+    constraint = NonlinearConstraint(constraint_function, constraint_limit_lower, np.inf)
     pbounds = {'x': (0, 6), 'y': (0, 6)}
 
     optimizer = BayesianOptimization(
         f=target_function,
-        constraint=conmod,
+        constraint=constraint,
         pbounds=pbounds,
         verbose=0,
         random_state=1,
@@ -58,45 +61,55 @@ def constraint_function(x, y):
         n_iter=10,
     )
 
-    assert constraint_function(**optimizer.max["params"]) <= constraint_limit
+    assert constraint_function(**optimizer.max["params"]) >= constraint_limit_lower
 
 
-def test_accurate_approximation_when_known():
+def test_single_constraint_lower_upper():
 
     def target_function(x, y):
         return np.cos(2 * x) * np.cos(y) + np.sin(x)
 
     def constraint_function(x, y):
         return np.cos(x) * np.cos(y) - np.sin(x) * np.sin(y)
 
-    constraint_limit = 0.5
+    constraint_limit_lower = -0.5
+    constraint_limit_upper = 0.5
 
-    conmod = ConstraintModel(constraint_function, constraint_limit)
+    constraint = NonlinearConstraint(constraint_function, constraint_limit_lower, constraint_limit_upper)
     pbounds = {'x': (0, 6), 'y': (0, 6)}
 
     optimizer = BayesianOptimization(
         f=target_function,
-        constraint=conmod,
+        constraint=constraint,
         pbounds=pbounds,
         verbose=0,
         random_state=1,
     )
 
+    assert optimizer.constraint.lb == constraint.lb
+    assert optimizer.constraint.ub == constraint.ub
+
     optimizer.maximize(
         init_points=2,
         n_iter=10,
     )
 
+    # Check limits
+    assert constraint_function(**optimizer.max["params"]) <= constraint_limit_upper
+    assert constraint_function(**optimizer.max["params"]) >= constraint_limit_lower
+
+
     # Exclude the last sampled point, because the constraint is not fitted on that.
     res = np.array([[r['target'], r['constraint'], r['params']['x'], r['params']['y']] for r in optimizer.res[:-1]])
 
     xy = res[:, [2, 3]]
     x = res[:, 2]
     y = res[:, 3]
 
-    assert constraint_function(x, y) == approx(conmod.approx(xy), rel=1e-5, abs=1e-5)
+    # Check accuracy of approximation for sampled points
+    assert constraint_function(x, y) == approx(optimizer.constraint.approx(xy), rel=1e-5, abs=1e-5)
     assert constraint_function(x, y) == approx(optimizer.space.constraint_values[:-1], rel=1e-5, abs=1e-5)
-
+    
 
 def test_multiple_constraints():
 
@@ -109,9 +122,10 @@ def constraint_function_2_dim(x, y):
             -np.cos(x) * np.cos(-y) + np.sin(x) * np.sin(-y)
         ])
 
-    constraint_limit = np.array([0.6, 0.6])
+    constraint_limit_lower = np.array([-np.inf, -np.inf])
+    constraint_limit_upper = np.array([0.6, 0.6])
 
-    conmod = ConstraintModel(constraint_function_2_dim, constraint_limit)
+    conmod = NonlinearConstraint(constraint_function_2_dim, constraint_limit_lower, constraint_limit_upper)
     pbounds = {'x': (0, 6), 'y': (0, 6)}
 
     optimizer = BayesianOptimization(
@@ -127,14 +141,13 @@ def constraint_function_2_dim(x, y):
         n_iter=10,
     )
 
-    assert np.all(
-        constraint_function_2_dim(
-            **optimizer.max["params"]) <= constraint_limit)
+    constraint_at_max = constraint_function_2_dim(**optimizer.max["params"])
+    assert np.all((constraint_at_max <= constraint_limit_upper) & (constraint_at_max >= constraint_limit_lower))
 
     params = optimizer.res[0]["params"]
     x, y = params['x'], params['y']
 
-    assert constraint_function_2_dim(x, y) == approx(conmod.approx(np.array([x, y])), rel=1e-5, abs=1e-5)
+    assert constraint_function_2_dim(x, y) == approx(optimizer.constraint.approx(np.array([x, y])), rel=1e-5, abs=1e-5)
 
 
 def test_kwargs_not_the_same():
@@ -145,14 +158,14 @@ def target_function(x, y):
     def constraint_function(a, b):
         return np.cos(a) * np.cos(b) - np.sin(a) * np.sin(b)
 
-    constraint_limit = 0.5
+    constraint_limit_upper = 0.5
 
-    conmod = ConstraintModel(constraint_function, constraint_limit)
+    constraint = NonlinearConstraint(constraint_function, -np.inf, constraint_limit_upper)
     pbounds = {'x': (0, 6), 'y': (0, 6)}
 
     optimizer = BayesianOptimization(
         f=target_function,
-        constraint=conmod,
+        constraint=constraint,
         pbounds=pbounds,
         verbose=0,
         random_state=1,
@@ -161,4 +174,4 @@ def constraint_function(a, b):
         optimizer.maximize(
             init_points=2,
             n_iter=10,
-        )
+        )