Callback

GPflowOpt/acquisition/acquisition.py

@@ -168,7 +168,7 @@ def models(self):
 
         :return: list of GPflow models 
         """
-        return self._models
+        return self._models.sorted_params
 
     @property
     def data(self):
@@ -296,7 +296,7 @@ def _optimize_models(self):
 
     @Acquisition.models.getter
     def models(self):
-        return ParamList([model for acq in self.operands for model in acq.models.sorted_params])
+        return [model for acq in self.operands for model in acq.models]
 
     def enable_scaling(self, domain):
         for oper in self.operands:
@@ -375,27 +375,31 @@ class MCMCAcquistion(AcquisitionSum):
     """
     Apply MCMC over the hyperparameters of an acquisition function (= over the hyperparameters of the contained models).
 
-    The models passed into an object of this class are optimized with MLE, and then further sampled with HMC.
-    These hyperparameter samples are then set in copies of the acquisition.
+    The models passed into an object of this class are optimized with MLE (fast burn-in), and then further sampled with
+    HMC. These hyperparameter samples are then set in copies of the acquisition.
 
     For evaluating the underlying acquisition function, the predictions of the acquisition copies are averaged.
     """
     def __init__(self, acquisition, n_slices, **kwargs):
         assert isinstance(acquisition, Acquisition)
         assert n_slices > 0
-
-        copies = [copy.deepcopy(acquisition) for _ in range(n_slices - 1)]
-        for c in copies:
-            c.optimize_restarts = 0
-
         # the call to the constructor of the parent classes, will optimize acquisition, so it obtains the MLE solution.
-        super(MCMCAcquistion, self).__init__([acquisition] + copies)
+        super(MCMCAcquistion, self).__init__([acquisition]*n_slices)
+        self._needs_new_copies = True
         self._sample_opt = kwargs
 
     def _optimize_models(self):
         # Optimize model #1
         self.operands[0]._optimize_models()
 
+        # Copy it again if needed due to changed free state
+        if self._needs_new_copies:
+            new_copies = [copy.deepcopy(self.operands[0]) for _ in range(len(self.operands) - 1)]
+            for c in new_copies:
+                c.optimize_restarts = 0
+            self.operands = ParamList([self.operands[0]] + new_copies)
+            self._needs_new_copies = False
+
         # 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])
@@ -419,3 +423,11 @@ def set_data(self, X, Y):
     def build_acquisition(self, Xcand):
         # Average the predictions of the copies.
         return 1. / len(self.operands) * super(MCMCAcquistion, self).build_acquisition(Xcand)
+
+    def _kill_autoflow(self):
+        """
+        Following the recompilation of models, the free state might have changed. This means updating the samples can
+        cause inconsistencies and errors. Flag for recreation on next optimize
+        """
+        super(MCMCAcquistion, self)._kill_autoflow()
+        self._needs_new_copies = True

GPflowOpt/bo.py

@@ -18,9 +18,9 @@
 from scipy.optimize import OptimizeResult
 
 from .acquisition import Acquisition, MCMCAcquistion
-from .optim import Optimizer, SciPyOptimizer
-from .objective import ObjectiveWrapper
 from .design import Design, EmptyDesign
+from .objective import ObjectiveWrapper
+from .optim import Optimizer, SciPyOptimizer
 from .pareto import non_dominated_sort
 
 
@@ -32,7 +32,8 @@ class BayesianOptimizer(Optimizer):
     Additionally, it is configured with a separate optimizer for the acquisition function.
     """
 
-    def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=True, hyper_draws=None):
+    def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=True, hyper_draws=None,
+                 callback=None):
         """
         :param Domain domain: The optimization space.
         :param Acquisition acquisition: The acquisition function to optimize over the domain.
@@ -51,6 +52,12 @@ def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=Tr
             are obtained using Hamiltonian MC.
             (see `GPflow documentation <https://gpflow.readthedocs.io/en/latest//>`_ for details) for each model.
             The acquisition score is computed for each draw, and averaged.
+        :param callable callback: (optional) this function or object will be called after each evaluate, after the
+            data of all models has been updated with all models as retrieved by acquisition.models as argument without
+            the wrapping model handling any scaling . This allows custom model optimization strategies to be implemented.
+            All manipulations of GPflow models are permitted. Combined with the optimize_restarts parameter of
+            :class:`~.Acquisition` this allows several scenarios: do the optimization manually from the callback
+            (optimize_restarts equals zero),  orchoose the starting point + some random restarts (optimize_restarts > 0).
         """
         assert isinstance(acquisition, Acquisition)
         assert hyper_draws is None or hyper_draws > 0
@@ -69,6 +76,8 @@ def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=Tr
         initial = initial or EmptyDesign(domain)
         self.set_initial(initial.generate())
 
+        self._iter_callback = callback
+
     @Optimizer.domain.setter
     def domain(self, dom):
         assert (self.domain.size == dom.size)
@@ -86,6 +95,8 @@ def _update_model_data(self, newX, newY):
         assert self.acquisition.data[0].shape[1] == newX.shape[-1]
         assert self.acquisition.data[1].shape[1] == newY.shape[-1]
         assert newX.shape[0] == newY.shape[0]
+        if newX.size == 0:
+            return
         X = np.vstack((self.acquisition.data[0], newX))
         Y = np.vstack((self.acquisition.data[1], newY))
         self.acquisition.set_data(X, Y)
@@ -175,8 +186,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()
         values = fx(initial)
@@ -190,6 +200,10 @@ def inverse_acquisition(x):
 
         # Optimization loop
         for i in range(n_iter):
+            # If callback specified, and acquisition has the setup flag enabled (indicating an upcoming compilation,
+            # run the callback.
+            if self._iter_callback and self.acquisition._needs_setup:
+                self._iter_callback([m.wrapped for m in self.acquisition.models])
             result = self.optimizer.optimize(inverse_acquisition)
             self._update_model_data(result.x, fx(result.x))
 

GPflowOpt/models.py

@@ -15,6 +15,35 @@
 from GPflow.model import Model
 
 
+class ParentHook(object):
+    """
+    Temporary solution for fixing the recompilation issues (#37, GPflow issue #442).
+
+    An object of this class is returned when highest_parent is called on a model, which holds references to the highest
+    parentable, as well as the highest model class. When setting the needs recompile flag, this is intercepted and
+    performed on the model. At the same time, kill autoflow is called on the highest parent.
+    """
+    def __init__(self, highest_parent, highest_model):
+        self._hp = highest_parent
+        self._hm = highest_model
+
+    def __getattr__(self, item):
+        if item is '_needs_recompile':
+            return getattr(self._hm, item)
+        return getattr(self._hp, item)
+
+    def __setattr__(self, key, value):
+        if key in ['_hp', '_hm']:
+            object.__setattr__(self, key, value)
+            return
+        if key is '_needs_recompile':
+            setattr(self._hm, key, value)
+            if value:
+                self._hp._kill_autoflow()
+        else:
+            setattr(self._hp, key, value)
+
+
 class ModelWrapper(Parameterized):
     """
     Class for fast implementation of a wrapper for models defined in GPflow.
@@ -25,6 +54,7 @@ class ModelWrapper(Parameterized):
     AutoFlow methods are influenced following this pattern, the original AF storage (if existing) is unaffected and a
     new storage is added to the subclass.
     """
+
     def __init__(self, model):
         """
         :param model: model to be wrapped
@@ -45,6 +75,7 @@ def __getattr__(self, item):
             method = item[1:].rstrip('_AF_storage')
             if method in dir(self):
                 raise AttributeError("{0} has no attribute {1}".format(self.__class__.__name__, item))
+
         return getattr(self.wrapped, item)
 
     def __setattr__(self, key, value):
@@ -90,3 +121,11 @@ def __eq__(self, other):
     def name(self):
         name = super(ModelWrapper, self).name
         return ".".join([name, str.lower(self.__class__.__name__)])
+
+    @Parameterized.highest_parent.getter
+    def highest_parent(self):
+        """
+        Returns an instance of the ParentHook instead of the usual reference to a Parentable.
+        """
+        original_hp = super(ModelWrapper, self).highest_parent
+        return original_hp if isinstance(original_hp, ParentHook) else ParentHook(original_hp, self)

GPflowOpt/scaling.py

@@ -13,7 +13,6 @@
 # limitations under the License.
 
 from GPflow.param import DataHolder, AutoFlow
-from GPflow.model import GPModel
 from GPflow import settings
 import numpy as np
 from .transforms import LinearTransform, DataTransform
@@ -53,6 +52,7 @@ class DataScaler(ModelWrapper):
     required, it is the responsibility of the implementation to rescale the hyperparameters. Additionally, applying
     hyperpriors should anticipate for the scaled data.
     """
+
     def __init__(self, model, domain=None, normalize_Y=False):
         """
         :param model: model to be wrapped
@@ -90,7 +90,7 @@ def input_transform(self, t):
 
         :param t: :class:`.DataTransform` object: the new input transform.
         """
-        assert(isinstance(t, DataTransform))
+        assert (isinstance(t, DataTransform))
         X = self.X.value
         self._input_transform.assign(t)
         self.X = X

testing/test_acquisition.py

@@ -146,7 +146,8 @@ 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, GPflowOpt.models.ModelWrapper) for m in acquisition.models.sorted_params))
+            self.assertTrue(all(isinstance(m, GPflowOpt.models.ModelWrapper) for m in acquisition.models))
+
 
     @parameterized.expand(list(zip(aggregations)))
     def test_data(self, acquisition):
@@ -217,10 +218,23 @@ def test_marginalized_score(self, acquisition):
         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)
-
+    def test_mcmc_acq(self):
+        acquisition = GPflowOpt.acquisition.MCMCAcquistion(
+            GPflowOpt.acquisition.ExpectedImprovement(create_parabola_model(domain)), 10)
+        for oper in acquisition.operands:
+            self.assertListEqual(acquisition.models, oper.models)
+            self.assertEqual(acquisition.operands[0], oper)
+        self.assertTrue(acquisition._needs_new_copies)
+        acquisition._optimize_models()
+        self.assertListEqual(acquisition.models, acquisition.operands[0].models)
+        for oper in acquisition.operands[1:]:
+            self.assertNotEqual(acquisition.operands[0], oper)
+        self.assertFalse(acquisition._needs_new_copies)
+        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)
 
 class TestJointAcquisition(unittest.TestCase):
 
@@ -298,3 +312,26 @@ def test_multi_aggr(self):
         joint = first * second
         self.assertIsInstance(joint, GPflowOpt.acquisition.AcquisitionProduct)
         self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3, acq4])
+
+
+class TestRecompile(unittest.TestCase):
+    """
+    Regression test for #37
+    """
+    def test_vgp(self):
+        domain = GPflowOpt.domain.UnitCube(2)
+        X = GPflowOpt.design.RandomDesign(10, domain).generate()
+        Y = np.sin(X[:,[0]])
+        m = GPflow.vgp.VGP(X, Y, GPflow.kernels.RBF(2), GPflow.likelihoods.Gaussian())
+        m._compile()
+        acq = GPflowOpt.acquisition.ExpectedImprovement(m)
+        self.assertFalse(m._needs_recompile)
+        acq.evaluate(GPflowOpt.design.RandomDesign(10, domain).generate())
+        self.assertTrue(hasattr(acq, '_evaluate_AF_storage'))
+
+        Xnew = GPflowOpt.design.RandomDesign(5, domain).generate()
+        Ynew = np.sin(Xnew[:,[0]])
+        acq.set_data(np.vstack((X, Xnew)), np.vstack((Y, Ynew)))
+        self.assertFalse(hasattr(acq, '_needs_recompile'))
+        self.assertFalse(hasattr(acq, '_evaluate_AF_storage'))
+        acq.evaluate(GPflowOpt.design.RandomDesign(10, domain).generate())

testing/test_datascaler.py

@@ -106,3 +106,4 @@ def test_predict_scaling(self):
         fs = n.predict_density(Xt, Yt)
         np.testing.assert_allclose(fr, fs, rtol=1e-2)
 
+

testing/test_modelwrapper.py

@@ -111,6 +111,31 @@ def test_name(self):
         n = MethodOverride(create_parabola_model(GPflowOpt.domain.UnitCube(2)))
         self.assertEqual(n.name, 'unnamed.methodoverride')
 
+    def test_parent_hook(self):
+        self.m.optimize(maxiter=5)
+        w = GPflowOpt.models.ModelWrapper(self.m)
+        self.assertTrue(isinstance(self.m.highest_parent, GPflowOpt.models.ParentHook))
+        self.assertEqual(self.m.highest_parent._hp, w)
+        self.assertEqual(self.m.highest_parent._hm, w)
+
+        w2 = GPflowOpt.models.ModelWrapper(w)
+        self.assertEqual(self.m.highest_parent._hp, w2)
+        self.assertEqual(self.m.highest_parent._hm, w2)
+
+        p = GPflow.param.Parameterized()
+        p.model = w2
+        self.assertEqual(self.m.highest_parent._hp, p)
+        self.assertEqual(self.m.highest_parent._hm, w2)
 
+        p.predictor = create_parabola_model(GPflowOpt.domain.UnitCube(2))
+        p.predictor.predict_f(p.predictor.X.value)
+        self.assertTrue(hasattr(p.predictor, '_predict_f_AF_storage'))
+        self.assertFalse(self.m._needs_recompile)
+        self.m.highest_parent._needs_recompile = True
+        self.assertFalse('_needs_recompile' in p.__dict__)
+        self.assertFalse('_needs_recompile' in w.__dict__)
+        self.assertFalse('_needs_recompile' in w2.__dict__)
+        self.assertTrue(self.m._needs_recompile)
+        self.assertFalse(hasattr(p.predictor, '_predict_f_AF_storage'))