diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py
index 69d4b80..6050aa7 100644
--- a/GPflowOpt/acquisition/acquisition.py
+++ b/GPflowOpt/acquisition/acquisition.py
@@ -24,24 +24,55 @@
 import tensorflow as tf
 
 import copy
+from functools import wraps
 
 float_type = settings.dtypes.float_type
 
 
+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:
+            # 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
+            hp._setup()
+        results = method(instance, *args, **kwargs)
+        return results
+
+    return runnable
+
+
 class Acquisition(Parameterized):
     """
     An acquisition function maps the belief represented by a Bayesian model into a
     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.
-
-    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 functions can be combined through addition or multiplication to construct joint criteria. 
-    For instance, for constrained optimization.
+    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.
+
+    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 :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):
@@ -56,14 +87,14 @@ 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):
         """
         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
@@ -97,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.
         """
@@ -105,12 +139,14 @@ 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):
         """
-        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
@@ -128,11 +164,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
@@ -187,14 +219,31 @@ 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 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
 
+    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):
         """
@@ -206,6 +255,7 @@ def evaluate_with_gradients(self, Xcand):
         acq = self.build_acquisition(Xcand)
         return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0]
 
+    @setup_required
     @AutoFlow((float_type, [None, None]))
     def evaluate(self, Xcand):
         """
@@ -241,6 +291,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.highest_parent._needs_setup = True
+
 
 class AcquisitionAggregation(Acquisition):
     """
@@ -256,10 +311,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):
@@ -273,13 +328,22 @@ 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:  # Small optimization, skip subtrees with objectives only
+                oper._setup_constraints()
+
+    def _setup_objectives(self):
+        for oper in self.operands:
+            oper._setup_objectives()
+
+    def _setup(self):
+        # Important: First setup acquisitions involving constraints
+        self._setup_constraints()
+        # Then objectives as these might depend on the constraint acquisition
+        self._setup_objectives()
 
     def constraint_indices(self):
         offset = [0]
@@ -304,7 +368,6 @@ class AcquisitionSum(AcquisitionAggregation):
     """
     Sum of acquisition functions
     """
-
     def __init__(self, operands):
         super(AcquisitionSum, self).__init__(operands, tf.reduce_sum)
 
@@ -319,7 +382,6 @@ class AcquisitionProduct(AcquisitionAggregation):
     """
     Product of acquisition functions
     """
-
     def __init__(self, operands):
         super(AcquisitionProduct, self).__init__(operands, tf.reduce_prod)
 
@@ -350,14 +412,16 @@ 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):
+        # Optimize model #1
+        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])
 
-        # 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, :])
 
@@ -371,9 +435,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/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 da40acc..ef7dba7 100644
--- a/GPflowOpt/acquisition/hvpoi.py
+++ b/GPflowOpt/acquisition/hvpoi.py
@@ -69,23 +69,27 @@ 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 empty for now - it is 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
         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
-        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..aff1089 100644
--- a/GPflowOpt/acquisition/poi.py
+++ b/GPflowOpt/acquisition/poi.py
@@ -37,11 +37,12 @@ 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()
-        samples_mean, _ = self.models[0].predict_f(self.data[0])
+    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))
 
     def build_acquisition(self, Xcand):
diff --git a/GPflowOpt/bo.py b/GPflowOpt/bo.py
index ec95058..8071104 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.
@@ -57,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)
@@ -66,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.
@@ -165,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()
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
         """
diff --git a/testing/data/vlmop.npz b/testing/data/vlmop.npz
index e92eb85..7e360cf 100644
Binary files a/testing/data/vlmop.npz and b/testing/data/vlmop.npz differ
diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py
index 4a2ebf3..999f2ae 100644
--- a/testing/test_acquisition.py
+++ b/testing/test_acquisition.py
@@ -3,90 +3,57 @@
 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)])
 
-def parabola2d(X):
-    return np.atleast_2d(np.sum(X ** 2, axis=1)).T
 
-def plane(X):
-    return X[:, [0]] - 0.5
+class SimpleAcquisition(GPflowOpt.acquisition.Acquisition):
+    def __init__(self, model):
+        super(SimpleAcquisition, self).__init__(model)
+        self.counter = 0
 
+    def _setup(self):
+        super(SimpleAcquisition, self)._setup()
+        self.counter += 1
 
-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 build_acquisition(self, Xcand):
+        return self.models[0].build_predict(Xcand)[0]
 
 
-
-class _TestAcquisition(object):
-    """
-    Defines some basic verifications for all acquisition functions. Test classes can derive from this
-    """
+class TestAcquisition(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)])
-
-    def load_data(self, file):
-        path = os.path.dirname(os.path.realpath(__file__))
-        return np.load(os.path.join(path, 'data', file))
-
-    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(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
+        self.model = create_parabola_model(domain)
+        self.acquisition = SimpleAcquisition(self.model)
 
-    def test_result_shape(self):
-        # Verify the returned shape of evaluate
-        design = GPflowOpt.design.RandomDesign(50, self.domain)
+    def run_setup(self):
+        # Optimize models & perform acquisition setup call.
+        self.acquisition._optimize_models()
+        self.acquisition._setup()
 
-        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)
-
-        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 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(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, 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, domain).generate())
+        self.assertFalse(self.acquisition._needs_setup)
+        self.assertEqual(self.acquisition.counter, 2)
 
     def test_data(self):
         # Test the data property
@@ -100,280 +67,167 @@ def test_data(self):
                                 msg="data property should return Tensors")
 
     def test_data_update(self):
-        # Verify a data update
-        design = GPflowOpt.design.RandomDesign(10, self.domain)
+        # Verify the effect of setting the data
+        design = GPflowOpt.design.RandomDesign(10, 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")
+        # 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.enable_scaling(self.domain)
+            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(domain)
         self.assertTrue(
-            all(m.wrapped.X.value in GPflowOpt.domain.UnitCube(self.domain.size) for m in self.acquisition.models))
-
-
-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")
-
+            all(m.wrapped.X.value in GPflowOpt.domain.UnitCube(domain.size) for m in self.acquisition.models))
+        self.assertTrue(self.acquisition._needs_setup)
 
-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_result_shape_tf(self):
+        # Verify the returned shape of evaluate
+        design = GPflowOpt.design.RandomDesign(50, domain)
 
-    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)
+        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_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)
+    def test_result_shape_np(self):
+        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, 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()))
 
-class TestHVProbabilityOfImprovement(_TestAcquisition, unittest.TestCase):
+        self.acquisition.optimize_restarts = 1
+        self.acquisition._optimize_models()
+        self.assertFalse(np.allclose(state, self.acquisition.get_free_state()))
 
-    @_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)
+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)
+)
 
-    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")
 
-    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(unittest.TestCase):
 
+    _multiprocess_can_split_ = True
 
-class _TestAcquisitionAggregation(_TestAcquisition):
-    def test_object_integrity(self):
-        for oper in self.acquisition.operands:
+    @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.assertListEqual(self.acquisition.models.sorted_params, self.models)
+        self.assertTrue(all(isinstance(m, GPflowOpt.models.ModelWrapper) for m in acquisition.models.sorted_params))
 
-    def test_data(self):
-        super(_TestAcquisitionAggregation, self).test_data()
-        np.testing.assert_allclose(self.acquisition.data[0], self.acquisition[0].data[0],
+    @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(self.acquisition.data[0], self.acquisition[1].data[0],
+        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, 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()
+        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 = self.acquisition.evaluate(design.generate())
+        p1 = 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")
+        np.testing.assert_allclose(p2, p1 / 2, rtol=1e-3)
 
-    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()
+    @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 = self.acquisition.evaluate(design.generate())
+        p1 = 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")
+        np.testing.assert_allclose(p2, np.sqrt(p1), rtol=1e-3)
 
-    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))
+    @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_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")
+    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))
 
-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:]]):
+    @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))
 
-    def test_marginalized_score(self):
+    @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 = self.acquisition.operands[0].evaluate(Xt)
-        ei_mcmc = self.acquisition.evaluate(Xt)
+        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
 
-    @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 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
@@ -383,18 +237,23 @@ 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")
 
     def test_hierarchy(self):
-        design = GPflowOpt.design.LatinHyperCube(16, self.domain)
+        design = GPflowOpt.design.LatinHyperCube(16, domain)
         X = design.generate()
         Yc = plane(X)
-        m1 = self.create_parabola_model()
-        m2 = self.create_parabola_model()
+        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)
@@ -404,8 +263,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(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])
@@ -428,7 +287,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
diff --git a/testing/test_datascaler.py b/testing/test_datascaler.py
index 07ffb9c..51f6670 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
new file mode 100644
index 0000000..f907d25
--- /dev/null
+++ b/testing/test_implementations.py
@@ -0,0 +1,186 @@
+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)
+        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):
+
+    _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)
+        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")
+
+    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)
+        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):
+
+    _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)
+        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):
+
+    _multiprocess_can_split_ = True
+
+    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")
+
+    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/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 28f9fbd..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))
 
@@ -253,6 +235,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))
@@ -261,15 +244,22 @@ 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):
         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)
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