This solves an issue reported in #10 which requires an acquisition to…

… be aware of what points are feasibly or not.

GPflowOpt/acquisition.py

@@ -42,7 +42,7 @@ def __init__(self, models=[], optimize_restarts=5):
         self.models = ParamList(np.atleast_1d(models).tolist())
         self._default_params = list(map(lambda m: m.get_free_state(), self.models))
 
-        assert(optimize_restarts >= 0)
+        assert (optimize_restarts >= 0)
         self._optimize_restarts = optimize_restarts
         self._optimize_all()
 
@@ -69,6 +69,9 @@ def _build_acquisition_wrapper(self, Xcand, gradients=True):
         else:
             return acq
 
+    def build_acquisition(self):
+        raise NotImplementedError
+
     def set_data(self, X, Y):
         num_outputs_sum = 0
         for model in self.models:
@@ -80,7 +83,8 @@ def set_data(self, X, Y):
             model.Y = Ypart
 
         self._optimize_all()
-        self.setup()
+        if self.highest_parent == self:
+            self.setup()
         return num_outputs_sum
 
     @property
@@ -96,6 +100,14 @@ def constraint_indices(self):
     def objective_indices(self):
         return np.setdiff1d(np.arange(self.data[1].shape[1]), self.constraint_indices())
 
+    def feasible_data_index(self):
+        """
+        Returns a boolean array indicating which data points are considered feasible (according to the acquisition 
+        function(s) ) and which not.
+        :return: boolean ndarray, N 
+        """
+        return np.ones(self.data[0].shape[0], dtype=bool)
+
     def setup(self):
         """
         Method triggered after calling set_data(). Override for pre-calculation of quantities used later in
@@ -132,11 +144,10 @@ def __init__(self, model):
 
     def setup(self):
         super(ExpectedImprovement, self).setup()
-        # Obtain the lowest posterior mean for the previous evaluations
-        samples_mean, _ = self.models[0].predict_f(self.data[0])
+        # Obtain the lowest posterior mean for the previous - feasible - evaluations
+        feasible_samples = self.data[0][self.highest_parent.feasible_data_index(), :]
+        samples_mean, _ = self.models[0].predict_f(feasible_samples)
         self.fmin.set_data(np.min(samples_mean, axis=0))
-        # samples_mean, _ = self.models[0].build_predict(self.data[0])
-        # self.fmin = tf.reduce_min(samples_mean, axis=0)
 
     def build_acquisition(self, Xcand):
         # Obtain predictive distributions for candidates
@@ -161,6 +172,10 @@ def __init__(self, model):
     def constraint_indices(self):
         return np.arange(self.data[1].shape[1])
 
+    def feasible_data_index(self):
+        pred = self.evaluate(self.data[0])
+        return pred.ravel() > 0.5
+
     def build_acquisition(self, Xcand):
         candidate_mean, candidate_var = self.models[0].build_predict(Xcand)
         candidate_var = tf.maximum(candidate_var, stability)
@@ -180,6 +195,7 @@ def __init__(self, lhs, rhs, oper):
         self.lhs = lhs
         self.rhs = rhs
         self._oper = oper
+        self.setup()
 
     @Acquisition.data.getter
     def data(self):
@@ -199,10 +215,17 @@ def constraint_indices(self):
         offset = self.lhs.data[1].shape[1]
         return np.hstack((self.lhs.constraint_indices(), offset + self.rhs.constraint_indices()))
 
+    def feasible_data_index(self):
+        return np.all(np.vstack((self.lhs.feasible_data_index(), self.rhs.feasible_data_index())), axis=0)
+
     def build_acquisition(self, Xcand):
         return self._oper(self.lhs.build_acquisition(Xcand), self.rhs.build_acquisition(Xcand),
                           name=self.__class__.__name__)
 
+    def setup(self):
+        self.lhs.setup()
+        self.rhs.setup()
+
 
 class AcquisitionSum(AcquisitionBinaryOperator):
     def __init__(self, lhs, rhs):

GPflowOpt/bo.py

@@ -76,9 +76,7 @@ def _create_result(self, success, message):
         X, Y = self.acquisition.data
 
         # Filter on constraints
-        # Extend with valid column - in case of no constrains
-        Yext = np.hstack((-np.ones((Y.shape[0], 1)), Y[:, self.acquisition.constraint_indices()]))
-        valid = np.all(Yext < 0, axis=1)
+        valid = self.acquisition.feasible_data_index()
 
         if not np.any(valid):
             return OptimizeResult(success=False,

testing/test_acquisition.py

@@ -84,20 +84,23 @@ def test_data_update(self):
         np.testing.assert_allclose(self.acquisition.data[0], X, err_msg="Samples not updated")
         np.testing.assert_allclose(self.acquisition.data[1], Y, err_msg="Values not updated")
 
+    def test_data_indices(self):
+        self.assertTupleEqual(self.acquisition.feasible_data_index().shape, (self.acquisition.data[0].shape[0],),
+                              msg="Incorrect shape returned.")
+
 
 class TestExpectedImprovement(_TestAcquisition, unittest.TestCase):
     def setUp(self):
         super(TestExpectedImprovement, self).setUp()
         self.model = self.create_parabola_model()
-        print(self.model)
         self.acquisition = GPflowOpt.acquisition.ExpectedImprovement(self.model)
 
     def test_object_integrity(self):
         self.assertEqual(len(self.acquisition.models), 1, msg="Model list has incorrect length.")
         self.assertEqual(self.acquisition.models[0], self.model, msg="Incorrect model stored in ExpectedImprovement")
         self.assertEqual(len(self.acquisition._default_params), 1)
-        print(self.acquisition._default_params[0])
-        self.assertTrue(np.allclose(np.sort(self.acquisition._default_params[0]), np.sort(np.array([0.5413, -0.2918, -0.2918, 0.5413])), atol=1e-2),
+        self.assertTrue(np.allclose(np.sort(self.acquisition._default_params[0]),
+                                    np.sort(np.array([0.5413, -0.2918, -0.2918, 0.5413])), atol=1e-2),
                         msg="Initial hypers improperly stored")
         self.assertEqual(self.acquisition.objective_indices(), np.arange(1, dtype=int),
                          msg="ExpectedImprovement returns all objectives")
@@ -122,6 +125,7 @@ def test_EI_validity(self):
         ei1 = self.acquisition.evaluate(Xborder)
         ei2 = self.acquisition.evaluate(Xcenter)
         self.assertGreater(np.min(ei2), np.max(ei1))
+        self.assertTrue(np.all(self.acquisition.feasible_data_index()), msg="EI does never invalidate points")
 
 
 class TestProbabilityOfFeasibility(_TestAcquisition, unittest.TestCase):
@@ -134,7 +138,8 @@ def test_object_integrity(self):
         self.assertEqual(len(self.acquisition.models), 1, msg="Model list has incorrect length.")
         self.assertEqual(self.acquisition.models[0], self.model, msg="Incorrect model stored in PoF")
         self.assertEqual(len(self.acquisition._default_params), 1)
-        self.assertTrue(np.allclose(np.sort(self.acquisition._default_params[0]), np.sort(np.array([0.5413, 0.0277,  0.0277,  0.5413])), atol=1e-2),
+        self.assertTrue(np.allclose(np.sort(self.acquisition._default_params[0]),
+                                    np.sort(np.array([0.5413, 0.0277, 0.0277, 0.5413])), atol=1e-2),
                         msg="Initial hypers improperly stored")
         self.assertEqual(self.acquisition.constraint_indices(), np.arange(1, dtype=int),
                          msg="PoF returns all constraints")
@@ -144,6 +149,8 @@ def test_PoF_validity(self):
         X2 = np.random.rand(10, 2) / 2 + 0.5
         self.assertTrue(np.allclose(self.acquisition.evaluate(X1), 1), msg="Left half of plane is feasible")
         self.assertTrue(np.allclose(self.acquisition.evaluate(X2), 0), msg="Right half of plane is not feasible")
+        self.assertTrue(np.array_equal(self.acquisition.feasible_data_index(), self.acquisition.data[0][:, 0] <= 0.0),
+                        msg='Incorrect feasible indices returned')
 
 
 class _TestAcquisitionBinaryOperator(_TestAcquisition):
@@ -233,13 +240,18 @@ def test_constrained_EI(self):
         design = GPflowOpt.design.FactorialDesign(4, self.domain)
         X = design.generate()
         Yo = parabola2d(X)
-        Yc = plane(X)
+        Yc = -parabola2d(X) + 0.5
         m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0)))
-        m2 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0) / 2))
-        joint = GPflowOpt.acquisition.ExpectedImprovement(m1) * GPflowOpt.acquisition.ProbabilityOfFeasibility(m2)
+        m2 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0)))
+        ei = GPflowOpt.acquisition.ExpectedImprovement(m1)
+        pof = GPflowOpt.acquisition.ProbabilityOfFeasibility(m2)
+        joint = ei * pof
 
         np.testing.assert_allclose(joint.objective_indices(), np.array([0], dtype=int))
         np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int))
+        self.assertGreater(ei.fmin.value, np.min(Yo), msg="The best objective value is in an infeasible area")
+        self.assertTrue(np.allclose(ei.fmin.value, np.min(Yo[pof.feasible_data_index(), :]), atol=1e-3),
+                        msg="fmin computed incorrectly")
 
     def test_hierarchy(self):
         design = GPflowOpt.design.FactorialDesign(4, self.domain)