From da1d1d8bdc3f85d3eb6e6ba2d6ce13fb52bfa504 Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Sun, 13 Aug 2017 04:55:52 +0200
Subject: [PATCH 01/15] First working version of the lazy setup

---
 GPflowOpt/acquisition/acquisition.py |  71 +++--
 GPflowOpt/bo.py                      |   1 +
 testing/test_acquisition.py          | 419 +++++++--------------------
 3 files changed, 153 insertions(+), 338 deletions(-)

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index 907ef8a..0878832 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -22,10 +22,30 @@
 import tensorflow as tf
 
 import copy
+from functools import wraps
 
 float_type = settings.dtypes.float_type
 
 
+class Setup(object):
+    def __call__(self, af_method):
+        @wraps(af_method)
+        def runnable(*args, **kwargs):
+            hp = args[0].highest_parent
+            if hp._needs_setup:
+                # 1 - optimize
+                hp._optimize_models()
+                # 2 - setup
+                # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition
+                # e.g. through feasible_data_index.
+                hp._needs_setup = False
+                hp.setup()
+            results = af_method(*args, **kwargs)
+            return results
+
+        return runnable
+
+
 class Acquisition(Parameterized):
     """
     An acquisition function maps the belief represented by a Bayesian model into a
@@ -53,7 +73,7 @@ def __init__(self, models=[], optimize_restarts=5):
 
         assert (optimize_restarts >= 0)
         self.optimize_restarts = optimize_restarts
-        self._optimize_models()
+        self._needs_setup = True
 
     def _optimize_models(self):
         """
@@ -101,7 +121,7 @@ def enable_scaling(self, domain):
         for m in self.models:
             m.input_transform = domain >> UnitCube(n_inputs)
             m.normalize_output = True
-        self._optimize_models()
+        self.highest_parent._needs_setup = True
 
     def set_data(self, X, Y):
         """
@@ -124,11 +144,7 @@ def set_data(self, X, Y):
             model.X = X
             model.Y = Ypart
 
-        self._optimize_models()
-
-        # Only call setup for the high-level acquisition function
-        if self.highest_parent == self:
-            self.setup()
+        self.highest_parent._needs_setup = True
         return num_outputs_sum
 
     @property
@@ -191,6 +207,7 @@ def setup(self):
         """
         pass
 
+    @Setup()
     @AutoFlow((float_type, [None, None]))
     def evaluate_with_gradients(self, Xcand):
         """
@@ -202,6 +219,7 @@ def evaluate_with_gradients(self, Xcand):
         acq = self.build_acquisition(Xcand)
         return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0]
 
+    @Setup()
     @AutoFlow((float_type, [None, None]))
     def evaluate(self, Xcand):
         """
@@ -237,6 +255,11 @@ def __mul__(self, other):
             return AcquisitionProduct([self] + other.operands.sorted_params)
         return AcquisitionProduct([self, other])
 
+    def __setattr__(self, key, value):
+        super(Acquisition, self).__setattr__(key, value)
+        if key is '_parent':
+            self._needs_setup = True
+
 
 class AcquisitionAggregation(Acquisition):
     """
@@ -252,10 +275,10 @@ def __init__(self, operands, oper):
         assert (all([isinstance(x, Acquisition) for x in operands]))
         self.operands = ParamList(operands)
         self._oper = oper
-        self.setup()
 
     def _optimize_models(self):
-        pass
+        for oper in self.operands:
+            oper._optimize_models()
 
     @Acquisition.models.getter
     def models(self):
@@ -269,13 +292,23 @@ def set_data(self, X, Y):
         offset = 0
         for operand in self.operands:
             offset += operand.set_data(X, Y[:, offset:])
-        if self.highest_parent == self:
-            self.setup()
         return offset
 
-    def setup(self):
+    def _setup_constraints(self):
         for oper in self.operands:
-            oper.setup()
+            if oper.constraint_indices().size > 0:
+                oper._setup_constraints() if isinstance(oper, AcquisitionAggregation) else oper.setup()
+
+    def _setup_objectives(self):
+        for oper in self.operands:
+            if oper.constraint_indices().size == 0:
+                oper._setup_objectives() if isinstance(oper, AcquisitionAggregation) else oper.setup()
+
+    def setup(self):
+        # First setup acquisitions involving constraints
+        self._setup_constraints()
+        # Then objectives
+        self._setup_objectives()
 
     def constraint_indices(self):
         offset = [0]
@@ -300,7 +333,6 @@ class AcquisitionSum(AcquisitionAggregation):
     """
     Sum of acquisition functions
     """
-
     def __init__(self, operands):
         super(AcquisitionSum, self).__init__(operands, tf.reduce_sum)
 
@@ -315,7 +347,6 @@ class AcquisitionProduct(AcquisitionAggregation):
     """
     Product of acquisition functions
     """
-
     def __init__(self, operands):
         super(AcquisitionProduct, self).__init__(operands, tf.reduce_prod)
 
@@ -346,10 +377,11 @@ def __init__(self, acquisition, n_slices, **kwargs):
         # the call to the constructor of the parent classes, will optimize acquisition, so it obtains the MLE solution.
         super(MCMCAcquistion, self).__init__([acquisition] + copies)
         self._sample_opt = kwargs
-        # Update the hyperparameters of the copies using HMC
-        self._update_hyper_draws()
 
-    def _update_hyper_draws(self):
+    def _optimize_models(self):
+        self.operands[0]._optimize_models()
+
+        # Draw samples using HMC
         # Sample each model of the acquisition function - results in a list of 2D ndarrays.
         hypers = np.hstack([model.sample(len(self.operands), **self._sample_opt) for model in self.models])
 
@@ -367,9 +399,6 @@ def set_data(self, X, Y):
             # This triggers model.optimize() on self.operands[0]
             # All copies have optimization disabled, but must have update data.
             offset = operand.set_data(X, Y)
-        self._update_hyper_draws()
-        if self.highest_parent == self:
-            self.setup()
         return offset
 
     def build_acquisition(self, Xcand):
diff --git a/GPflowOpt/bo.py b/GPflowOpt/bo.py
index ec95058..c54f4f9 100644
--- a/GPflowOpt/bo.py
+++ b/GPflowOpt/bo.py
@@ -23,6 +23,7 @@
 from .design import Design, EmptyDesign
 from .pareto import non_dominated_sort
 
+
 class BayesianOptimizer(Optimizer):
     """
     A traditional Bayesian optimization framework implementation.
diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py
index 22151a8..428fbfa 100644
--- a/testing/test_acquisition.py
+++ b/testing/test_acquisition.py
@@ -6,6 +6,7 @@
 import os
 
 
+## CONVENIENT FUNCTIONS ##
 def parabola2d(X):
     return np.atleast_2d(np.sum(X ** 2, axis=1)).T
 
@@ -23,71 +24,80 @@ def vlmop2(x):
     return np.hstack((y1, y2))
 
 
+def load_data(file):
+    path = os.path.dirname(os.path.realpath(__file__))
+    return np.load(os.path.join(path, 'data', file))
 
-class _TestAcquisition(object):
-    """
-    Defines some basic verifications for all acquisition functions. Test classes can derive from this
-    """
 
-    _multiprocess_can_split_ = True
+def create_parabola_model(domain, design=None):
+    if design is None:
+        design = GPflowOpt.design.LatinHyperCube(16, domain)
+    X, Y = design.generate(), parabola2d(design.generate())
+    m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
+    return m
 
-    @property
-    def domain(self):
-        return np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
 
-    def load_data(self, file):
-        path = os.path.dirname(os.path.realpath(__file__))
-        return np.load(os.path.join(path, 'data', file))
+def create_plane_model(domain, design=None):
+    if design is None:
+        design = GPflowOpt.design.LatinHyperCube(25, domain)
+    X, Y = design.generate(), plane(design.generate())
+    m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
+    return m
 
-    def create_parabola_model(self, design=None):
-        if design is None:
-            design = GPflowOpt.design.LatinHyperCube(16, self.domain)
-        X, Y = design.generate(), parabola2d(design.generate())
-        m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
-        return m
 
-    def create_plane_model(self, design=None):
-        if design is None:
-            design = GPflowOpt.design.LatinHyperCube(25, self.domain)
-        X, Y = design.generate(), plane(design.generate())
-        m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
-        return m
+def create_vlmop2_model():
+    data = load_data('vlmop.npz')
+    m1 = GPflow.gpr.GPR(data['X'], data['Y'][:, [0]], kern=GPflow.kernels.Matern32(2))
+    m2 = GPflow.gpr.GPR(data['X'], data['Y'][:, [1]], kern=GPflow.kernels.Matern32(2))
+    return [m1, m2]
 
-    def create_vlmop2_model(self):
-        data = self.load_data('vlmop.npz')
-        m1 = GPflow.gpr.GPR(data['X'], data['Y'][:, [0]], kern=GPflow.kernels.Matern32(2))
-        m2 = GPflow.gpr.GPR(data['X'], data['Y'][:, [1]], kern=GPflow.kernels.Matern32(2))
-        return [m1, m2]
 
-    def setUp(self):
-        self.acquisition = None
-        self.model = None
+## TESTS ##
 
-    def test_result_shape(self):
-        # Verify the returned shape of evaluate
-        design = GPflowOpt.design.RandomDesign(50, self.domain)
+class SimpleAcquisition(GPflowOpt.acquisition.Acquisition):
+    def __init__(self, model):
+        super(SimpleAcquisition, self).__init__(model)
+        self.counter = 0
 
-        with tf.Graph().as_default():
-            free_vars = tf.placeholder(tf.float64, [None])
-            l = self.acquisition.make_tf_array(free_vars)
-            x_tf = tf.placeholder(tf.float64, shape=(50, 2))
-            with self.acquisition.tf_mode():
-                tens = self.acquisition.build_acquisition(x_tf)
-                self.assertTrue(isinstance(tens, tf.Tensor), msg="no Tensor was returned")
-                # tf_shape = tens.get_shape().as_list()
-                # self.assertEqual(tf_shape[0], 50, msg="Tensor of incorrect shape returned")
-                # self.assertTrue(tf_shape[1] == 1 or tf_shape[1] is None)
+    def setup(self):
+        super(SimpleAcquisition, self).setup()
+        self.counter += 1
 
-        res = self.acquisition.evaluate(design.generate())
-        self.assertTupleEqual(res.shape, (50, 1),
-                              msg="Incorrect shape returned for evaluate of {0}".format(self.__class__.__name__))
-        res = self.acquisition.evaluate_with_gradients(design.generate())
-        self.assertTrue(isinstance(res, tuple))
-        self.assertTrue(len(res), 2)
-        self.assertTupleEqual(res[0].shape, (50, 1),
-                              msg="Incorrect shape returned for evaluate of {0}".format(self.__class__.__name__))
-        self.assertTupleEqual(res[1].shape, (50, self.domain.size),
-                              msg="Incorrect shape returned for gradient of {0}".format(self.__class__.__name__))
+    def build_acquisition(self, Xcand):
+        return self.models[0].build_predict(Xcand)[0]
+
+
+class TestAcquisition(unittest.TestCase):
+
+    def setUp(self):
+        self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
+        self.model = create_parabola_model(self.domain)
+        self.acquisition = SimpleAcquisition(self.model)
+
+    def run_setup(self):
+        # Optimize models & perform acquisition setup call.
+        self.acquisition._optimize_models()
+        self.acquisition.setup()
+
+    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.")
+
+    def test_setup_trigger(self):
+        m = create_parabola_model(self.domain)
+        self.assertTrue(np.allclose(m.get_free_state(), self.acquisition.models[0].get_free_state()))
+        self.assertTrue(self.acquisition._needs_setup)
+        self.assertEqual(self.acquisition.counter, 0)
+        self.acquisition.evaluate(GPflowOpt.design.RandomDesign(10, self.domain).generate())
+        self.assertFalse(self.acquisition._needs_setup)
+        self.assertEqual(self.acquisition.counter, 1)
+        self.assertFalse(np.allclose(m.get_free_state(), self.acquisition.models[0].get_free_state()))
+
+        self.acquisition._needs_setup = True
+        self.acquisition.models[0].set_state(m.get_free_state())
+        self.acquisition.evaluate_with_gradients(GPflowOpt.design.RandomDesign(10, self.domain).generate())
+        self.assertFalse(self.acquisition._needs_setup)
+        self.assertEqual(self.acquisition.counter, 2)
 
     def test_data(self):
         # Test the data property
@@ -101,287 +111,58 @@ def test_data(self):
                                 msg="data property should return Tensors")
 
     def test_data_update(self):
-        # Verify a data update
+        # Verify the effect of setting the data
         design = GPflowOpt.design.RandomDesign(10, self.domain)
         X = np.vstack((self.acquisition.data[0], design.generate()))
-        Y = np.hstack([parabola2d(X)] * self.acquisition.data[1].shape[1])
+        Y = parabola2d(X)
+        self.acquisition._needs_setup = False
         self.acquisition.set_data(X, Y)
         np.testing.assert_allclose(self.acquisition.data[0], X, atol=1e-5, err_msg="Samples not updated")
         np.testing.assert_allclose(self.acquisition.data[1], Y, atol=1e-5, err_msg="Values not updated")
+        self.assertTrue(self.acquisition._needs_setup)
 
     def test_data_indices(self):
-        self.assertTupleEqual(self.acquisition.feasible_data_index().shape, (self.acquisition.data[0].shape[0],),
-                              msg="Incorrect shape returned.")
-
-    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)
-        self.assertTrue(
-            np.allclose(np.sort(self.acquisition._default_params[0]), np.sort(np.array([0.5413] * 4)), atol=1e-2),
-            msg="Initial hypers improperly stored")
+        # Return all data as feasible.
+        self.assertTupleEqual(self.acquisition.feasible_data_index().shape, (self.acquisition.data[0].shape[0],))
 
     def test_enable_scaling(self):
         self.assertFalse(
             any(m.wrapped.X.value in GPflowOpt.domain.UnitCube(self.domain.size) for m in self.acquisition.models))
+        self.acquisition._needs_setup = False
         self.acquisition.enable_scaling(self.domain)
-        print(self.acquisition.models[0].wrapped.X.value)
         self.assertTrue(
             all(m.wrapped.X.value in GPflowOpt.domain.UnitCube(self.domain.size) for m in self.acquisition.models))
+        self.assertTrue(self.acquisition._needs_setup)
 
+    def test_result_shape_tf(self):
+        # Verify the returned shape of evaluate
+        design = GPflowOpt.design.RandomDesign(50, self.domain)
 
-class TestExpectedImprovement(_TestAcquisition, unittest.TestCase):
-    def setUp(self):
-        super(TestExpectedImprovement, self).setUp()
-        self.model = self.create_parabola_model()
-        self.acquisition = GPflowOpt.acquisition.ExpectedImprovement(self.model)
-
-    def test_objective_indices(self):
-        self.assertEqual(self.acquisition.objective_indices(), np.arange(1, dtype=int),
-                         msg="ExpectedImprovement returns all objectives")
-
-    def test_setup(self):
-        fmin = np.min(self.acquisition.data[1])
-        self.assertGreater(self.acquisition.fmin.value, 0, msg="The minimum (0) is not amongst the design.")
-        self.assertTrue(np.allclose(self.acquisition.fmin.value, fmin, atol=1e-2), msg="fmin computed incorrectly")
-
-        # Now add the actual minimum
-        p = np.array([[0.0, 0.0]])
-        self.acquisition.set_data(np.vstack((self.acquisition.data[0], p)),
-                                  np.vstack((self.acquisition.data[1], parabola2d(p))))
-        self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated")
-
-    def test_EI_validity(self):
-        Xcenter = np.random.rand(20, 2) * 0.25 - 0.125
-        X = np.random.rand(100, 2) * 2 - 1
-        hor_idx = np.abs(X[:, 0]) > 0.8
-        ver_idx = np.abs(X[:, 1]) > 0.8
-        Xborder = np.vstack((X[hor_idx, :], X[ver_idx, :]))
-        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 TestProbabilityOfImprovement(_TestAcquisition, unittest.TestCase):
-    def setUp(self):
-        super(TestProbabilityOfImprovement, self).setUp()
-        self.model = self.create_parabola_model()
-        self.acquisition = GPflowOpt.acquisition.ProbabilityOfImprovement(self.model)
-
-    def test_objective_indices(self):
-        self.assertEqual(self.acquisition.objective_indices(), np.arange(1, dtype=int),
-                         msg="PoI returns all objectives")
-
-    def test_setup(self):
-        fmin = np.min(self.acquisition.data[1])
-        self.assertGreater(self.acquisition.fmin.value, 0, msg="The minimum (0) is not amongst the design.")
-        self.assertTrue(np.allclose(self.acquisition.fmin.value, fmin, atol=1e-2), msg="fmin computed incorrectly")
-
-        # Now add the actual minimum
-        p = np.array([[0.0, 0.0]])
-        self.acquisition.set_data(np.vstack((self.acquisition.data[0], p)),
-                                  np.vstack((self.acquisition.data[1], parabola2d(p))))
-        self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated")
-
-
-class TestProbabilityOfFeasibility(_TestAcquisition, unittest.TestCase):
-    def setUp(self):
-        super(TestProbabilityOfFeasibility, self).setUp()
-        self.model = self.create_plane_model()
-        self.acquisition = GPflowOpt.acquisition.ProbabilityOfFeasibility(self.model)
-
-    def test_constraint_indices(self):
-        self.assertEqual(self.acquisition.constraint_indices(), np.arange(1, dtype=int),
-                         msg="PoF returns all constraints")
-
-    def test_PoF_validity(self):
-        X1 = np.random.rand(10, 2) / 4
-        X2 = np.random.rand(10, 2) / 4 + 0.75
-        self.assertTrue(np.all(self.acquisition.evaluate(X1) > 0.85), msg="Left half of plane is feasible")
-        self.assertTrue(np.all(self.acquisition.evaluate(X2) < 0.15), msg="Right half of plane is feasible")
-        self.assertTrue(np.all(self.acquisition.evaluate(X1) > self.acquisition.evaluate(X2).T))
-
-
-class TestLowerConfidenceBound(_TestAcquisition, unittest.TestCase):
-    def setUp(self):
-        super(TestLowerConfidenceBound, self).setUp()
-        self.model = self.create_plane_model()
-        self.acquisition = GPflowOpt.acquisition.LowerConfidenceBound(self.model, 3.2)
-
-    def test_objective_indices(self):
-        self.assertEqual(self.acquisition.objective_indices(), np.arange(1, dtype=int),
-                         msg="LCB returns all objectives")
-
-    def test_object_integrity(self):
-        super(TestLowerConfidenceBound, self).test_object_integrity()
-        self.assertEqual(self.acquisition.sigma, 3.2)
-
-    def test_LCB_validity(self):
-        design = GPflowOpt.design.RandomDesign(200, self.domain).generate()
-        p = self.acquisition.models[0].predict_f(design)[0]
-        q = self.acquisition.evaluate(design)
-        np.testing.assert_array_less(q, p)
-
-    def test_LCB_validity_2(self):
-        design = GPflowOpt.design.RandomDesign(200, self.domain).generate()
-        self.acquisition.sigma = 0
-        p = self.acquisition.models[0].predict_f(design)[0]
-        q = self.acquisition.evaluate(design)
-        np.testing.assert_allclose(q, p)
-
-
-class TestHVProbabilityOfImprovement(_TestAcquisition, unittest.TestCase):
-
-    @_TestAcquisition.domain.getter
-    def domain(self):
-        return np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -2, 2) for i in range(1, 3)])
-
-    def setUp(self):
-        super(TestHVProbabilityOfImprovement, self).setUp()
-        self.model = self.create_vlmop2_model()
-        data = self.load_data('vlmop.npz')
-        self.candidates = data['candidates']
-        self.acquisition = GPflowOpt.acquisition.HVProbabilityOfImprovement(self.model)
-
-    def test_object_integrity(self):
-        self.assertEqual(len(self.acquisition.models), 2, msg="Model list has incorrect length.")
-        for m1, m2 in zip(self.acquisition.models, self.model):
-            self.assertEqual(m1, m2, msg="Incorrect model stored in ExpectedImprovement")
-        self.assertEqual(len(self.acquisition._default_params), 2)
-        for i in np.arange(2):
-            self.assertTrue(np.allclose(np.sort(self.acquisition._default_params[i]), np.sort(np.array([0.5413] * 3)),
-                                        atol=1e-2), msg="Initial hypers improperly stored")
-
-    def test_hvpoi_validity(self):
-        scores = self.acquisition.evaluate(self.candidates)
-        np.testing.assert_almost_equal(scores.ravel(), np.array(
-            [2.23723742e-03, 1.00906739e-03, 1.21152340e-02, 6.51774004e-03, 4.42413300e-03, 3.99320061e-02,
-             6.24365778e-04, 1.67279166e-02, 3.70006497e-03, 2.79794264e-02, 1.33966839e-02, 5.08016917e-03,
-             7.85351395e-04, 1.69967446e-02, 5.16896760e-03, 3.87581677e-05, 2.59530418e-03, 1.42613142e-02,
-             4.71508049e-02, 1.01988869e-02, 4.27149696e-04, 2.20649794e-02, 0.00000000e+00]), decimal=2,
-                                       err_msg="hvPoI ranker produced the wrong candidate scores")
-
-
-class _TestAcquisitionAggregation(_TestAcquisition):
-    def test_object_integrity(self):
-        for oper in self.acquisition.operands:
-            self.assertTrue(isinstance(oper, GPflowOpt.acquisition.Acquisition),
-                            msg="All operands should be an acquisition object")
-        self.assertEqual(len(self.acquisition._default_params), 0)
-        self.assertListEqual(self.acquisition.models.sorted_params, self.models)
-
-    def test_data(self):
-        super(_TestAcquisitionAggregation, self).test_data()
-        np.testing.assert_allclose(self.acquisition.data[0], self.acquisition[0].data[0],
-                                   err_msg="Samples should be equal for all operands")
-        np.testing.assert_allclose(self.acquisition.data[0], self.acquisition[1].data[0],
-                                   err_msg="Samples should be equal for all operands")
-
-        Y = np.hstack(map(lambda model: model.Y.value, self.acquisition.models))
-        np.testing.assert_allclose(self.acquisition.data[1], Y,
-                                   err_msg="Value should be horizontally concatenated")
-
-    def test_enable_scaling(self):
-        for oper in self.acquisition.operands:
-            self.assertFalse(any(m.wrapped.X.value in GPflowOpt.domain.UnitCube(self.domain.size) for m in oper.models))
-        self.acquisition.enable_scaling(self.domain)
-        for oper in self.acquisition.operands:
-            self.assertTrue(all(m.wrapped.X.value in GPflowOpt.domain.UnitCube(self.domain.size) for m in oper.models))
-
-
-class TestAcquisitionSum(_TestAcquisitionAggregation, unittest.TestCase):
-    def setUp(self):
-        super(TestAcquisitionSum, self).setUp()
-        self.models = [self.create_parabola_model(), self.create_parabola_model()]
-        self.acquisition = GPflowOpt.acquisition.AcquisitionSum([
-            GPflowOpt.acquisition.ExpectedImprovement(self.models[0]),
-            GPflowOpt.acquisition.ExpectedImprovement(self.models[1])
-        ])
-
-    def test_sum_validity(self):
-        design = GPflowOpt.design.FactorialDesign(4, self.domain)
-        m = self.create_parabola_model()
-        single_ei = GPflowOpt.acquisition.ExpectedImprovement(m)
-        p1 = self.acquisition.evaluate(design.generate())
-        p2 = single_ei.evaluate(design.generate())
-        np.testing.assert_allclose(p2, p1 / 2, rtol=1e-3, err_msg="The sum of 2 EI should be the double of only EI")
-
-    def test_generating_operator(self):
-        joint = GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model()) + \
-                GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model())
-        self.assertTrue(isinstance(joint, GPflowOpt.acquisition.AcquisitionSum))
-
-    def test_indices(self):
-        np.testing.assert_allclose(self.acquisition.objective_indices(), np.arange(2, dtype=int),
-                                   err_msg="Sum of two EI should return all objectives")
-        np.testing.assert_allclose(self.acquisition.constraint_indices(), np.arange(0, dtype=int),
-                                   err_msg="Sum of two EI should return no constraints")
-
-
-class TestAcquisitionProduct(_TestAcquisitionAggregation, unittest.TestCase):
-    def setUp(self):
-        super(TestAcquisitionProduct, self).setUp()
-        self.models = [self.create_parabola_model(), self.create_parabola_model()]
-        self.acquisition = GPflowOpt.acquisition.AcquisitionProduct([
-            GPflowOpt.acquisition.ExpectedImprovement(self.models[0]),
-            GPflowOpt.acquisition.ExpectedImprovement(self.models[1])
-        ])
-
-    def test_product_validity(self):
-        design = GPflowOpt.design.FactorialDesign(4, self.domain)
-        m = self.create_parabola_model()
-        single_ei = GPflowOpt.acquisition.ExpectedImprovement(m)
-        p1 = self.acquisition.evaluate(design.generate())
-        p2 = single_ei.evaluate(design.generate())
-        np.testing.assert_allclose(p2, np.sqrt(p1), rtol=1e-3,
-                                   err_msg="The product of 2 EI should be the square of one EI")
-
-    def test_generating_operator(self):
-        joint = GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model()) * \
-                GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model())
-        self.assertTrue(isinstance(joint, GPflowOpt.acquisition.AcquisitionProduct))
-
-    def test_indices(self):
-        np.testing.assert_allclose(self.acquisition.objective_indices(), np.arange(2, dtype=int),
-                                   err_msg="Product of two EI should return all objectives")
-        np.testing.assert_allclose(self.acquisition.constraint_indices(), np.arange(0, dtype=int),
-                                   err_msg="Product of two EI should return no constraints")
-
-
-class TestMCMCAcquisition(_TestAcquisitionAggregation, unittest.TestCase):
-    def setUp(self):
-        super(TestMCMCAcquisition, self).setUp()
-        self.models = [self.create_parabola_model()]
-        self.acquisition = GPflowOpt.acquisition.MCMCAcquistion(GPflowOpt.acquisition.ExpectedImprovement(self.models[0]), 5)
-
-    def test_hyper_updates(self):
-        orig_hypers = [c.get_free_state() for c in self.acquisition.operands[1:]]
-        self.acquisition._update_hyper_draws()
-        for co, cn in zip(orig_hypers, [c.get_free_state() for c in self.acquisition.operands[1:]]):
-            self.assertFalse(np.allclose(co, cn))
+        with tf.Graph().as_default():
+            free_vars = tf.placeholder(tf.float64, [None])
+            l = self.acquisition.make_tf_array(free_vars)
+            x_tf = tf.placeholder(tf.float64, shape=(50, 2))
+            with self.acquisition.tf_mode():
+                tens = self.acquisition.build_acquisition(x_tf)
+                self.assertTrue(isinstance(tens, tf.Tensor), msg="no Tensor was returned")
 
-    def test_marginalized_score(self):
-        Xt = np.random.rand(20, 2) * 2 - 1
-        ei_mle = self.acquisition.operands[0].evaluate(Xt)
-        ei_mcmc = self.acquisition.evaluate(Xt)
-        np.testing.assert_almost_equal(ei_mle, ei_mcmc, decimal=5)
+    def test_result_shape_np(self):
+        design = GPflowOpt.design.RandomDesign(50, self.domain)
+        res = self.acquisition.evaluate(design.generate())
+        self.assertTupleEqual(res.shape, (50, 1))
+        res = self.acquisition.evaluate_with_gradients(design.generate())
+        self.assertTrue(isinstance(res, tuple))
+        self.assertTrue(len(res), 2)
+        self.assertTupleEqual(res[0].shape, (50, 1))
+        self.assertTupleEqual(res[1].shape, (50, self.domain.size))
 
 
 class TestJointAcquisition(unittest.TestCase):
 
     _multiprocessing_can_split_ = True
 
-    @property
-    def domain(self):
-        return np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
-
-    def create_parabola_model(self, design=None):
-        if design is None:
-            design = GPflowOpt.design.LatinHyperCube(16, self.domain)
-        X, Y = design.generate(), parabola2d(design.generate())
-        m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
-        return m
+    def setUp(self):
+        self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
 
     def test_constrained_EI(self):
         design = GPflowOpt.design.LatinHyperCube(16, self.domain)
@@ -394,8 +175,13 @@ def test_constrained_EI(self):
         pof = GPflowOpt.acquisition.ProbabilityOfFeasibility(m2)
         joint = ei * pof
 
+        # Test output indices
         np.testing.assert_allclose(joint.objective_indices(), np.array([0], dtype=int))
         np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int))
+
+        # Test proper setup
+        joint._optimize_models()
+        joint.setup()
         self.assertGreater(ei.fmin.value, np.min(ei.data[1]), msg="The best objective value is in an infeasible area")
         self.assertTrue(np.allclose(ei.fmin.value, np.min(ei.data[1][pof.feasible_data_index(), :]), atol=1e-3),
                         msg="fmin computed incorrectly")
@@ -404,8 +190,8 @@ def test_hierarchy(self):
         design = GPflowOpt.design.LatinHyperCube(16, self.domain)
         X = design.generate()
         Yc = plane(X)
-        m1 = self.create_parabola_model()
-        m2 = self.create_parabola_model()
+        m1 = create_parabola_model(self.domain, design)
+        m2 = create_parabola_model(self.domain, design)
         m3 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True))
         joint = GPflowOpt.acquisition.ExpectedImprovement(m1) * \
                 (GPflowOpt.acquisition.ProbabilityOfFeasibility(m3)
@@ -415,8 +201,8 @@ def test_hierarchy(self):
         np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int))
 
     def test_multi_aggr(self):
-        models = [self.create_parabola_model(), self.create_parabola_model(), self.create_parabola_model()]
-        acq1, acq2, acq3 = tuple(map(lambda m: GPflowOpt.acquisition.ExpectedImprovement(m), models))
+        acq = [GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(self.domain)) for i in range(4)]
+        acq1, acq2, acq3, acq4 = acq
         joint = acq1 + acq2 + acq3
         self.assertIsInstance(joint, GPflowOpt.acquisition.AcquisitionSum)
         self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3])
@@ -439,7 +225,6 @@ def test_multi_aggr(self):
         self.assertIsInstance(joint, GPflowOpt.acquisition.AcquisitionProduct)
         self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3])
 
-        acq4 = GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model())
         first = acq1 + acq2
         second = acq3 + acq4
         joint = first + second

From b78db1575ead93124b1e129ffc3cca330172696e Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Mon, 14 Aug 2017 11:14:48 +0200
Subject: [PATCH 02/15] Refactoring test files, tests became a bit tangled

---
 GPflowOpt/acquisition/acquisition.py |   1 +
 testing/data/vlmop.npz               | Bin 1834 -> 2194 bytes
 testing/test_acquisition.py          | 184 +++++++++++++++++----------
 testing/test_implementations.py      | 133 +++++++++++++++++++
 testing/utility.py                   |  49 +++++++
 5 files changed, 299 insertions(+), 68 deletions(-)
 create mode 100644 testing/test_implementations.py
 create mode 100644 testing/utility.py

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index 0878832..56b4471 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -379,6 +379,7 @@ def __init__(self, acquisition, n_slices, **kwargs):
         self._sample_opt = kwargs
 
     def _optimize_models(self):
+        # Optimize model #1
         self.operands[0]._optimize_models()
 
         # Draw samples using HMC
diff --git a/testing/data/vlmop.npz b/testing/data/vlmop.npz
index e92eb85cad2f03735a3227493b4e144854866948..7e360cf3dfbbb2e34ce3816f62059b35dab70c58 100644
GIT binary patch
delta 470
zcmZ3*H%X8$z?+#xgaHJ8hw*uD<g;gz0x`;D7d+x%WMBYcE}&#_a(+>2@kDv2dSe|0
zLroopS_K&3S{b%x<K^Ha`v>n7?u2#(*y~?x)VqJa!~XC^p4Ny>RrW2mM>V+Qi|l*$
zwVn80xWvA7Z;?-bsf)e->K9e=iWBV}LLD7v`sdm6nWeTaUN*;Gv*m2%j{I(W(f#7>
zAJ3K88!eA;XWZjye_#EC#f-@GiT2uu_UQ2Gl-fHiDQ(~gVzJk|R1mb^Db@aj(k25=
z(_Z^{3EnC-l@<1{6WVjtAGg?tPU!ua_|Mut`@rJO#xthbgS`tRQ2oBymc^Wj2RSr0
z2Qoci<YQzKVMc_|<R&(K8*Cv2Gqj=6nF$tNdU*wv0p4f^d}Nb~E@oz6h`=UWrBy!X
t05Hr!m=9(okZx$)#>&8uoS2uAnUYwNS`5+4$_7%x0))4K^dfc;4*-lYhTH%E

delta 158
zcmbOvxQdT2z?+#xgaHH|Z+YXjk<Xrq1;p5Vg!KU<53<zeMrLy+Ax0(<W`uI5X}{|x
zD{_cUKEo~w);2kSLrTe+iGd+GF)t-EC9x#6STC=jGQb;E-(*=1DS3Ss28IZTEGrwx
N0n9+S2uK%!c>qH}DOmsj

diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py
index 428fbfa..a7e8b0a 100644
--- a/testing/test_acquisition.py
+++ b/testing/test_acquisition.py
@@ -3,56 +3,11 @@
 import numpy as np
 import GPflow
 import tensorflow as tf
-import os
+from parameterized import parameterized
+from .utility import create_parabola_model, parabola2d, plane
 
+domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
 
-## CONVENIENT FUNCTIONS ##
-def parabola2d(X):
-    return np.atleast_2d(np.sum(X ** 2, axis=1)).T
-
-
-def plane(X):
-    return X[:, [0]] - 0.5
-
-
-def vlmop2(x):
-    transl = 1 / np.sqrt(2)
-    part1 = (x[:, [0]] - transl) ** 2 + (x[:, [1]] - transl) ** 2
-    part2 = (x[:, [0]] + transl) ** 2 + (x[:, [1]] + transl) ** 2
-    y1 = 1 - np.exp(-1 * part1)
-    y2 = 1 - np.exp(-1 * part2)
-    return np.hstack((y1, y2))
-
-
-def load_data(file):
-    path = os.path.dirname(os.path.realpath(__file__))
-    return np.load(os.path.join(path, 'data', file))
-
-
-def create_parabola_model(domain, design=None):
-    if design is None:
-        design = GPflowOpt.design.LatinHyperCube(16, domain)
-    X, Y = design.generate(), parabola2d(design.generate())
-    m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
-    return m
-
-
-def create_plane_model(domain, design=None):
-    if design is None:
-        design = GPflowOpt.design.LatinHyperCube(25, domain)
-    X, Y = design.generate(), plane(design.generate())
-    m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
-    return m
-
-
-def create_vlmop2_model():
-    data = load_data('vlmop.npz')
-    m1 = GPflow.gpr.GPR(data['X'], data['Y'][:, [0]], kern=GPflow.kernels.Matern32(2))
-    m2 = GPflow.gpr.GPR(data['X'], data['Y'][:, [1]], kern=GPflow.kernels.Matern32(2))
-    return [m1, m2]
-
-
-## TESTS ##
 
 class SimpleAcquisition(GPflowOpt.acquisition.Acquisition):
     def __init__(self, model):
@@ -70,8 +25,7 @@ def build_acquisition(self, Xcand):
 class TestAcquisition(unittest.TestCase):
 
     def setUp(self):
-        self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
-        self.model = create_parabola_model(self.domain)
+        self.model = create_parabola_model(domain)
         self.acquisition = SimpleAcquisition(self.model)
 
     def run_setup(self):
@@ -84,18 +38,18 @@ def test_object_integrity(self):
         self.assertEqual(self.acquisition.models[0], self.model, msg="Incorrect model stored.")
 
     def test_setup_trigger(self):
-        m = create_parabola_model(self.domain)
+        m = create_parabola_model(domain)
         self.assertTrue(np.allclose(m.get_free_state(), self.acquisition.models[0].get_free_state()))
         self.assertTrue(self.acquisition._needs_setup)
         self.assertEqual(self.acquisition.counter, 0)
-        self.acquisition.evaluate(GPflowOpt.design.RandomDesign(10, self.domain).generate())
+        self.acquisition.evaluate(GPflowOpt.design.RandomDesign(10, domain).generate())
         self.assertFalse(self.acquisition._needs_setup)
         self.assertEqual(self.acquisition.counter, 1)
         self.assertFalse(np.allclose(m.get_free_state(), self.acquisition.models[0].get_free_state()))
 
         self.acquisition._needs_setup = True
         self.acquisition.models[0].set_state(m.get_free_state())
-        self.acquisition.evaluate_with_gradients(GPflowOpt.design.RandomDesign(10, self.domain).generate())
+        self.acquisition.evaluate_with_gradients(GPflowOpt.design.RandomDesign(10, domain).generate())
         self.assertFalse(self.acquisition._needs_setup)
         self.assertEqual(self.acquisition.counter, 2)
 
@@ -112,7 +66,7 @@ def test_data(self):
 
     def test_data_update(self):
         # Verify the effect of setting the data
-        design = GPflowOpt.design.RandomDesign(10, self.domain)
+        design = GPflowOpt.design.RandomDesign(10, domain)
         X = np.vstack((self.acquisition.data[0], design.generate()))
         Y = parabola2d(X)
         self.acquisition._needs_setup = False
@@ -127,16 +81,16 @@ def test_data_indices(self):
 
     def test_enable_scaling(self):
         self.assertFalse(
-            any(m.wrapped.X.value in GPflowOpt.domain.UnitCube(self.domain.size) for m in self.acquisition.models))
+            any(m.wrapped.X.value in GPflowOpt.domain.UnitCube(domain.size) for m in self.acquisition.models))
         self.acquisition._needs_setup = False
-        self.acquisition.enable_scaling(self.domain)
+        self.acquisition.enable_scaling(domain)
         self.assertTrue(
-            all(m.wrapped.X.value in GPflowOpt.domain.UnitCube(self.domain.size) for m in self.acquisition.models))
+            all(m.wrapped.X.value in GPflowOpt.domain.UnitCube(domain.size) for m in self.acquisition.models))
         self.assertTrue(self.acquisition._needs_setup)
 
     def test_result_shape_tf(self):
         # Verify the returned shape of evaluate
-        design = GPflowOpt.design.RandomDesign(50, self.domain)
+        design = GPflowOpt.design.RandomDesign(50, domain)
 
         with tf.Graph().as_default():
             free_vars = tf.placeholder(tf.float64, [None])
@@ -147,25 +101,119 @@ def test_result_shape_tf(self):
                 self.assertTrue(isinstance(tens, tf.Tensor), msg="no Tensor was returned")
 
     def test_result_shape_np(self):
-        design = GPflowOpt.design.RandomDesign(50, self.domain)
+        design = GPflowOpt.design.RandomDesign(50, domain)
         res = self.acquisition.evaluate(design.generate())
         self.assertTupleEqual(res.shape, (50, 1))
         res = self.acquisition.evaluate_with_gradients(design.generate())
         self.assertTrue(isinstance(res, tuple))
         self.assertTrue(len(res), 2)
         self.assertTupleEqual(res[0].shape, (50, 1))
-        self.assertTupleEqual(res[1].shape, (50, self.domain.size))
+        self.assertTupleEqual(res[1].shape, (50, domain.size))
+
+
+aggregations = list()
+aggregations.append(GPflowOpt.acquisition.AcquisitionSum([
+            GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain)),
+            GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain))
+        ]))
+aggregations.append(GPflowOpt.acquisition.AcquisitionProduct([
+            GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain)),
+            GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain))
+        ]))
+aggregations.append(GPflowOpt.acquisition.MCMCAcquistion(
+    GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain)), 5)
+)
+
+
+class TestAcquisitionAggregation(unittest.TestCase):
+    @parameterized.expand(list(zip(aggregations)))
+    def test_object_integrity(self, acquisition):
+        for oper in acquisition.operands:
+            self.assertTrue(isinstance(oper, GPflowOpt.acquisition.Acquisition),
+                            msg="All operands should be an acquisition object")
+        self.assertEqual(len(acquisition._default_params), 0)
+        self.assertTrue(all(isinstance(m, GPflow.model.GPModel) for m in acquisition.models.sorted_params))
+
+    @parameterized.expand(list(zip(aggregations)))
+    def test_data(self, acquisition):
+        np.testing.assert_allclose(acquisition.data[0], acquisition[0].data[0],
+                                   err_msg="Samples should be equal for all operands")
+        np.testing.assert_allclose(acquisition.data[0], acquisition[1].data[0],
+                                   err_msg="Samples should be equal for all operands")
+
+        Y = np.hstack(map(lambda model: model.Y.value, acquisition.models))
+        np.testing.assert_allclose(acquisition.data[1], Y, err_msg="Value should be horizontally concatenated")
+
+    @parameterized.expand(list(zip(aggregations)))
+    def test_enable_scaling(self, acquisition):
+        for oper in acquisition.operands:
+            self.assertFalse(any(m.wrapped.X.value in GPflowOpt.domain.UnitCube(2) for m in oper.models))
+        acquisition.enable_scaling(domain)
+        for oper in acquisition.operands:
+            self.assertTrue(all(m.wrapped.X.value in GPflowOpt.domain.UnitCube(2) for m in oper.models))
+
+    @parameterized.expand(list(zip([aggregations[0]])))
+    def test_sum_validity(self, acquisition):
+        design = GPflowOpt.design.FactorialDesign(4, domain)
+        m = create_parabola_model(domain)
+        single_ei = GPflowOpt.acquisition.ExpectedImprovement(m)
+        p1 = acquisition.evaluate(design.generate())
+        p2 = single_ei.evaluate(design.generate())
+        np.testing.assert_allclose(p2, p1 / 2, rtol=1e-3)
+
+    @parameterized.expand(list(zip([aggregations[1]])))
+    def test_product_validity(self, acquisition):
+        design = GPflowOpt.design.FactorialDesign(4, domain)
+        m = create_parabola_model(domain)
+        single_ei = GPflowOpt.acquisition.ExpectedImprovement(m)
+        p1 = acquisition.evaluate(design.generate())
+        p2 = single_ei.evaluate(design.generate())
+        np.testing.assert_allclose(p2, np.sqrt(p1), rtol=1e-3)
+
+    @parameterized.expand(list(zip(aggregations[0:2])))
+    def test_indices(self, acquisition):
+        np.testing.assert_allclose(acquisition.objective_indices(), np.arange(2, dtype=int))
+        np.testing.assert_allclose(acquisition.constraint_indices(), np.arange(0, dtype=int))
+
+    def test_generating_operators(self):
+        joint = GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain)) + \
+                GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain))
+        self.assertTrue(isinstance(joint, GPflowOpt.acquisition.AcquisitionSum))
+
+        joint = GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain)) * \
+                GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain))
+        self.assertTrue(isinstance(joint, GPflowOpt.acquisition.AcquisitionProduct))
+
+    @parameterized.expand(list(zip([aggregations[2]])))
+    def test_hyper_updates(self, acquisition):
+        orig_hypers = [c.get_free_state() for c in acquisition.operands[1:]]
+        lik_start = acquisition.operands[0].models[0].compute_log_likelihood()
+        acquisition._optimize_models()
+        self.assertGreater(acquisition.operands[0].models[0].compute_log_likelihood(), lik_start)
+
+        for co, cn in zip(orig_hypers, [c.get_free_state() for c in acquisition.operands[1:]]):
+            self.assertFalse(np.allclose(co, cn))
+
+    @parameterized.expand(list(zip([aggregations[2]])))
+    def test_marginalized_score(self, acquisition):
+        acquisition._optimize_models()
+        acquisition.setup()
+        Xt = np.random.rand(20, 2) * 2 - 1
+        ei_mle = acquisition.operands[0].evaluate(Xt)
+        ei_mcmc = acquisition.evaluate(Xt)
+        np.testing.assert_almost_equal(ei_mle, ei_mcmc, decimal=5)
+
+    @parameterized.expand(list(zip([aggregations[2]])))
+    def test_mcmc_acq_models(self, acquisition):
+        self.assertListEqual(acquisition.models.sorted_params, acquisition.operands[0].models.sorted_params)
 
 
 class TestJointAcquisition(unittest.TestCase):
 
     _multiprocessing_can_split_ = True
 
-    def setUp(self):
-        self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
-
     def test_constrained_EI(self):
-        design = GPflowOpt.design.LatinHyperCube(16, self.domain)
+        design = GPflowOpt.design.LatinHyperCube(16, domain)
         X = design.generate()
         Yo = parabola2d(X)
         Yc = -parabola2d(X) + 0.5
@@ -187,11 +235,11 @@ def test_constrained_EI(self):
                         msg="fmin computed incorrectly")
 
     def test_hierarchy(self):
-        design = GPflowOpt.design.LatinHyperCube(16, self.domain)
+        design = GPflowOpt.design.LatinHyperCube(16, domain)
         X = design.generate()
         Yc = plane(X)
-        m1 = create_parabola_model(self.domain, design)
-        m2 = create_parabola_model(self.domain, design)
+        m1 = create_parabola_model(domain, design)
+        m2 = create_parabola_model(domain, design)
         m3 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True))
         joint = GPflowOpt.acquisition.ExpectedImprovement(m1) * \
                 (GPflowOpt.acquisition.ProbabilityOfFeasibility(m3)
@@ -201,7 +249,7 @@ def test_hierarchy(self):
         np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int))
 
     def test_multi_aggr(self):
-        acq = [GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(self.domain)) for i in range(4)]
+        acq = [GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain)) for i in range(4)]
         acq1, acq2, acq3, acq4 = acq
         joint = acq1 + acq2 + acq3
         self.assertIsInstance(joint, GPflowOpt.acquisition.AcquisitionSum)
diff --git a/testing/test_implementations.py b/testing/test_implementations.py
new file mode 100644
index 0000000..3715852
--- /dev/null
+++ b/testing/test_implementations.py
@@ -0,0 +1,133 @@
+import unittest
+import GPflowOpt
+import numpy as np
+from .utility import create_parabola_model, create_plane_model, create_vlmop2_model, parabola2d, load_data
+
+
+class TestExpectedImprovement(unittest.TestCase):
+    def setUp(self):
+        self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
+        self.model = create_parabola_model(self.domain)
+        self.acquisition = GPflowOpt.acquisition.ExpectedImprovement(self.model)
+
+    def test_objective_indices(self):
+        self.assertEqual(self.acquisition.objective_indices(), np.arange(1, dtype=int),
+                         msg="ExpectedImprovement returns all objectives")
+
+    def test_setup(self):
+        self.acquisition._optimize_models()
+        self.acquisition.setup()
+        fmin = np.min(self.acquisition.data[1])
+        self.assertGreater(self.acquisition.fmin.value, 0, msg="The minimum (0) is not amongst the design.")
+        self.assertTrue(np.allclose(self.acquisition.fmin.value, fmin, atol=1e-2), msg="fmin computed incorrectly")
+
+        # Now add the actual minimum
+        p = np.array([[0.0, 0.0]])
+        self.acquisition.set_data(np.vstack((self.acquisition.data[0], p)),
+                                  np.vstack((self.acquisition.data[1], parabola2d(p))))
+        self.acquisition._optimize_models()
+        self.acquisition.setup()
+        self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated")
+
+    def test_EI_validity(self):
+        Xcenter = np.random.rand(20, 2) * 0.25 - 0.125
+        X = np.random.rand(100, 2) * 2 - 1
+        hor_idx = np.abs(X[:, 0]) > 0.8
+        ver_idx = np.abs(X[:, 1]) > 0.8
+        Xborder = np.vstack((X[hor_idx, :], X[ver_idx, :]))
+        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 TestProbabilityOfImprovement(unittest.TestCase):
+    def setUp(self):
+        self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
+        self.model = create_parabola_model(self.domain)
+        self.acquisition = GPflowOpt.acquisition.ProbabilityOfImprovement(self.model)
+
+    def test_objective_indices(self):
+        self.assertEqual(self.acquisition.objective_indices(), np.arange(1, dtype=int),
+                         msg="PoI returns all objectives")
+
+    def test_setup(self):
+        self.acquisition._optimize_models()
+        self.acquisition.setup()
+        fmin = np.min(self.acquisition.data[1])
+        self.assertGreater(self.acquisition.fmin.value, 0, msg="The minimum (0) is not amongst the design.")
+        self.assertTrue(np.allclose(self.acquisition.fmin.value, fmin, atol=1e-2), msg="fmin computed incorrectly")
+
+        # Now add the actual minimum
+        p = np.array([[0.0, 0.0]])
+        self.acquisition.set_data(np.vstack((self.acquisition.data[0], p)),
+                                  np.vstack((self.acquisition.data[1], parabola2d(p))))
+        self.acquisition._optimize_models()
+        self.acquisition.setup()
+        self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated")
+
+
+class TestProbabilityOfFeasibility(unittest.TestCase):
+    def setUp(self):
+        self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
+        self.model = create_plane_model(self.domain)
+        self.acquisition = GPflowOpt.acquisition.ProbabilityOfFeasibility(self.model)
+
+    def test_constraint_indices(self):
+        self.assertEqual(self.acquisition.constraint_indices(), np.arange(1, dtype=int),
+                         msg="PoF returns all constraints")
+
+    def test_PoF_validity(self):
+        X1 = np.random.rand(10, 2) / 4
+        X2 = np.random.rand(10, 2) / 4 + 0.75
+        self.assertTrue(np.all(self.acquisition.evaluate(X1) > 0.85), msg="Left half of plane is feasible")
+        self.assertTrue(np.all(self.acquisition.evaluate(X2) < 0.15), msg="Right half of plane is feasible")
+        self.assertTrue(np.all(self.acquisition.evaluate(X1) > self.acquisition.evaluate(X2).T))
+
+
+class TestLowerConfidenceBound(unittest.TestCase):
+    def setUp(self):
+        self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
+        self.model = create_plane_model(self.domain)
+        self.acquisition = GPflowOpt.acquisition.LowerConfidenceBound(self.model, 3.2)
+
+    def test_objective_indices(self):
+        self.assertEqual(self.acquisition.objective_indices(), np.arange(1, dtype=int),
+                         msg="LCB returns all objectives")
+
+    def test_object_integrity(self):
+        self.assertEqual(self.acquisition.sigma, 3.2)
+
+    def test_LCB_validity(self):
+        design = GPflowOpt.design.RandomDesign(200, self.domain).generate()
+        q = self.acquisition.evaluate(design)
+        p = self.acquisition.models[0].predict_f(design)[0]
+        np.testing.assert_array_less(q, p)
+
+    def test_LCB_validity_2(self):
+        design = GPflowOpt.design.RandomDesign(200, self.domain).generate()
+        self.acquisition.sigma = 0
+        q = self.acquisition.evaluate(design)
+        p = self.acquisition.models[0].predict_f(design)[0]
+        np.testing.assert_allclose(q, p)
+
+
+class TestHVProbabilityOfImprovement(unittest.TestCase):
+
+    def setUp(self):
+        self.model = create_vlmop2_model()
+        self.data = load_data('vlmop.npz')
+        self.acquisition = GPflowOpt.acquisition.HVProbabilityOfImprovement(self.model)
+
+    def test_object_integrity(self):
+        self.assertEqual(len(self.acquisition.models), 2, msg="Model list has incorrect length.")
+        for m1, m2 in zip(self.acquisition.models, self.model):
+            self.assertEqual(m1, m2, msg="Incorrect model stored in ExpectedImprovement")
+        self.assertEqual(len(self.acquisition._default_params), 2)
+        for i in np.arange(2):
+            self.assertTrue(np.allclose(np.sort(self.acquisition._default_params[i]), np.sort(np.array([0.5413] * 3)),
+                                        atol=1e-2), msg="Initial hypers improperly stored")
+
+    def test_HvPoI_validity(self):
+        scores = self.acquisition.evaluate(self.data['candidates'])
+        np.testing.assert_almost_equal(scores, self.data['scores'], decimal=2)
diff --git a/testing/utility.py b/testing/utility.py
new file mode 100644
index 0000000..986aedb
--- /dev/null
+++ b/testing/utility.py
@@ -0,0 +1,49 @@
+import numpy as np
+import GPflow
+import GPflowOpt
+import os
+
+
+def parabola2d(X):
+    return np.atleast_2d(np.sum(X ** 2, axis=1)).T
+
+
+def plane(X):
+    return X[:, [0]] - 0.5
+
+
+def vlmop2(x):
+    transl = 1 / np.sqrt(2)
+    part1 = (x[:, [0]] - transl) ** 2 + (x[:, [1]] - transl) ** 2
+    part2 = (x[:, [0]] + transl) ** 2 + (x[:, [1]] + transl) ** 2
+    y1 = 1 - np.exp(-1 * part1)
+    y2 = 1 - np.exp(-1 * part2)
+    return np.hstack((y1, y2))
+
+
+def load_data(file):
+    path = os.path.dirname(os.path.realpath(__file__))
+    return np.load(os.path.join(path, 'data', file))
+
+
+def create_parabola_model(domain, design=None):
+    if design is None:
+        design = GPflowOpt.design.LatinHyperCube(16, domain)
+    X, Y = design.generate(), parabola2d(design.generate())
+    m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
+    return m
+
+
+def create_plane_model(domain, design=None):
+    if design is None:
+        design = GPflowOpt.design.LatinHyperCube(25, domain)
+    X, Y = design.generate(), plane(design.generate())
+    m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
+    return m
+
+
+def create_vlmop2_model():
+    data = load_data('vlmop.npz')
+    m1 = GPflow.gpr.GPR(data['X'], data['Y'][:, [0]], kern=GPflow.kernels.Matern32(2))
+    m2 = GPflow.gpr.GPR(data['X'], data['Y'][:, [1]], kern=GPflow.kernels.Matern32(2))
+    return [m1, m2]
\ No newline at end of file

From 06499512a889d3dfe38e4a7e705ff09c7acd759d Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Mon, 14 Aug 2017 13:24:43 +0200
Subject: [PATCH 03/15] Account for feasible samples in PoI and HvPoI

---
 GPflowOpt/acquisition/hvpoi.py | 12 ++++++++----
 GPflowOpt/acquisition/poi.py   |  3 ++-
 2 files changed, 10 insertions(+), 5 deletions(-)

diff --git a/GPflowOpt/acquisition/hvpoi.py b/GPflowOpt/acquisition/hvpoi.py
index da40acc..1c310ea 100644
--- a/GPflowOpt/acquisition/hvpoi.py
+++ b/GPflowOpt/acquisition/hvpoi.py
@@ -69,9 +69,12 @@ def __init__(self, models):
         :param models: A list of (possibly multioutput) GPflow representing our belief of the objectives.
         """
         super(HVProbabilityOfImprovement, self).__init__(models)
-        assert self.data[1].shape[1] > 1
-        self.pareto = Pareto(np.hstack((m.predict_f(self.data[0])[0] for m in self.models)))
-        self.reference = DataHolder(self._estimate_reference())
+        num_objectives = self.data[1].shape[1]
+        assert num_objectives > 1
+
+        # Keep pareto empty for now - its updated in setup()
+        self.pareto = Pareto(np.empty((0, num_objectives)))
+        self.reference = DataHolder(np.ones((1, num_objectives)))
 
     def _estimate_reference(self):
         pf = self.pareto.front.value
@@ -85,7 +88,8 @@ def setup(self):
         super(HVProbabilityOfImprovement, self).setup()
 
         # Obtain hypervolume cell bounds, use prediction mean
-        F = np.hstack((m.predict_f(self.data[0])[0] for m in self.models))
+        feasible_samples = self.data[0][self.highest_parent.feasible_data_index(), :]
+        F = np.hstack((m.predict_f(feasible_samples)[0] for m in self.models))
         self.pareto.update(F)
 
         # Calculate reference point.
diff --git a/GPflowOpt/acquisition/poi.py b/GPflowOpt/acquisition/poi.py
index cef81ed..707a59e 100644
--- a/GPflowOpt/acquisition/poi.py
+++ b/GPflowOpt/acquisition/poi.py
@@ -41,7 +41,8 @@ def __init__(self, model):
 
     def setup(self):
         super(ProbabilityOfImprovement, self).setup()
-        samples_mean, _ = self.models[0].predict_f(self.data[0])
+        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))
 
     def build_acquisition(self, Xcand):

From 605f56941be5a71fd19f882afeb2371aa2077a75 Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Tue, 15 Aug 2017 00:39:36 +0200
Subject: [PATCH 04/15] Fix failsafe test

---
 testing/test_optimizers.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/testing/test_optimizers.py b/testing/test_optimizers.py
index 88757b4..17380f6 100644
--- a/testing/test_optimizers.py
+++ b/testing/test_optimizers.py
@@ -247,6 +247,7 @@ def test_failsafe(self):
         with self.assertRaises(RuntimeError) as e:
             with self.optimizer.failsafe():
                 self.optimizer.acquisition.set_data(X, Y)
+                self.optimizer.acquisition.evaluate(X)
 
         fname = 'failed_bopt_{0}.npz'.format(id(e.exception))
         self.assertTrue(os.path.isfile(fname))

From e4dd2991b6f87f0978f7b4f60f97d3ccb67a915e Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Tue, 15 Aug 2017 14:04:28 +0200
Subject: [PATCH 05/15] Add setting domain functionality in BayesianOptimizer

---
 GPflowOpt/bo.py            | 11 ++++++++++-
 testing/test_optimizers.py | 10 ++++++++++
 2 files changed, 20 insertions(+), 1 deletion(-)

diff --git a/GPflowOpt/bo.py b/GPflowOpt/bo.py
index c54f4f9..d153d90 100644
--- a/GPflowOpt/bo.py
+++ b/GPflowOpt/bo.py
@@ -58,8 +58,10 @@ def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=Tr
         assert initial is None or isinstance(initial, Design)
         super(BayesianOptimizer, self).__init__(domain, exclude_gradient=True)
 
-        if scaling:
+        self._scaling = scaling
+        if self._scaling:
             acquisition.enable_scaling(domain)
+
         self.acquisition = acquisition if hyper_draws is None else MCMCAcquistion(acquisition, hyper_draws)
 
         self.optimizer = optimizer or SciPyOptimizer(domain)
@@ -67,6 +69,13 @@ def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=Tr
         initial = initial or EmptyDesign(domain)
         self.set_initial(initial.generate())
 
+    @Optimizer.domain.setter
+    def domain(self, dom):
+        assert (self.domain.size == dom.size)
+        super(BayesianOptimizer, self.__class__).domain.fset(self, dom)
+        if self._scaling:
+            self.acquisition.enable_scaling(dom)
+
     def _update_model_data(self, newX, newY):
         """
         Update the underlying models of the acquisition function with new data.
diff --git a/testing/test_optimizers.py b/testing/test_optimizers.py
index 2bd6599..47fa29b 100644
--- a/testing/test_optimizers.py
+++ b/testing/test_optimizers.py
@@ -262,6 +262,16 @@ def test_failsafe(self):
             np.testing.assert_almost_equal(data['Y'], Y)
         os.remove(fname)
 
+    def test_set_domain(self):
+        with self.assertRaises(AssertionError):
+            super(TestBayesianOptimizer, self).test_set_domain()
+
+        domain = GPflowOpt.domain.ContinuousParameter("x1", -2.0, 2.0) + \
+                 GPflowOpt.domain.ContinuousParameter("x2", -2.0, 2.0)
+        self.optimizer.domain = domain
+        expected = GPflowOpt.design.LatinHyperCube(16, self.domain).generate() / 4 + 0.5
+        self.assertTrue(np.allclose(expected, self.optimizer.acquisition.models[0].wrapped.X.value))
+
 
 class TestBayesianOptimizerConfigurations(unittest.TestCase):
     def setUp(self):

From 73ef916588a87b48bf5d09549ad3d5a4ec4232a9 Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Tue, 15 Aug 2017 15:11:31 +0200
Subject: [PATCH 06/15] Increasing test coverage (some lines were missing after
 the refactor)

---
 testing/test_acquisition.py     | 15 +++++++++
 testing/test_datascaler.py      |  3 ++
 testing/test_implementations.py | 57 +++++++++++++++++++++++++++++++++
 3 files changed, 75 insertions(+)

diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py
index 7adf882..c7a54b3 100644
--- a/testing/test_acquisition.py
+++ b/testing/test_acquisition.py
@@ -24,6 +24,8 @@ def build_acquisition(self, Xcand):
 
 class TestAcquisition(unittest.TestCase):
 
+    _multiprocess_can_split_ = True
+
     def setUp(self):
         self.model = create_parabola_model(domain)
         self.acquisition = SimpleAcquisition(self.model)
@@ -110,6 +112,16 @@ def test_result_shape_np(self):
         self.assertTupleEqual(res[0].shape, (50, 1))
         self.assertTupleEqual(res[1].shape, (50, domain.size))
 
+    def test_optimize(self):
+        self.acquisition.optimize_restarts = 0
+        state = self.acquisition.get_free_state()
+        self.acquisition._optimize_models()
+        self.assertTrue(np.allclose(state, self.acquisition.get_free_state()))
+
+        self.acquisition.optimize_restarts = 1
+        self.acquisition._optimize_models()
+        self.assertFalse(np.allclose(state, self.acquisition.get_free_state()))
+
 
 aggregations = list()
 aggregations.append(GPflowOpt.acquisition.AcquisitionSum([
@@ -126,6 +138,9 @@ def test_result_shape_np(self):
 
 
 class TestAcquisitionAggregation(unittest.TestCase):
+
+    _multiprocess_can_split_ = True
+
     @parameterized.expand(list(zip(aggregations)))
     def test_object_integrity(self, acquisition):
         for oper in acquisition.operands:
diff --git a/testing/test_datascaler.py b/testing/test_datascaler.py
index ba91577..9927c03 100644
--- a/testing/test_datascaler.py
+++ b/testing/test_datascaler.py
@@ -10,6 +10,9 @@ def parabola2d(X):
 
 
 class TestDataScaler(unittest.TestCase):
+
+    _multiprocess_can_split_ = True
+
     @property
     def domain(self):
         return np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
diff --git a/testing/test_implementations.py b/testing/test_implementations.py
index fa850e0..e65eeb4 100644
--- a/testing/test_implementations.py
+++ b/testing/test_implementations.py
@@ -1,10 +1,45 @@
 import unittest
 import GPflowOpt
 import numpy as np
+from parameterized import parameterized
 from .utility import create_parabola_model, create_plane_model, create_vlmop2_model, parabola2d, load_data
 
+domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
+acquisitions = [GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain)),
+                GPflowOpt.acquisition.ProbabilityOfImprovement(create_parabola_model(domain)),
+                GPflowOpt.acquisition.ProbabilityOfFeasibility(create_parabola_model(domain)),
+                GPflowOpt.acquisition.LowerConfidenceBound(create_parabola_model(domain)),
+                GPflowOpt.acquisition.HVProbabilityOfImprovement([create_parabola_model(domain),
+                                                                  create_parabola_model(domain)])
+                ]
+
+
+class TestAcquisitionEvaluate(unittest.TestCase):
+
+    _multiprocess_can_split_ = True
+
+    @parameterized.expand(list(zip(acquisitions)))
+    def test_evaluate(self, acquisition):
+        X = GPflowOpt.design.RandomDesign(10, domain).generate()
+        p = acquisition.evaluate(X)
+        self.assertTrue(isinstance(p, np.ndarray))
+        self.assertTupleEqual(p.shape, (10, 1))
+
+        q = acquisition.evaluate_with_gradients(X)
+        self.assertTrue(isinstance(q, tuple))
+        self.assertTrue(len(q), 2)
+        for i in q:
+            self.assertTrue(isinstance(i, np.ndarray))
+
+        self.assertTupleEqual(q[0].shape, (10, 1))
+        self.assertTrue(np.allclose(p, q[0]))
+        self.assertTupleEqual(q[1].shape, (10, 2))
+
 
 class TestExpectedImprovement(unittest.TestCase):
+
+    _multiprocess_can_split_ = True
+
     def setUp(self):
         self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
         self.model = create_parabola_model(self.domain)
@@ -42,6 +77,9 @@ def test_EI_validity(self):
 
 
 class TestProbabilityOfImprovement(unittest.TestCase):
+
+    _multiprocess_can_split_ = True
+
     def setUp(self):
         self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
         self.model = create_parabola_model(self.domain)
@@ -66,8 +104,22 @@ def test_setup(self):
         self.acquisition.setup()
         self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated")
 
+    def test_PoI_validity(self):
+        Xcenter = np.random.rand(20, 2) * 0.25 - 0.125
+        X = np.random.rand(100, 2) * 2 - 1
+        hor_idx = np.abs(X[:, 0]) > 0.8
+        ver_idx = np.abs(X[:, 1]) > 0.8
+        Xborder = np.vstack((X[hor_idx, :], X[ver_idx, :]))
+        poi1 = self.acquisition.evaluate(Xborder)
+        poi2 = self.acquisition.evaluate(Xcenter)
+        self.assertGreater(np.min(poi2), np.max(poi1))
+        self.assertTrue(np.all(self.acquisition.feasible_data_index()), msg="EI does never invalidate points")
+
 
 class TestProbabilityOfFeasibility(unittest.TestCase):
+
+    _multiprocess_can_split_ = True
+
     def setUp(self):
         self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
         self.model = create_plane_model(self.domain)
@@ -86,6 +138,9 @@ def test_PoF_validity(self):
 
 
 class TestLowerConfidenceBound(unittest.TestCase):
+
+    _multiprocess_can_split_ = True
+
     def setUp(self):
         self.domain = np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)])
         self.model = create_plane_model(self.domain)
@@ -114,6 +169,8 @@ def test_LCB_validity_2(self):
 
 class TestHVProbabilityOfImprovement(unittest.TestCase):
 
+    _multiprocess_can_split_ = True
+
     def setUp(self):
         self.model = create_vlmop2_model()
         self.data = load_data('vlmop.npz')

From 1dbfe746d235879f41e4696faa823b6c11a514de Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Sat, 19 Aug 2017 15:19:16 +0200
Subject: [PATCH 07/15] Fixing failing tests after merge

---
 GPflowOpt/acquisition/acquisition.py | 1 +
 testing/test_acquisition.py          | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index b71d604..8b3ce68 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -30,6 +30,7 @@
 
 
 class Setup(object):
+        
     def __call__(self, af_method):
         @wraps(af_method)
         def runnable(*args, **kwargs):
diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py
index c7a54b3..cb38ddc 100644
--- a/testing/test_acquisition.py
+++ b/testing/test_acquisition.py
@@ -146,7 +146,7 @@ def test_object_integrity(self, acquisition):
         for oper in acquisition.operands:
             self.assertTrue(isinstance(oper, GPflowOpt.acquisition.Acquisition),
                             msg="All operands should be an acquisition object")
-        self.assertTrue(all(isinstance(m, GPflow.model.GPModel) for m in acquisition.models.sorted_params))
+        self.assertTrue(all(isinstance(m, GPflowOpt.models.ModelWrapper) for m in acquisition.models.sorted_params))
 
     @parameterized.expand(list(zip(aggregations)))
     def test_data(self, acquisition):

From 1129daad02c9de4621e1a6243a427caa9b27e55f Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Sat, 19 Aug 2017 16:10:07 +0200
Subject: [PATCH 08/15] Fixed decorator to a normal function

---
 GPflowOpt/acquisition/acquisition.py | 38 +++++++++++++---------------
 1 file changed, 18 insertions(+), 20 deletions(-)

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index 8b3ce68..63a84a3 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -29,24 +29,22 @@
 float_type = settings.dtypes.float_type
 
 
-class Setup(object):
-        
-    def __call__(self, af_method):
-        @wraps(af_method)
-        def runnable(*args, **kwargs):
-            hp = args[0].highest_parent
-            if hp._needs_setup:
-                # 1 - optimize
-                hp._optimize_models()
-                # 2 - setup
-                # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition
-                # e.g. through feasible_data_index.
-                hp._needs_setup = False
-                hp.setup()
-            results = af_method(*args, **kwargs)
-            return results
-
-        return runnable
+def setup_required(af_method):
+    @wraps(af_method)
+    def runnable(*args, **kwargs):
+        hp = args[0].highest_parent
+        if hp._needs_setup:
+            # 1 - optimize
+            hp._optimize_models()
+            # 2 - setup
+            # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition
+            # e.g. through feasible_data_index.
+            hp._needs_setup = False
+            hp.setup()
+        results = af_method(*args, **kwargs)
+        return results
+
+    return runnable
 
 
 class Acquisition(Parameterized):
@@ -212,7 +210,7 @@ def setup(self):
         """
         pass
 
-    @Setup()
+    @setup_required
     @AutoFlow((float_type, [None, None]))
     def evaluate_with_gradients(self, Xcand):
         """
@@ -224,7 +222,7 @@ def evaluate_with_gradients(self, Xcand):
         acq = self.build_acquisition(Xcand)
         return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0]
 
-    @Setup()
+    @setup_required
     @AutoFlow((float_type, [None, None]))
     def evaluate(self, Xcand):
         """

From 961e09c4e98831fc2c856fe160f155dcf43d660b Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Sat, 19 Aug 2017 22:34:16 +0200
Subject: [PATCH 09/15] Added some documentation, + small fix

---
 GPflowOpt/acquisition/acquisition.py | 35 ++++++++++++++++++----------
 1 file changed, 23 insertions(+), 12 deletions(-)

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index 63a84a3..80a986a 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -29,10 +29,15 @@
 float_type = settings.dtypes.float_type
 
 
-def setup_required(af_method):
-    @wraps(af_method)
-    def runnable(*args, **kwargs):
-        hp = args[0].highest_parent
+def setup_required(method):
+    """
+    Decorator function to mark methods in Acquisition classes which require running setup if indicated by _needs_setup
+    :param method: acquisition method
+    """
+    @wraps(method)
+    def runnable(instance, *args, **kwargs):
+        assert isinstance(instance, Acquisition)
+        hp = instance.highest_parent
         if hp._needs_setup:
             # 1 - optimize
             hp._optimize_models()
@@ -41,7 +46,7 @@ def runnable(*args, **kwargs):
             # e.g. through feasible_data_index.
             hp._needs_setup = False
             hp.setup()
-        results = af_method(*args, **kwargs)
+        results = method(instance, *args, **kwargs)
         return results
 
     return runnable
@@ -53,14 +58,20 @@ class Acquisition(Parameterized):
     score indicating how promising a point is for evaluation.
 
     In Bayesian Optimization this function is typically optimized over the optimization domain
-    to determine the next point for evaluation.
+    to determine the next point for evaluation. An object of this class holds a list of GPflow models. Subclasses
+    implement a build_acquisition function which computes the acquisition function (usually from the predictive
+    distribution) using TensorFlow. Optionally, a method setup can be implemented which computes some quantities which
+    are used to compute the acquisition, but do not depend on candidate points.
+
+    Acquisition functions can be combined through addition or multiplication to construct joint criteria. For instance,
+    for constrained optimization. The objects then form a tree hierarchy.
 
-    An object of this class holds a list of GPflow models. Subclasses implement a build_acquisition function
-    which computes the acquisition function (usually from the predictive distribution) using TensorFlow.
-    Each model is automatically optimized when an acquisition object is constructed or when set_data is called.
+    Acquisition models implement a lazy strategy to optimize models and run setup. This is implemented by a _needs_setup
+    attribute (similar to the _needs_recompile in GPflow). Calling set_data sets this flag to true. Calling methods
+    marked with the setup_require decorator (such as evaluate) optimize all models, then call setup if this flag is set.
+    In hierarchies, first acquisition objects handling constraint objectives are set up, then the objects handling
+    objectives.
 
-    Acquisition functions can be combined through addition or multiplication to construct joint criteria. 
-    For instance, for constrained optimization.
     """
 
     def __init__(self, models=[], optimize_restarts=5):
@@ -261,7 +272,7 @@ def __mul__(self, other):
     def __setattr__(self, key, value):
         super(Acquisition, self).__setattr__(key, value)
         if key is '_parent':
-            self._needs_setup = True
+            self.highest_parent._needs_setup = True
 
 
 class AcquisitionAggregation(Acquisition):

From 0479807f6a3636d0040f6e8cacbafc85cec0dfcc Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Mon, 21 Aug 2017 14:29:53 +0200
Subject: [PATCH 10/15] Reworked some testing code

---
 testing/test_modelwrapper.py | 49 ++++++++++++++++--------------------
 testing/test_optimizers.py   | 37 ++++++---------------------
 2 files changed, 29 insertions(+), 57 deletions(-)

diff --git a/testing/test_modelwrapper.py b/testing/test_modelwrapper.py
index 439858e..0046627 100644
--- a/testing/test_modelwrapper.py
+++ b/testing/test_modelwrapper.py
@@ -2,6 +2,7 @@
 import unittest
 import GPflow
 import numpy as np
+from .utility import create_parabola_model
 
 float_type = GPflow.settings.dtypes.float_type
 
@@ -39,55 +40,47 @@ def foo(self, val):
 
 class TestModelWrapper(unittest.TestCase):
 
-    def simple_model(self):
-        x = np.random.rand(10,2) * 2 * np.pi
-        y = np.sin(x[:,[0]])
-        m = GPflow.gpr.GPR(x,y, kern=GPflow.kernels.RBF(1))
-        return m
+    def setUp(self):
+        self.m = create_parabola_model(GPflowOpt.domain.UnitCube(2))
 
     def test_object_integrity(self):
-        m = self.simple_model()
-        w = GPflowOpt.models.ModelWrapper(m)
-        self.assertEqual(w.wrapped, m)
-        self.assertEqual(m._parent, w)
-        self.assertEqual(w.optimize, m.optimize)
+        w = GPflowOpt.models.ModelWrapper(self.m)
+        self.assertEqual(w.wrapped, self.m)
+        self.assertEqual(self.m._parent, w)
+        self.assertEqual(w.optimize, self.m.optimize)
 
     def test_optimize(self):
-        m = self.simple_model()
-        w = GPflowOpt.models.ModelWrapper(m)
-        logL = m.compute_log_likelihood()
+        w = GPflowOpt.models.ModelWrapper(self.m)
+        logL = self.m.compute_log_likelihood()
         self.assertTrue(np.allclose(logL, w.compute_log_likelihood()))
 
         # Check if compiled & optimized, verify attributes are set in the right object.
         w.optimize(maxiter=5)
-        self.assertTrue(hasattr(m, '_minusF'))
+        self.assertTrue(hasattr(self.m, '_minusF'))
         self.assertFalse('_minusF' in w.__dict__)
-        self.assertGreater(m.compute_log_likelihood(), logL)
+        self.assertGreater(self.m.compute_log_likelihood(), logL)
 
     def test_af_storage_detection(self):
         # Regression test for a bug with predict_f/predict_y... etc.
-        m = self.simple_model()
         x = np.random.rand(10,2)
-        m.predict_f(x)
-        self.assertTrue(hasattr(m, '_predict_f_AF_storage'))
-        w = MethodOverride(m)
+        self.m.predict_f(x)
+        self.assertTrue(hasattr(self.m, '_predict_f_AF_storage'))
+        w = MethodOverride(self.m)
         self.assertFalse(hasattr(w, '_predict_f_AF_storage'))
         w.predict_f(x)
         self.assertTrue(hasattr(w, '_predict_f_AF_storage'))
 
     def test_set_wrapped_attributes(self):
         # Regression test for setting certain keys in the right object
-        m = self.simple_model()
-        w = GPflowOpt.models.ModelWrapper(m)
+        w = GPflowOpt.models.ModelWrapper(self.m)
         w._needs_recompile = False
         self.assertFalse('_needs_recompile' in w.__dict__)
-        self.assertTrue('_needs_recompile' in m.__dict__)
+        self.assertTrue('_needs_recompile' in self.m.__dict__)
         self.assertFalse(w._needs_recompile)
-        self.assertFalse(m._needs_recompile)
+        self.assertFalse(self.m._needs_recompile)
 
     def test_double_wrap(self):
-        m = self.simple_model()
-        n = GPflowOpt.models.ModelWrapper(MethodOverride(m))
+        n = GPflowOpt.models.ModelWrapper(MethodOverride(self.m))
         n.optimize(maxiter=10)
         Xt = np.random.rand(10, 2)
         n.predict_f(Xt)
@@ -95,7 +88,7 @@ def test_double_wrap(self):
         self.assertTrue('_predict_f_AF_storage' in n.wrapped.__dict__)
         self.assertFalse('_predict_f_AF_storage' in n.wrapped.wrapped.__dict__)
 
-        n = MethodOverride(GPflowOpt.models.ModelWrapper(m))
+        n = MethodOverride(GPflowOpt.models.ModelWrapper(self.m))
         Xn = np.random.rand(10, 2)
         Yn = np.random.rand(10, 1)
         n.X = Xn
@@ -110,12 +103,12 @@ def test_double_wrap(self):
         self.assertFalse('foo' in n.wrapped.wrapped.__dict__)
 
     def test_name(self):
-        n = GPflowOpt.models.ModelWrapper(self.simple_model())
+        n = GPflowOpt.models.ModelWrapper(self.m)
         self.assertEqual(n.name, 'unnamed.modelwrapper')
         p = GPflow.param.Parameterized()
         p.model = n
         self.assertEqual(n.name, 'model.modelwrapper')
-        n = MethodOverride(self.simple_model())
+        n = MethodOverride(create_parabola_model(GPflowOpt.domain.UnitCube(2)))
         self.assertEqual(n.name, 'unnamed.methodoverride')
 
 
diff --git a/testing/test_optimizers.py b/testing/test_optimizers.py
index 47fa29b..87ce3db 100644
--- a/testing/test_optimizers.py
+++ b/testing/test_optimizers.py
@@ -8,21 +8,13 @@
 import warnings
 from contextlib import contextmanager
 from scipy.optimize import OptimizeResult
+from .utility import vlmop2, create_parabola_model, create_vlmop2_model
 
 
 def parabola2d(X):
     return np.atleast_2d(np.sum(X ** 2, axis=1)).T, 2 * X
 
 
-def vlmop2(x):
-    transl = 1 / np.sqrt(2)
-    part1 = (x[:, [0]] - transl) ** 2 + (x[:, [1]] - transl) ** 2
-    part2 = (x[:, [0]] + transl) ** 2 + (x[:, [1]] + transl) ** 2
-    y1 = 1 - np.exp(-1 * part1)
-    y2 = 1 - np.exp(-1 * part2)
-    return np.hstack((y1, y2))
-
-
 class KeyboardRaiser:
     """
     This wraps a function and makes it raise a KeyboardInterrupt after some number of calls
@@ -202,21 +194,9 @@ def test_set_domain(self):
 class TestBayesianOptimizer(_TestOptimizer, unittest.TestCase):
     def setUp(self):
         super(TestBayesianOptimizer, self).setUp()
-        design = GPflowOpt.design.LatinHyperCube(16, self.domain)
-        X, Y = design.generate(), parabola2d(design.generate())[0]
-        model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True))
-        acquisition = GPflowOpt.acquisition.ExpectedImprovement(model)
+        acquisition = GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(self.domain))
         self.optimizer = GPflowOpt.BayesianOptimizer(self.domain, acquisition)
 
-    def setup_multi_objective(self):
-        design = GPflowOpt.design.LatinHyperCube(16, self.domain)
-        X = design.generate()
-        Y = vlmop2(X)
-        m1 = GPflow.gpr.GPR(X, Y[:,[0]], GPflow.kernels.RBF(2, ARD=True))
-        m2 = GPflow.gpr.GPR(X.copy(), Y[:,[1]], GPflow.kernels.RBF(2, ARD=True))
-        acquisition = GPflowOpt.acquisition.ExpectedImprovement(m1) + GPflowOpt.acquisition.ExpectedImprovement(m2)
-        return GPflowOpt.BayesianOptimizer(self.domain, acquisition)
-
     def test_default_initial(self):
         self.assertTupleEqual(self.optimizer._initial.shape, (0, 2), msg="Invalid shape of initial points array")
 
@@ -228,12 +208,14 @@ def test_optimize(self):
         self.assertTrue(np.allclose(result.fun, 0), msg="Incorrect function value returned")
 
     def test_optimize_multi_objective(self):
-        optimizer = self.setup_multi_objective()
+        m1, m2 = create_vlmop2_model()
+        acquisition = GPflowOpt.acquisition.ExpectedImprovement(m1) + GPflowOpt.acquisition.ExpectedImprovement(m2)
+        optimizer = GPflowOpt.BayesianOptimizer(self.domain, acquisition)
         result = optimizer.optimize(vlmop2, n_iter=2)
         self.assertTrue(result.success)
         self.assertEqual(result.nfev, 2, "Only 2 evaluations permitted")
-        self.assertTupleEqual(result.x.shape, (8, 2))
-        self.assertTupleEqual(result.fun.shape, (8, 2))
+        self.assertTupleEqual(result.x.shape, (9, 2))
+        self.assertTupleEqual(result.fun.shape, (9, 2))
         _, dom = GPflowOpt.pareto.non_dominated_sort(result.fun)
         self.assertTrue(np.all(dom==0))
 
@@ -277,10 +259,7 @@ class TestBayesianOptimizerConfigurations(unittest.TestCase):
     def setUp(self):
         self.domain = GPflowOpt.domain.ContinuousParameter("x1", 0.0, 1.0) + \
                       GPflowOpt.domain.ContinuousParameter("x2", 0.0, 1.0)
-        design = GPflowOpt.design.LatinHyperCube(16, self.domain)
-        X, Y = design.generate(), parabola2d(design.generate())[0]
-        model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.var(axis=0)))
-        self.acquisition = GPflowOpt.acquisition.ExpectedImprovement(model)
+        self.acquisition = GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(self.domain))
 
     def test_initial_design(self):
         design = GPflowOpt.design.RandomDesign(5, self.domain)

From 8cd806dbb1ff38183296cf0129d3fb3dd7914208 Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Tue, 22 Aug 2017 14:47:38 +0200
Subject: [PATCH 11/15] Improved design for setup objectives vs setup
 constraints

---
 GPflowOpt/acquisition/acquisition.py | 21 +++++++++++++++++----
 GPflowOpt/acquisition/hvpoi.py       |  2 +-
 2 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index 80a986a..6b44a3d 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -221,6 +221,20 @@ def setup(self):
         """
         pass
 
+    def _setup_constraints(self):
+        """
+        Run only if some outputs handled by this acquisition are constraints. Used in aggregation.
+        """
+        if self.constraint_indices().size > 0:
+            self.setup()
+
+    def _setup_objectives(self):
+        """
+        Run only if all outputs handled by this acquisition are objectives. Used in aggregation.
+        """
+        if self.constraint_indices().size == 0:
+            self.setup()
+
     @setup_required
     @AutoFlow((float_type, [None, None]))
     def evaluate_with_gradients(self, Xcand):
@@ -310,13 +324,12 @@ def set_data(self, X, Y):
 
     def _setup_constraints(self):
         for oper in self.operands:
-            if oper.constraint_indices().size > 0:
-                oper._setup_constraints() if isinstance(oper, AcquisitionAggregation) else oper.setup()
+            if oper.constraint_indices().size > 0:  # Small optimization, skip subtrees with objectives only
+                oper._setup_constraints()
 
     def _setup_objectives(self):
         for oper in self.operands:
-            if oper.constraint_indices().size == 0:
-                oper._setup_objectives() if isinstance(oper, AcquisitionAggregation) else oper.setup()
+            oper._setup_objectives()
 
     def setup(self):
         # First setup acquisitions involving constraints
diff --git a/GPflowOpt/acquisition/hvpoi.py b/GPflowOpt/acquisition/hvpoi.py
index 1c310ea..7d68bd6 100644
--- a/GPflowOpt/acquisition/hvpoi.py
+++ b/GPflowOpt/acquisition/hvpoi.py
@@ -72,7 +72,7 @@ def __init__(self, models):
         num_objectives = self.data[1].shape[1]
         assert num_objectives > 1
 
-        # Keep pareto empty for now - its updated in setup()
+        # Keep pareto empty for now - it is updated in setup()
         self.pareto = Pareto(np.empty((0, num_objectives)))
         self.reference = DataHolder(np.ones((1, num_objectives)))
 

From 9786bfa5438ad2cd8d64b5ae18fc75244d2ad514 Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Tue, 22 Aug 2017 14:53:37 +0200
Subject: [PATCH 12/15] Moving reset of setup flag

---
 GPflowOpt/acquisition/acquisition.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index 6b44a3d..c544174 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -39,12 +39,12 @@ def runnable(instance, *args, **kwargs):
         assert isinstance(instance, Acquisition)
         hp = instance.highest_parent
         if hp._needs_setup:
+            hp._needs_setup = False
             # 1 - optimize
             hp._optimize_models()
             # 2 - setup
             # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition
             # e.g. through feasible_data_index.
-            hp._needs_setup = False
             hp.setup()
         results = method(instance, *args, **kwargs)
         return results

From ad83c12a192eb7337edddf5e31571201b471323d Mon Sep 17 00:00:00 2001
From: Joachim van der Herten <joachim.vanderherten@ugent.be>
Date: Tue, 22 Aug 2017 23:19:02 +0200
Subject: [PATCH 13/15] Made setup method in acquisition private

---
 GPflowOpt/acquisition/acquisition.py | 15 +++++++++------
 GPflowOpt/acquisition/ei.py          |  6 +++---
 GPflowOpt/acquisition/hvpoi.py       |  4 ++--
 GPflowOpt/acquisition/poi.py         |  6 +++---
 testing/test_acquisition.py          | 10 +++++-----
 testing/test_implementations.py      |  8 ++++----
 6 files changed, 26 insertions(+), 23 deletions(-)

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index c544174..68e1137 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -45,7 +45,7 @@ def runnable(instance, *args, **kwargs):
             # 2 - setup
             # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition
             # e.g. through feasible_data_index.
-            hp.setup()
+            hp._setup()
         results = method(instance, *args, **kwargs)
         return results
 
@@ -213,11 +213,14 @@ def feasible_data_index(self):
         """
         return np.ones(self.data[0].shape[0], dtype=bool)
 
-    def setup(self):
+    def _setup(self):
         """
         Pre-calculation of quantities used later in the evaluation of the acquisition function for candidate points.
         
-        Automatically triggered by :meth:`~.Acquisition.set_data`.
+        Subclasses can implement this method to compute quantites (such as fmin). The decision when to run this function
+        is governed by the acquisition class, based on the setup_required decorator on methods and methods which require
+        setup to be run (e.g. set_data). This method shouldn't be called directly as the results are likely to be
+        overwritten at a later point when the Acquisition object decides to run the method.
         """
         pass
 
@@ -226,14 +229,14 @@ def _setup_constraints(self):
         Run only if some outputs handled by this acquisition are constraints. Used in aggregation.
         """
         if self.constraint_indices().size > 0:
-            self.setup()
+            self._setup()
 
     def _setup_objectives(self):
         """
         Run only if all outputs handled by this acquisition are objectives. Used in aggregation.
         """
         if self.constraint_indices().size == 0:
-            self.setup()
+            self._setup()
 
     @setup_required
     @AutoFlow((float_type, [None, None]))
@@ -331,7 +334,7 @@ def _setup_objectives(self):
         for oper in self.operands:
             oper._setup_objectives()
 
-    def setup(self):
+    def _setup(self):
         # First setup acquisitions involving constraints
         self._setup_constraints()
         # Then objectives
diff --git a/GPflowOpt/acquisition/ei.py b/GPflowOpt/acquisition/ei.py
index 2c9ce87..77daf0a 100644
--- a/GPflowOpt/acquisition/ei.py
+++ b/GPflowOpt/acquisition/ei.py
@@ -58,10 +58,10 @@ def __init__(self, model):
         """
         super(ExpectedImprovement, self).__init__(model)
         self.fmin = DataHolder(np.zeros(1))
-        self.setup()
+        self._setup()
 
-    def setup(self):
-        super(ExpectedImprovement, self).setup()
+    def _setup(self):
+        super(ExpectedImprovement, self)._setup()
         # 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)
diff --git a/GPflowOpt/acquisition/hvpoi.py b/GPflowOpt/acquisition/hvpoi.py
index 7d68bd6..e5e5033 100644
--- a/GPflowOpt/acquisition/hvpoi.py
+++ b/GPflowOpt/acquisition/hvpoi.py
@@ -81,11 +81,11 @@ def _estimate_reference(self):
         f = np.max(pf, axis=0, keepdims=True) - np.min(pf, axis=0, keepdims=True)
         return np.max(pf, axis=0, keepdims=True) + 2 * f / pf.shape[0]
 
-    def setup(self):
+    def _setup(self):
         """
         Pre-computes the Pareto set and cell bounds for integrating over the non-dominated region.
         """
-        super(HVProbabilityOfImprovement, self).setup()
+        super(HVProbabilityOfImprovement, self)._setup()
 
         # Obtain hypervolume cell bounds, use prediction mean
         feasible_samples = self.data[0][self.highest_parent.feasible_data_index(), :]
diff --git a/GPflowOpt/acquisition/poi.py b/GPflowOpt/acquisition/poi.py
index 707a59e..aff1089 100644
--- a/GPflowOpt/acquisition/poi.py
+++ b/GPflowOpt/acquisition/poi.py
@@ -37,10 +37,10 @@ def __init__(self, model):
         """
         super(ProbabilityOfImprovement, self).__init__(model)
         self.fmin = DataHolder(np.zeros(1))
-        self.setup()
+        self._setup()
 
-    def setup(self):
-        super(ProbabilityOfImprovement, self).setup()
+    def _setup(self):
+        super(ProbabilityOfImprovement, self)._setup()
         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))
diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py
index cb38ddc..999f2ae 100644
--- a/testing/test_acquisition.py
+++ b/testing/test_acquisition.py
@@ -14,8 +14,8 @@ def __init__(self, model):
         super(SimpleAcquisition, self).__init__(model)
         self.counter = 0
 
-    def setup(self):
-        super(SimpleAcquisition, self).setup()
+    def _setup(self):
+        super(SimpleAcquisition, self)._setup()
         self.counter += 1
 
     def build_acquisition(self, Xcand):
@@ -33,7 +33,7 @@ def setUp(self):
     def run_setup(self):
         # Optimize models & perform acquisition setup call.
         self.acquisition._optimize_models()
-        self.acquisition.setup()
+        self.acquisition._setup()
 
     def test_object_integrity(self):
         self.assertEqual(len(self.acquisition.models), 1, msg="Model list has incorrect length.")
@@ -211,7 +211,7 @@ def test_hyper_updates(self, acquisition):
     @parameterized.expand(list(zip([aggregations[2]])))
     def test_marginalized_score(self, acquisition):
         acquisition._optimize_models()
-        acquisition.setup()
+        acquisition._setup()
         Xt = np.random.rand(20, 2) * 2 - 1
         ei_mle = acquisition.operands[0].evaluate(Xt)
         ei_mcmc = acquisition.evaluate(Xt)
@@ -243,7 +243,7 @@ def test_constrained_EI(self):
 
         # Test proper setup
         joint._optimize_models()
-        joint.setup()
+        joint._setup()
         self.assertGreater(ei.fmin.value, np.min(ei.data[1]), msg="The best objective value is in an infeasible area")
         self.assertTrue(np.allclose(ei.fmin.value, np.min(ei.data[1][pof.feasible_data_index(), :]), atol=1e-3),
                         msg="fmin computed incorrectly")
diff --git a/testing/test_implementations.py b/testing/test_implementations.py
index e65eeb4..f907d25 100644
--- a/testing/test_implementations.py
+++ b/testing/test_implementations.py
@@ -51,7 +51,7 @@ def test_objective_indices(self):
 
     def test_setup(self):
         self.acquisition._optimize_models()
-        self.acquisition.setup()
+        self.acquisition._setup()
         fmin = np.min(self.acquisition.data[1])
         self.assertGreater(self.acquisition.fmin.value, 0, msg="The minimum (0) is not amongst the design.")
         self.assertTrue(np.allclose(self.acquisition.fmin.value, fmin, atol=1e-2), msg="fmin computed incorrectly")
@@ -61,7 +61,7 @@ def test_setup(self):
         self.acquisition.set_data(np.vstack((self.acquisition.data[0], p)),
                                   np.vstack((self.acquisition.data[1], parabola2d(p))))
         self.acquisition._optimize_models()
-        self.acquisition.setup()
+        self.acquisition._setup()
         self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated")
 
     def test_EI_validity(self):
@@ -91,7 +91,7 @@ def test_objective_indices(self):
 
     def test_setup(self):
         self.acquisition._optimize_models()
-        self.acquisition.setup()
+        self.acquisition._setup()
         fmin = np.min(self.acquisition.data[1])
         self.assertGreater(self.acquisition.fmin.value, 0, msg="The minimum (0) is not amongst the design.")
         self.assertTrue(np.allclose(self.acquisition.fmin.value, fmin, atol=1e-2), msg="fmin computed incorrectly")
@@ -101,7 +101,7 @@ def test_setup(self):
         self.acquisition.set_data(np.vstack((self.acquisition.data[0], p)),
                                   np.vstack((self.acquisition.data[1], parabola2d(p))))
         self.acquisition._optimize_models()
-        self.acquisition.setup()
+        self.acquisition._setup()
         self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated")
 
     def test_PoI_validity(self):

From cd8c392ee254fcf01e3da30c3ecc1b28da686fd5 Mon Sep 17 00:00:00 2001
From: icouckuy <ivo.couckuyt@ugent.be>
Date: Wed, 23 Aug 2017 21:25:07 +0200
Subject: [PATCH 14/15] Doc update to reflect the new lazy setup better

---
 GPflowOpt/acquisition/acquisition.py | 35 ++++++++++++++++------------
 GPflowOpt/acquisition/hvpoi.py       |  2 +-
 GPflowOpt/bo.py                      |  4 ++--
 3 files changed, 23 insertions(+), 18 deletions(-)

diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index 68e1137..6050aa7 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -39,12 +39,14 @@ def runnable(instance, *args, **kwargs):
         assert isinstance(instance, Acquisition)
         hp = instance.highest_parent
         if hp._needs_setup:
+            # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition
+            # e.g. through feasible_data_index.
             hp._needs_setup = False
+
             # 1 - optimize
             hp._optimize_models()
+
             # 2 - setup
-            # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition
-            # e.g. through feasible_data_index.
             hp._setup()
         results = method(instance, *args, **kwargs)
         return results
@@ -67,11 +69,10 @@ class Acquisition(Parameterized):
     for constrained optimization. The objects then form a tree hierarchy.
 
     Acquisition models implement a lazy strategy to optimize models and run setup. This is implemented by a _needs_setup
-    attribute (similar to the _needs_recompile in GPflow). Calling set_data sets this flag to true. Calling methods
+    attribute (similar to the _needs_recompile in GPflow). Calling :meth:`set_data` sets this flag to True. Calling methods
     marked with the setup_require decorator (such as evaluate) optimize all models, then call setup if this flag is set.
     In hierarchies, first acquisition objects handling constraint objectives are set up, then the objects handling
     objectives.
-
     """
 
     def __init__(self, models=[], optimize_restarts=5):
@@ -92,8 +93,8 @@ def _optimize_models(self):
         """
         Optimizes the hyperparameters of all models that the acquisition function is based on.
 
-        It is called automatically during initialization and each time set_data() is called.
-        When using the high-level :class:`..BayesianOptimizer` class calling set_data() is taken care of.
+        It is called automatically during initialization and each time :meth:`set_data` is called.
+        When using the high-level :class:`..BayesianOptimizer` class calling :meth:`set_data` is taken care of.
 
         For each model the hyperparameters of the model at the time it was passed to __init__() are used as initial
         point and optimized. If optimize_restarts is set to >1, additional randomization
@@ -127,6 +128,9 @@ def enable_scaling(self, domain):
         """
         Enables and configures the :class:`.DataScaler` objects wrapping the GP models.
         
+        Sets the _needs_setup attribute to True so the contained models are optimized and :meth:`setup` is run again
+        right before evaluating the :class:`Acquisition` function.
+
         :param domain: :class:`.Domain` object, the input transform of the data scalers is configured as a transform
             from domain to the unit cube with the same dimensionality.
         """
@@ -139,8 +143,10 @@ def enable_scaling(self, domain):
 
     def set_data(self, X, Y):
         """
-        Update the training data of the contained models. Automatically triggers a hyperparameter optimization
-        step by calling _optimize_all() and an update of pre-computed quantities by calling setup().
+        Update the training data of the contained models
+
+        Sets the _needs_setup attribute to True so the contained models are optimized and :meth:`setup` is run again
+        right before evaluating the :class:`Acquisition` function.
 
         Let Q be the the sum of the output dimensions of all contained models, Y should have a minimum of
         Q columns. Only the first Q columns of Y are used while returning the scalar Q
@@ -217,10 +223,9 @@ def _setup(self):
         """
         Pre-calculation of quantities used later in the evaluation of the acquisition function for candidate points.
         
-        Subclasses can implement this method to compute quantites (such as fmin). The decision when to run this function
-        is governed by the acquisition class, based on the setup_required decorator on methods and methods which require
-        setup to be run (e.g. set_data). This method shouldn't be called directly as the results are likely to be
-        overwritten at a later point when the Acquisition object decides to run the method.
+        Subclasses can implement this method to compute quantities (such as fmin). The decision when to run this function
+        is governed by :class:`Acquisition`, based on the setup_required decorator on methods which require
+        setup to be run (e.g. set_data).
         """
         pass
 
@@ -335,9 +340,9 @@ def _setup_objectives(self):
             oper._setup_objectives()
 
     def _setup(self):
-        # First setup acquisitions involving constraints
+        # Important: First setup acquisitions involving constraints
         self._setup_constraints()
-        # Then objectives
+        # Then objectives as these might depend on the constraint acquisition
         self._setup_objectives()
 
     def constraint_indices(self):
@@ -416,7 +421,7 @@ def _optimize_models(self):
         # Sample each model of the acquisition function - results in a list of 2D ndarrays.
         hypers = np.hstack([model.sample(len(self.operands), **self._sample_opt) for model in self.models])
 
-        # Now visit all copies, and set state
+        # Now visit all acquisition copies, and set state
         for idx, draw in enumerate(self.operands):
             draw.set_state(hypers[idx, :])
 
diff --git a/GPflowOpt/acquisition/hvpoi.py b/GPflowOpt/acquisition/hvpoi.py
index e5e5033..ef7dba7 100644
--- a/GPflowOpt/acquisition/hvpoi.py
+++ b/GPflowOpt/acquisition/hvpoi.py
@@ -72,7 +72,7 @@ def __init__(self, models):
         num_objectives = self.data[1].shape[1]
         assert num_objectives > 1
 
-        # Keep pareto empty for now - it is updated in setup()
+        # Keep empty for now - it is updated in _setup()
         self.pareto = Pareto(np.empty((0, num_objectives)))
         self.reference = DataHolder(np.ones((1, num_objectives)))
 
diff --git a/GPflowOpt/bo.py b/GPflowOpt/bo.py
index d153d90..8071104 100644
--- a/GPflowOpt/bo.py
+++ b/GPflowOpt/bo.py
@@ -71,7 +71,7 @@ def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=Tr
 
     @Optimizer.domain.setter
     def domain(self, dom):
-        assert (self.domain.size == dom.size)
+        assert self.domain.size == dom.size
         super(BayesianOptimizer, self.__class__).domain.fset(self, dom)
         if self._scaling:
             self.acquisition.enable_scaling(dom)
@@ -175,7 +175,7 @@ def _optimize(self, fx, n_iter):
         :return: OptimizeResult object
         """
 
-        assert(isinstance(fx, ObjectiveWrapper))
+        assert isinstance(fx, ObjectiveWrapper)
 
         # Evaluate and add the initial design (if any)
         initial = self.get_initial()

From 8b18be224ba6ce5bb4555b33b0aa45a38a08cabe Mon Sep 17 00:00:00 2001
From: icouckuy <ivo.couckuyt@ugent.be>
Date: Wed, 23 Aug 2017 21:48:25 +0200
Subject: [PATCH 15/15] Updated docs in scaling and transforms slightly. Added
 missing NotImplementedError for assign method

---
 GPflowOpt/scaling.py    | 14 ++++++++------
 GPflowOpt/transforms.py |  8 ++++++--
 2 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/GPflowOpt/scaling.py b/GPflowOpt/scaling.py
index a8dee0e..fb6d140 100644
--- a/GPflowOpt/scaling.py
+++ b/GPflowOpt/scaling.py
@@ -86,14 +86,16 @@ def input_transform(self):
     @input_transform.setter
     def input_transform(self, t):
         """
-        Configure a new input transform. Data in the model is automatically updated with the new transform.
-        
+        Configure a new input transform.
+
+        Data in the wrapped model is automatically updated with the new transform.
+
         :param t: :class:`.DataTransform` object: the new input transform.
         """
-        assert(isinstance(t, DataTransform))
-        X = self.X.value
+        assert isinstance(t, DataTransform)
+        X = self.X.value  # unscales the data
         self._input_transform.assign(t)
-        self.X = X
+        self.X = X  # scales the back using the new input transform
 
     @property
     def output_transform(self):
@@ -111,7 +113,7 @@ def output_transform(self, t):
         
         :param t: :class:`.DataTransform` object: the new output transform.
         """
-        assert (isinstance(t, DataTransform))
+        assert isinstance(t, DataTransform)
         Y = self.Y.value
         self._output_transform.assign(t)
         self.Y = Y
diff --git a/GPflowOpt/transforms.py b/GPflowOpt/transforms.py
index c87fac3..984a5b2 100644
--- a/GPflowOpt/transforms.py
+++ b/GPflowOpt/transforms.py
@@ -55,6 +55,9 @@ def backward(self, Y):
         """
         return (~self).forward(Y)
 
+    def assign(self, other):
+        raise NotImplementedError
+
     def __invert__(self):
         """
         Return a :class:`.DataTransform` object implementing the reverse transform V -> U
@@ -138,8 +141,9 @@ def build_backward_variance(self, Yvar):
 
     def assign(self, other):
         """
-        Assign the parameters of another :class:`LinearTransform`. Can be useful to avoid graph
-        re-compilation.
+        Assign the parameters of another :class:`LinearTransform`.
+
+        Useful to avoid graph re-compilation.
 
         :param other: :class:`.LinearTransform` object
         """