From ac559f37f65569a131a075a50260816cb22c92ab Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Sun, 30 Jul 2017 23:11:24 +0200 Subject: [PATCH 01/15] Adding basic support for batch optimisation to the framework --- GPflowOpt/acquisition/acquisition.py | 17 +++++++++-------- GPflowOpt/bo.py | 16 +++++++++++++--- 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py index b386cf8..6939069 100644 --- a/GPflowOpt/acquisition/acquisition.py +++ b/GPflowOpt/acquisition/acquisition.py @@ -43,7 +43,7 @@ class Acquisition(Parameterized): (for instance for constrained optimization) """ - def __init__(self, models=[], optimize_restarts=5): + def __init__(self, models=[], optimize_restarts=5, batch_size=1): super(Acquisition, self).__init__() self._models = ParamList([DataScaler(m) for m in np.atleast_1d(models).tolist()]) self._default_params = list(map(lambda m: m.get_free_state(), self._models)) @@ -51,6 +51,7 @@ def __init__(self, models=[], optimize_restarts=5): assert (optimize_restarts >= 0) self._optimize_restarts = optimize_restarts self._optimize_models() + self.batch_size = batch_size def _optimize_models(self): """ @@ -82,7 +83,7 @@ def _optimize_models(self): best_idx = np.argmin([r.fun for r in runs]) model.set_state(runs[best_idx].x) - def build_acquisition(self): + def build_acquisition(self, *args): raise NotImplementedError def enable_scaling(self, domain): @@ -180,7 +181,7 @@ def evaluate_with_gradients(self, Xcand): """ AutoFlow method to compute the acquisition scores for candidates, also returns the gradients. """ - acq = self.build_acquisition(Xcand) + acq = self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size)) return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] @AutoFlow((float_type, [None, None])) @@ -188,7 +189,7 @@ def evaluate(self, Xcand): """ AutoFlow method to compute the acquisition scores for candidates, without returning the gradients. """ - return self.build_acquisition(Xcand) + return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size)) def __add__(self, other): """ @@ -261,8 +262,8 @@ def constraint_indices(self): def feasible_data_index(self): return np.all(np.vstack(map(lambda o: o.feasible_data_index(), self.operands)), axis=0) - def build_acquisition(self, Xcand): - return self._oper(tf.concat(list(map(lambda operand: operand.build_acquisition(Xcand), self.operands)), 1), + def build_acquisition(self, *args): + return self._oper(tf.concat(list(map(lambda operand: operand.build_acquisition(*args), self.operands)), 1), axis=1, keep_dims=True, name=self.__class__.__name__) def __getitem__(self, item): @@ -344,6 +345,6 @@ def set_data(self, X, Y): self.setup() return offset - def build_acquisition(self, Xcand): + def build_acquisition(self, *args): # Average the predictions of the copies. - return 1. / len(self.operands) * super(MCMCAcquistion, self).build_acquisition(Xcand) + return 1. / len(self.operands) * super(MCMCAcquistion, self).build_acquisition(*args) diff --git a/GPflowOpt/bo.py b/GPflowOpt/bo.py index bf506d2..57a7432 100644 --- a/GPflowOpt/bo.py +++ b/GPflowOpt/bo.py @@ -13,6 +13,7 @@ # limitations under the License. from contextlib import contextmanager +import copy import numpy as np from scipy.optimize import OptimizeResult @@ -54,12 +55,20 @@ 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) + batch_domain = np.sum([copy.deepcopy(domain) for i in range(acquisition.batch_size)]) + + # Configure scaling if scaling: acquisition.enable_scaling(domain) + + # Configure MCMC self.acquisition = acquisition if hyper_draws is None else MCMCAcquistion(acquisition, hyper_draws) - self.optimizer = optimizer or SciPyOptimizer(domain) - self.optimizer.domain = domain + # Setup optimizer + self.optimizer = optimizer or SciPyOptimizer(batch_domain) + self.optimizer.domain = batch_domain + + # Setup initial evaluations initial = initial or EmptyDesign(domain) self.set_initial(initial.generate()) @@ -165,7 +174,8 @@ def inverse_acquisition(x): # Optimization loop for i in range(n_iter): result = self.optimizer.optimize(inverse_acquisition) - self._update_model_data(result.x, fx(result.x)) + Xnew = np.vstack(np.split(result.x, self.acquisition.batch_size)) + self._update_model_data(Xnew, fx(Xnew)) return self._create_bo_result(True, "OK") From 82642ddf01f1bd2e01a08349d44fea667ec4e724 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Wed, 23 Aug 2017 21:30:50 +0200 Subject: [PATCH 02/15] Add missing axis parameters & introduced a batch function in domain --- GPflowOpt/acquisition/acquisition.py | 4 ++-- GPflowOpt/bo.py | 4 ++-- GPflowOpt/domain.py | 4 ++++ testing/test_domain.py | 7 +++++++ 4 files changed, 15 insertions(+), 4 deletions(-) diff --git a/GPflowOpt/acquisition/acquisition.py b/GPflowOpt/acquisition/acquisition.py index a32a0b5..a3c3a41 100644 --- a/GPflowOpt/acquisition/acquisition.py +++ b/GPflowOpt/acquisition/acquisition.py @@ -204,7 +204,7 @@ def evaluate_with_gradients(self, Xcand): :return: acquisition scores, size N x 1 the gradients of the acquisition scores, size N x D """ - acq = self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size)) + acq = self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] @AutoFlow((float_type, [None, None])) @@ -214,7 +214,7 @@ def evaluate(self, Xcand): :return: acquisition scores, size N x 1 """ - return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size)) + return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) def __add__(self, other): """ diff --git a/GPflowOpt/bo.py b/GPflowOpt/bo.py index 74fed16..b6cee77 100644 --- a/GPflowOpt/bo.py +++ b/GPflowOpt/bo.py @@ -58,7 +58,7 @@ 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) - batch_domain = np.sum([copy.deepcopy(domain) for i in range(acquisition.batch_size)]) + batch_domain = domain.batch(acquisition.batch_size) # Configure scaling if scaling: @@ -190,7 +190,7 @@ def inverse_acquisition(x): # Optimization loop for i in range(n_iter): result = self.optimizer.optimize(inverse_acquisition) - Xnew = np.vstack(np.split(result.x, self.acquisition.batch_size)) + Xnew = np.vstack(np.split(result.x, self.acquisition.batch_size), axis=1) self._update_model_data(Xnew, fx(Xnew)) return self._create_bo_result(True, "OK") diff --git a/GPflowOpt/domain.py b/GPflowOpt/domain.py index f95275c..94e4855 100644 --- a/GPflowOpt/domain.py +++ b/GPflowOpt/domain.py @@ -13,6 +13,7 @@ # limitations under the License. import numpy as np +import copy from itertools import chain from GPflow.param import Parentable @@ -53,6 +54,9 @@ def size(self): """ return sum(map(lambda param: param.size, self._parameters)) + def batch(self, size): + return np.sum([copy.deepcopy(self) for i in range(size)]) + def __setattr__(self, key, value): super(Domain, self).__setattr__(key, value) if key is not '_parent': diff --git a/testing/test_domain.py b/testing/test_domain.py index 87c0c9c..fd3bd9a 100644 --- a/testing/test_domain.py +++ b/testing/test_domain.py @@ -70,6 +70,13 @@ def test_value(self): self.assertTupleEqual(p.value.shape, (1,), msg="Default value has incorrect shape.") self.assertTrue(np.allclose(p.value, 0.2), msg="Parameter has incorrect initialized value") + def test_batch(self): + p = GPflowOpt.domain.ContinuousParameter("x1", 0, 1) + b = p.batch(3) + self.assertEqual(b.size, 3) + self.assertTrue(np.allclose(b.lower, 0)) + self.assertTrue(np.allclose(b.upper, 1)) + class TestHypercubeDomain(unittest.TestCase): From 76b14e8fcf3e8e9e1ba4992838e1e251712b4756 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Wed, 23 Aug 2017 21:53:56 +0200 Subject: [PATCH 03/15] Put axis parameter in the wrong function --- GPflowOpt/bo.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPflowOpt/bo.py b/GPflowOpt/bo.py index b6cee77..2a4be9f 100644 --- a/GPflowOpt/bo.py +++ b/GPflowOpt/bo.py @@ -190,7 +190,7 @@ def inverse_acquisition(x): # Optimization loop for i in range(n_iter): result = self.optimizer.optimize(inverse_acquisition) - Xnew = np.vstack(np.split(result.x, self.acquisition.batch_size), axis=1) + Xnew = np.vstack(np.split(result.x, self.acquisition.batch_size, axis=1)) self._update_model_data(Xnew, fx(Xnew)) return self._create_bo_result(True, "OK") From 418e716d3ff24b92c2759c46f3bfc741ed64fd74 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Thu, 7 Sep 2017 13:48:31 +0200 Subject: [PATCH 04/15] Implementation of ParallelBatchAcquisition --- gpflowopt/acquisition/__init__.py | 5 ++- gpflowopt/acquisition/acquisition.py | 43 ++++++++++++++++++++--- gpflowopt/bo.py | 14 +++----- testing/test_acquisition.py | 52 ++++++++++++++++++++++++---- 4 files changed, 93 insertions(+), 21 deletions(-) diff --git a/gpflowopt/acquisition/__init__.py b/gpflowopt/acquisition/__init__.py index eca69fe..e75698f 100644 --- a/gpflowopt/acquisition/__init__.py +++ b/gpflowopt/acquisition/__init__.py @@ -13,7 +13,8 @@ # limitations under the License. # Framework components and interfaces -from .acquisition import Acquisition, AcquisitionAggregation, AcquisitionProduct, AcquisitionSum, MCMCAcquistion +from .acquisition import (Acquisition, ParallelBatchAcquisition, AcquisitionAggregation, AcquisitionProduct, \ + AcquisitionSum, MCMCAcquistion) # Single objective from .ei import ExpectedImprovement @@ -23,5 +24,7 @@ # Multiobjective from .hvpoi import HVProbabilityOfImprovement +# Batch + # Black-box constraint from .pof import ProbabilityOfFeasibility diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index 07e1968..81693d4 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -75,7 +75,7 @@ class Acquisition(Parameterized): objectives. """ - def __init__(self, models=[], optimize_restarts=5, batch_size=1): + def __init__(self, models=[], optimize_restarts=5): """ :param models: list of GPflow models representing our beliefs about the problem :param optimize_restarts: number of optimization restarts to use when training the models @@ -87,7 +87,6 @@ def __init__(self, models=[], optimize_restarts=5, batch_size=1): assert (optimize_restarts >= 0) self.optimize_restarts = optimize_restarts - self.batch_size = batch_size self._needs_setup = True def _optimize_models(self): @@ -122,6 +121,13 @@ def _optimize_models(self): best_idx = np.argmin([r.fun for r in runs]) model.set_state(runs[best_idx].x) + def get_suggestion(self, optimizer): + def inverse_acquisition(x): + return tuple(map(lambda r: -r, self.evaluate_with_gradients(np.atleast_2d(x)))) + + result = optimizer.optimize(inverse_acquisition) + return result.x + def build_acquisition(self, *args): raise NotImplementedError @@ -253,7 +259,7 @@ def evaluate_with_gradients(self, Xcand): :return: acquisition scores, size N x 1 the gradients of the acquisition scores, size N x D """ - acq = self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + acq = self.build_acquisition(Xcand) return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] @setup_required @@ -264,7 +270,7 @@ def evaluate(self, Xcand): :return: acquisition scores, size N x 1 """ - return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + return self.build_acquisition(Xcand) def __add__(self, other): """ @@ -298,6 +304,35 @@ def __setattr__(self, key, value): self.highest_parent._needs_setup = True +class ParallelBatchAcquisition(Acquisition): + + def __init__(self, models=[], optimize_restarts=5, batch_size=1): + super(ParallelBatchAcquisition, self).__init__(models, optimize_restarts=optimize_restarts) + self.batch_size = batch_size + + def get_suggestion(self, optimizer): + opt = copy.deepcopy(optimizer) + batch_domain = optimizer.domain.batch(self.batch_size) + opt.domain = batch_domain + + def inverse_acquisition(x): + return tuple(map(lambda r: -r, self.evaluate_with_gradients(np.atleast_2d(x)))) + + result = opt.optimize(inverse_acquisition) + return np.vstack(np.split(result.x, self.batch_size, axis=1)) + + @setup_required + @AutoFlow((float_type, [None, None])) + def evaluate_with_gradients(self, Xcand): + acq = self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] + + @setup_required + @AutoFlow((float_type, [None, None])) + def evaluate(self, Xcand): + return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + + class AcquisitionAggregation(Acquisition): """ Aggregates multiple acquisition functions, using a TensorFlow reduce operation. diff --git a/gpflowopt/bo.py b/gpflowopt/bo.py index 685bf99..069aca6 100644 --- a/gpflowopt/bo.py +++ b/gpflowopt/bo.py @@ -94,8 +94,6 @@ 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) - batch_domain = domain.batch(acquisition.batch_size) - # Configure MCMC self._scaling = scaling if self._scaling: @@ -104,8 +102,8 @@ def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=Tr self.acquisition = acquisition if hyper_draws is None else MCMCAcquistion(acquisition, hyper_draws) # Setup optimizer - self.optimizer = optimizer or SciPyOptimizer(batch_domain) - self.optimizer.domain = batch_domain + self.optimizer = optimizer or SciPyOptimizer(domain) + self.optimizer.domain = domain # Setup initial evaluations initial = initial or EmptyDesign(domain) @@ -230,17 +228,13 @@ def _optimize(self, fx, n_iter): # Remove initial design for additional calls to optimize to proceed optimization self.set_initial(EmptyDesign(self.domain).generate()) - def inverse_acquisition(x): - return tuple(map(lambda r: -r, self.acquisition.evaluate_with_gradients(np.atleast_2d(x)))) - # Optimization loop for i in range(n_iter): # If a callback is specified, and acquisition has the setup flag enabled (indicating an upcoming - # compilation), run the callback. + # setup), run the callback. if self._model_callback and self.acquisition._needs_setup: self._model_callback([m.wrapped for m in self.acquisition.models]) - result = self.optimizer.optimize(inverse_acquisition) - Xnew = np.vstack(np.split(result.x, self.acquisition.batch_size, axis=1)) + Xnew = self.acquisition.get_suggestion(self.optimizer) self._update_model_data(Xnew, fx(Xnew)) return self._create_bo_result(True, "OK") diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index cc81e6c..717c9e4 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -19,7 +19,7 @@ def _setup(self): self.counter += 1 def build_acquisition(self, Xcand): - return self.models[0].build_predict(Xcand)[0] + return -self.models[0].build_predict(Xcand)[0] class TestAcquisition(unittest.TestCase): @@ -30,11 +30,6 @@ def setUp(self): self.model = create_parabola_model(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.") @@ -122,6 +117,48 @@ def test_optimize(self): self.acquisition._optimize_models() self.assertFalse(np.allclose(state, self.acquisition.get_free_state())) + def test_get_suggestion(self): + opt = gpflowopt.optim.SciPyOptimizer(domain) + Xnew = self.acquisition.get_suggestion(opt) + self.assertTupleEqual(Xnew.shape, (1, 2)) + self.assertTrue(np.allclose(Xnew, 0)) + + +class SimpleParallelBatch(gpflowopt.acquisition.ParallelBatchAcquisition): + def __init__(self, model): + super(SimpleParallelBatch, self).__init__(model, batch_size=2) + self.n_args = 0 + + def build_acquisition(self, *args): + self.n_args = len(args) + return -tf.add(self.models[0].build_predict(args[0])[0], self.models[0].build_predict(args[1]+1)[0]) + + +class TestParallelBatchAcquisition(unittest.TestCase): + + _multiprocess_can_split_ = True + + def setUp(self): + self.model = create_parabola_model(domain) + self.acquisition = SimpleParallelBatch(self.model) + + def test_arg_splitting(self): + self.acquisition.evaluate(np.random.rand(10, 4)) + self.assertEqual(self.acquisition.n_args, 2) + + self.acquisition.n_args = 0 + X = np.random.rand(10, 4) + A = self.acquisition.evaluate_with_gradients(X) + self.assertEqual(self.acquisition.n_args, 2) + + def test_get_suggestion(self): + opt = gpflowopt.optim.StagedOptimizer([gpflowopt.optim.MCOptimizer(domain, 200), + gpflowopt.optim.SciPyOptimizer(domain)]) + Xnew = self.acquisition.get_suggestion(opt) + self.assertTupleEqual(Xnew.shape, (2, 2)) + self.assertTrue(np.allclose(Xnew[0, :], 0, atol=1e-4)) + self.assertTrue(np.allclose(Xnew[1, :], -1, atol=1e-4)) + aggregations = list() aggregations.append(gpflowopt.acquisition.AcquisitionSum([ @@ -335,3 +372,6 @@ def test_vgp(self): self.assertFalse(hasattr(acq, '_needs_recompile')) self.assertFalse(hasattr(acq, '_evaluate_AF_storage')) acq.evaluate(gpflowopt.design.RandomDesign(10, domain).generate()) + + + From a0b6d0ee836bb373324fc5536b77c2bf84588669 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Thu, 7 Sep 2017 14:25:17 +0200 Subject: [PATCH 05/15] Making inverse acquisition available as a private method --- gpflowopt/acquisition/acquisition.py | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index 81693d4..e285e35 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -121,11 +121,11 @@ def _optimize_models(self): best_idx = np.argmin([r.fun for r in runs]) model.set_state(runs[best_idx].x) - def get_suggestion(self, optimizer): - def inverse_acquisition(x): - return tuple(map(lambda r: -r, self.evaluate_with_gradients(np.atleast_2d(x)))) + def _inverse_acquisition(self, x): + return tuple(map(lambda r: -r, self.evaluate_with_gradients(np.atleast_2d(x)))) - result = optimizer.optimize(inverse_acquisition) + def get_suggestion(self, optimizer): + result = optimizer.optimize(self._inverse_acquisition) return result.x def build_acquisition(self, *args): @@ -314,11 +314,7 @@ def get_suggestion(self, optimizer): opt = copy.deepcopy(optimizer) batch_domain = optimizer.domain.batch(self.batch_size) opt.domain = batch_domain - - def inverse_acquisition(x): - return tuple(map(lambda r: -r, self.evaluate_with_gradients(np.atleast_2d(x)))) - - result = opt.optimize(inverse_acquisition) + result = opt.optimize(self._inverse_acquisition) return np.vstack(np.split(result.x, self.batch_size, axis=1)) @setup_required From 226b9eb127f5d08365b8b86231c58da0d1cd39c0 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Thu, 7 Sep 2017 15:21:36 +0200 Subject: [PATCH 06/15] Minor interface fix --- gpflowopt/acquisition/acquisition.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index e285e35..fb95ea8 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -128,7 +128,7 @@ def get_suggestion(self, optimizer): result = optimizer.optimize(self._inverse_acquisition) return result.x - def build_acquisition(self, *args): + def build_acquisition(self, candidates): raise NotImplementedError def enable_scaling(self, domain): @@ -317,6 +317,9 @@ def get_suggestion(self, optimizer): result = opt.optimize(self._inverse_acquisition) return np.vstack(np.split(result.x, self.batch_size, axis=1)) + def build_acquisition(self, *args): + raise NotImplementedError + @setup_required @AutoFlow((float_type, [None, None])) def evaluate_with_gradients(self, Xcand): From c29c35b4d278b8b4c240444025c8a1364e185c04 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Thu, 7 Sep 2017 23:14:59 +0200 Subject: [PATCH 07/15] Following the re-design, ParallelBatch is now the parent class. --- gpflowopt/acquisition/__init__.py | 2 +- gpflowopt/acquisition/acquisition.py | 121 +++++++++++++-------------- testing/test_acquisition.py | 77 +++++++++-------- 3 files changed, 98 insertions(+), 102 deletions(-) diff --git a/gpflowopt/acquisition/__init__.py b/gpflowopt/acquisition/__init__.py index e75698f..4cda5f5 100644 --- a/gpflowopt/acquisition/__init__.py +++ b/gpflowopt/acquisition/__init__.py @@ -14,7 +14,7 @@ # Framework components and interfaces from .acquisition import (Acquisition, ParallelBatchAcquisition, AcquisitionAggregation, AcquisitionProduct, \ - AcquisitionSum, MCMCAcquistion) + AcquisitionSum, MCMCAcquistion, setup_required) # Single objective from .ei import ExpectedImprovement diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index fb95ea8..3e96b2a 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -36,7 +36,7 @@ def setup_required(method): """ @wraps(method) def runnable(instance, *args, **kwargs): - assert isinstance(instance, Acquisition) + assert isinstance(instance, ParallelBatchAcquisition) hp = instance.highest_parent if hp._needs_setup: # Avoid infinite loops, caused by setup() somehow invoking the evaluate on another acquisition @@ -54,33 +54,14 @@ def runnable(instance, *args, **kwargs): 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. 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. - """ +class ParallelBatchAcquisition(Parameterized): - def __init__(self, models=[], optimize_restarts=5): + def __init__(self, models=[], optimize_restarts=5, batch_size=1): """ :param models: list of GPflow models representing our beliefs about the problem :param optimize_restarts: number of optimization restarts to use when training the models """ - super(Acquisition, self).__init__() + super(ParallelBatchAcquisition, self).__init__() models = np.atleast_1d(models) assert all(isinstance(model, (Model, ModelWrapper)) for model in models) self._models = ParamList([DataScaler(m) for m in models]) @@ -88,6 +69,7 @@ def __init__(self, models=[], optimize_restarts=5): assert (optimize_restarts >= 0) self.optimize_restarts = optimize_restarts self._needs_setup = True + self.batch_size = batch_size def _optimize_models(self): """ @@ -125,16 +107,19 @@ def _inverse_acquisition(self, x): return tuple(map(lambda r: -r, self.evaluate_with_gradients(np.atleast_2d(x)))) def get_suggestion(self, optimizer): - result = optimizer.optimize(self._inverse_acquisition) - return result.x + opt = copy.deepcopy(optimizer) + batch_domain = optimizer.domain.batch(self.batch_size) + opt.domain = batch_domain + result = opt.optimize(self._inverse_acquisition) + return np.vstack(np.split(result.x, self.batch_size, axis=1)) - def build_acquisition(self, candidates): + def build_acquisition(self, *args): raise NotImplementedError 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. @@ -178,8 +163,8 @@ def set_data(self, X, Y): def models(self): """ The GPflow models representing our beliefs of the optimization problem. - - :return: list of GPflow models + + :return: list of GPflow models """ return self._models.sorted_params @@ -208,9 +193,9 @@ def constraint_indices(self): def objective_indices(self): """ Method returning the indices of the model outputs which are objective functions. - + By default all outputs are objectives. - + :return: indices to the objectives, size R """ return np.setdiff1d(np.arange(self.data[1].shape[1]), self.constraint_indices()) @@ -219,9 +204,9 @@ def feasible_data_index(self): """ Returns a boolean array indicating which data points are considered feasible (according to the acquisition function(s) ) and which not. - + By default all data is considered feasible. - + :return: logical indices to the feasible data points, size N """ return np.ones(self.data[0].shape[0], dtype=bool) @@ -229,7 +214,7 @@ def feasible_data_index(self): 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 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). @@ -250,6 +235,46 @@ def _setup_objectives(self): if self.constraint_indices().size == 0: self._setup() + @setup_required + @AutoFlow((float_type, [None, None])) + def evaluate_with_gradients(self, Xcand): + acq = self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] + + @setup_required + @AutoFlow((float_type, [None, None])) + def evaluate(self, Xcand): + return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + + +class Acquisition(ParallelBatchAcquisition): + """ + 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. 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 get_suggestion(self, optimizer): + result = optimizer.optimize(self._inverse_acquisition) + return result.x + + def build_acquisition(self, candidates): + raise NotImplementedError + @setup_required @AutoFlow((float_type, [None, None])) def evaluate_with_gradients(self, Xcand): @@ -304,34 +329,6 @@ def __setattr__(self, key, value): self.highest_parent._needs_setup = True -class ParallelBatchAcquisition(Acquisition): - - def __init__(self, models=[], optimize_restarts=5, batch_size=1): - super(ParallelBatchAcquisition, self).__init__(models, optimize_restarts=optimize_restarts) - self.batch_size = batch_size - - def get_suggestion(self, optimizer): - opt = copy.deepcopy(optimizer) - batch_domain = optimizer.domain.batch(self.batch_size) - opt.domain = batch_domain - result = opt.optimize(self._inverse_acquisition) - return np.vstack(np.split(result.x, self.batch_size, axis=1)) - - def build_acquisition(self, *args): - raise NotImplementedError - - @setup_required - @AutoFlow((float_type, [None, None])) - def evaluate_with_gradients(self, Xcand): - acq = self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) - return acq, tf.gradients(acq, [Xcand], name="acquisition_gradient")[0] - - @setup_required - @AutoFlow((float_type, [None, None])) - def evaluate(self, Xcand): - return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) - - class AcquisitionAggregation(Acquisition): """ Aggregates multiple acquisition functions, using a TensorFlow reduce operation. diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index 717c9e4..1997684 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -10,43 +10,53 @@ class SimpleAcquisition(gpflowopt.acquisition.Acquisition): + + def build_acquisition(self, Xcand): + return -self.models[0].build_predict(Xcand)[0] + + +class SimpleParallelBatch(gpflowopt.acquisition.ParallelBatchAcquisition): def __init__(self, model): - super(SimpleAcquisition, self).__init__(model) + super(SimpleParallelBatch, self).__init__(model, batch_size=2) + self.n_args = 0 self.counter = 0 def _setup(self): - super(SimpleAcquisition, self)._setup() + super(SimpleParallelBatch, self)._setup() self.counter += 1 - def build_acquisition(self, Xcand): - return -self.models[0].build_predict(Xcand)[0] + def build_acquisition(self, *args): + self.n_args = len(args) + return -tf.add(self.models[0].build_predict(args[0])[0], self.models[0].build_predict(args[1]+1)[0]) -class TestAcquisition(unittest.TestCase): +class TestParallelBatchAcquisition(unittest.TestCase): _multiprocess_can_split_ = True def setUp(self): self.model = create_parabola_model(domain) - self.acquisition = SimpleAcquisition(self.model) + self.acquisition = SimpleParallelBatch(self.model) def test_object_integrity(self): self.assertEqual(len(self.acquisition.models), 1, msg="Model list has incorrect length.") self.assertEqual(self.acquisition.models[0], self.model, msg="Incorrect model stored.") def test_setup_trigger(self): + Xrand = np.hstack((gpflowopt.design.RandomDesign(10, domain).generate(), + gpflowopt.design.RandomDesign(10, domain).generate())) 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.acquisition.evaluate(Xrand) 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.acquisition.evaluate_with_gradients(Xrand) self.assertFalse(self.acquisition._needs_setup) self.assertEqual(self.acquisition.counter, 2) @@ -87,25 +97,24 @@ def test_enable_scaling(self): def test_result_shape_tf(self): # Verify the returned shape of evaluate - design = gpflowopt.design.RandomDesign(50, domain) - 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) + tens = self.acquisition.build_acquisition(x_tf, x_tf) self.assertTrue(isinstance(tens, tf.Tensor), msg="no Tensor was returned") def test_result_shape_np(self): design = gpflowopt.design.RandomDesign(50, domain) - res = self.acquisition.evaluate(design.generate()) + candidates = np.hstack((design.generate(), design.generate())) + res = self.acquisition.evaluate(candidates) self.assertTupleEqual(res.shape, (50, 1)) - res = self.acquisition.evaluate_with_gradients(design.generate()) + res = self.acquisition.evaluate_with_gradients(candidates) 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)) + self.assertTupleEqual(res[1].shape, (50, domain.size * 2)) def test_optimize(self): self.acquisition.optimize_restarts = 0 @@ -117,31 +126,6 @@ def test_optimize(self): self.acquisition._optimize_models() self.assertFalse(np.allclose(state, self.acquisition.get_free_state())) - def test_get_suggestion(self): - opt = gpflowopt.optim.SciPyOptimizer(domain) - Xnew = self.acquisition.get_suggestion(opt) - self.assertTupleEqual(Xnew.shape, (1, 2)) - self.assertTrue(np.allclose(Xnew, 0)) - - -class SimpleParallelBatch(gpflowopt.acquisition.ParallelBatchAcquisition): - def __init__(self, model): - super(SimpleParallelBatch, self).__init__(model, batch_size=2) - self.n_args = 0 - - def build_acquisition(self, *args): - self.n_args = len(args) - return -tf.add(self.models[0].build_predict(args[0])[0], self.models[0].build_predict(args[1]+1)[0]) - - -class TestParallelBatchAcquisition(unittest.TestCase): - - _multiprocess_can_split_ = True - - def setUp(self): - self.model = create_parabola_model(domain) - self.acquisition = SimpleParallelBatch(self.model) - def test_arg_splitting(self): self.acquisition.evaluate(np.random.rand(10, 4)) self.assertEqual(self.acquisition.n_args, 2) @@ -160,6 +144,21 @@ def test_get_suggestion(self): self.assertTrue(np.allclose(Xnew[1, :], -1, atol=1e-4)) +class TestAcquisition(unittest.TestCase): + + _multiprocess_can_split_ = True + + def setUp(self): + self.model = create_parabola_model(domain) + self.acquisition = SimpleAcquisition(self.model) + + def test_get_suggestion(self): + opt = gpflowopt.optim.SciPyOptimizer(domain) + Xnew = self.acquisition.get_suggestion(opt) + self.assertTupleEqual(Xnew.shape, (1, 2)) + self.assertTrue(np.allclose(Xnew, 0)) + + aggregations = list() aggregations.append(gpflowopt.acquisition.AcquisitionSum([ gpflowopt.acquisition.ExpectedImprovement(create_parabola_model(domain)), From 6bd2561f7550b6073e5d6dd3df798e2775ec0602 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Sun, 17 Sep 2017 00:29:30 +0200 Subject: [PATCH 08/15] Major rework of acquisition framework. Introduced an interface class, created a wrapper class and changed hierarchy accordingly --- gpflowopt/acquisition/__init__.py | 2 +- gpflowopt/acquisition/acquisition.py | 488 ++++++++++++++++++--------- gpflowopt/bo.py | 5 +- testing/test_acquisition.py | 18 +- 4 files changed, 346 insertions(+), 167 deletions(-) diff --git a/gpflowopt/acquisition/__init__.py b/gpflowopt/acquisition/__init__.py index 4cda5f5..a8f4dba 100644 --- a/gpflowopt/acquisition/__init__.py +++ b/gpflowopt/acquisition/__init__.py @@ -14,7 +14,7 @@ # Framework components and interfaces from .acquisition import (Acquisition, ParallelBatchAcquisition, AcquisitionAggregation, AcquisitionProduct, \ - AcquisitionSum, MCMCAcquistion, setup_required) + AcquisitionSum, MCMCAcquistion, IAcquisition, setup_required) # Single objective from .ei import ExpectedImprovement diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index d476ff8..72b8698 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -25,6 +25,7 @@ import copy from functools import wraps +from abc import ABCMeta, abstractmethod, abstractproperty float_type = settings.dtypes.float_type @@ -54,23 +55,26 @@ def runnable(instance, *args, **kwargs): return runnable -class ParallelBatchAcquisition(Parameterized): +class IAcquisition(Parameterized): # pragma: no cover + """ + Interface for Acquisition functions mapping the belief represented by a Bayesian model into a score indicating how + promising a point is for evaluation. - def __init__(self, models=[], optimize_restarts=5, batch_size=1): - """ - :param models: list of GPflow models representing our beliefs about the problem - :param optimize_restarts: number of optimization restarts to use when training the models - """ - super(ParallelBatchAcquisition, self).__init__() - models = np.atleast_1d(models) - assert all(isinstance(model, (Model, ModelWrapper)) for model in models) - self._models = ParamList([DataScaler(m) for m in models]) + In Bayesian Optimization this function is typically optimized over the optimization domain + to determine the next point for evaluation. An objects satisfying this interface hold a list of GPflow models. + Implementing build_acquisition yields the score of the acquisition function (usually obtained by transforming the + predictive distribution) using TensorFlow. Optionally, a method _setup can be implemented which computes some + quantities which are used to compute the acquisition, do not vary with the candidate points. This method is not + automatically called, classes implementing this interface may implement their own mechanism on when _setup is + called. - assert (optimize_restarts >= 0) - self.optimize_restarts = optimize_restarts - self._needs_setup = True - self.batch_size = batch_size + 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. + """ + __metaclass__ = ABCMeta + + @abstractmethod def _optimize_models(self): """ Optimizes the hyperparameters of all models that the acquisition function is based on. @@ -85,6 +89,184 @@ def _optimize_models(self): As a special case, if optimize_restarts is set to zero, the hyperparameters of the models are not optimized. This is useful when the hyperparameters are sampled using MCMC. """ + pass + + @abstractmethod + 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 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 + + @abstractmethod + def _setup_constraints(self): + """ + Run only if some outputs handled by this acquisition are constraints. Used in aggregation. + """ + pass + + @abstractmethod + def _setup_objectives(self): + """ + Run only if all outputs handled by this acquisition are objectives. Used in aggregation. + """ + pass + + @abstractmethod + def get_suggestion(self, optimizer): + pass + + @abstractmethod + def build_acquisition(self, *args): + pass + + @abstractmethod + 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 acquisition function. Note that the models are modified directly and + references to them outside of the object will also point to scaled instances. + + :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. + """ + pass + + @abstractproperty + def models(self): + """ + The GPflow models representing our beliefs of the optimization problem. + + :return: list of GPflow models + """ + pass + + @abstractproperty + def data(self): + """ + The training data of the models. + + Corresponds to the input data X which is the same for every model, + and column-wise concatenation of the Y data over all models + + :return: tuple X, Y of tensors (if in tf_mode) or numpy arrays. + """ + pass + + @abstractmethod + def set_data(self, X, Y): + """ + 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 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 + + :param X: input data N x D + :param Y: output data N x R (R >= Q) + :return: Q (sum of output dimensions of contained models) + """ + pass + + @abstractmethod + def constraint_indices(self): + """ + Method returning the indices of the model outputs which correspond to the (expensive) constraint functions. + By default there are no constraint functions + """ + pass + + @abstractmethod + def objective_indices(self): + """ + Method returning the indices of the model outputs which are objective functions. + + By default all outputs are objectives. + + :return: indices to the objectives, size R + """ + pass + + @abstractmethod + def feasible_data_index(self): + """ + Returns a boolean array indicating which data points are considered feasible (according to the acquisition + function(s) ) and which not. + + By default all data is considered feasible. + + :return: logical indices to the feasible data points, size N + """ + pass + + @abstractmethod + def evaluate(self, candidates): + pass + + @abstractmethod + def evaluate_with_gradients(self, candidates): + pass + + @abstractmethod + def __add__(self, other): + """ + Operator for adding acquisition functions. Example: + + >>> a1 = gpflowopt.acquisition.ExpectedImprovement(m1) + >>> a2 = gpflowopt.acquisition.ProbabilityOfFeasibility(m2) + >>> type(a1 + a2) + + """ + pass + + @abstractmethod + def __mul__(self, other): + """ + Operator for multiplying acquisition functions. Example: + + >>> a1 = gpflowopt.acquisition.ExpectedImprovement(m1) + >>> a2 = gpflowopt.acquisition.ProbabilityOfFeasibility(m2) + >>> type(a1 * a2) + + """ + pass + + +class ParallelBatchAcquisition(IAcquisition): + """ + ParallelBatchAcquisition objects 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. + """ + + __metaclass__ = ABCMeta + + def __init__(self, models=[], optimize_restarts=5, batch_size=1): + """ + :param models: list of GPflow models representing our beliefs about the problem + :param optimize_restarts: number of optimization restarts to use when training the models + """ + super(ParallelBatchAcquisition, self).__init__() + models = np.atleast_1d(models) + assert all(isinstance(model, (Model, ModelWrapper)) for model in models) + self._models = ParamList([DataScaler(m) for m in models]) + + assert (optimize_restarts >= 0) + self.optimize_restarts = optimize_restarts + self._needs_setup = True + self.batch_size = batch_size + + def _optimize_models(self): if self.optimize_restarts == 0: return @@ -113,20 +295,7 @@ def get_suggestion(self, optimizer): result = opt.optimize(self._inverse_acquisition) return np.vstack(np.split(result.x, self.batch_size, axis=1)) - def build_acquisition(self, *args): - raise NotImplementedError - 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. Note that the models are modified directly and - references to them outside of the object will also point to scaled instances. - - :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. - """ n_inputs = self.data[0].shape[1] assert (domain.size == n_inputs) for m in self.models: @@ -135,19 +304,6 @@ def enable_scaling(self, domain): self.highest_parent._needs_setup = True def set_data(self, X, Y): - """ - 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 - - :param X: input data N x D - :param Y: output data N x R (R >= Q) - :return: Q (sum of output dimensions of contained models) - """ num_outputs_sum = 0 for model in self.models: num_outputs = model.Y.shape[1] @@ -162,77 +318,32 @@ def set_data(self, X, Y): @property def models(self): - """ - The GPflow models representing our beliefs of the optimization problem. - - :return: list of GPflow models - """ return self._models.sorted_params @property def data(self): - """ - The training data of the models. - - Corresponds to the input data X which is the same for every model, - and column-wise concatenation of the Y data over all models - - :return: tuple X, Y of tensors (if in tf_mode) or numpy arrays. - """ if self._tf_mode: return self.models[0].X, tf.concat(list(map(lambda model: model.Y, self.models)), 1) else: return self.models[0].X.value, np.hstack(map(lambda model: model.Y.value, self.models)) def constraint_indices(self): - """ - Method returning the indices of the model outputs which correspond to the (expensive) constraint functions. - By default there are no constraint functions - """ return np.empty((0,), dtype=int) def objective_indices(self): - """ - Method returning the indices of the model outputs which are objective functions. - - By default all outputs are objectives. - - :return: indices to the objectives, size R - """ return np.setdiff1d(np.arange(self.data[1].shape[1]), self.constraint_indices()) def feasible_data_index(self): - """ - Returns a boolean array indicating which data points are considered feasible (according to the acquisition - function(s) ) and which not. - - By default all data is considered feasible. - - :return: logical indices to the feasible data points, size N - """ return np.ones(self.data[0].shape[0], dtype=bool) 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 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() @@ -247,34 +358,62 @@ def evaluate_with_gradients(self, Xcand): def evaluate(self, Xcand): return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) + @abstractmethod + def build_acquisition(self, *args): + pass -class Acquisition(ParallelBatchAcquisition): - """ - An acquisition function maps the belief represented by a Bayesian model into a - score indicating how promising a point is for evaluation. + def __add__(self, other): + if not isinstance(other, IAcquisition): + return NotImplemented - 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. 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. + if isinstance(other, Acquisition): + other = ParToSeqAcquisitionWrapper(other, self.batch_size) - 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. + return AcquisitionSum([self, other]) + + def __mul__(self, other): + if not isinstance(other, IAcquisition): + return NotImplemented + + if isinstance(other, Acquisition): + other = ParToSeqAcquisitionWrapper(other, self.batch_size) + + return AcquisitionProduct([self, other]) + + def __radd__(self, other): + if not isinstance(other, IAcquisition): + return NotImplemented + + if isinstance(other, Acquisition): + other = ParToSeqAcquisitionWrapper(other, self.batch_size) + + return AcquisitionSum([other, self]) - 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 __rmul__(self, other): + if not isinstance(other, IAcquisition): + return NotImplemented + + if isinstance(other, Acquisition): + other = ParToSeqAcquisitionWrapper(other, self.batch_size) + + return AcquisitionProduct([other, self]) + + + +class Acquisition(ParallelBatchAcquisition): """ + Class to be implemented for standard single-point acquisition functions like standard EI. + """ + + __metaclass__ = ABCMeta def get_suggestion(self, optimizer): result = optimizer.optimize(self._inverse_acquisition) return result.x + @abstractmethod def build_acquisition(self, candidates): - raise NotImplementedError + pass @setup_required @AutoFlow((float_type, [None, None])) @@ -299,29 +438,15 @@ def evaluate(self, Xcand): return self.build_acquisition(Xcand) def __add__(self, other): - """ - Operator for adding acquisition functions. Example: + if not isinstance(other, Acquisition): + return NotImplemented - >>> a1 = gpflowopt.acquisition.ExpectedImprovement(m1) - >>> a2 = gpflowopt.acquisition.ProbabilityOfFeasibility(m2) - >>> type(a1 + a2) - - """ - if isinstance(other, AcquisitionSum): - return AcquisitionSum([self] + other.operands.sorted_params) return AcquisitionSum([self, other]) def __mul__(self, other): - """ - Operator for multiplying acquisition functions. Example: + if not isinstance(other, Acquisition): + return NotImplemented - >>> a1 = gpflowopt.acquisition.ExpectedImprovement(m1) - >>> a2 = gpflowopt.acquisition.ProbabilityOfFeasibility(m2) - >>> type(a1 * a2) - - """ - if isinstance(other, AcquisitionProduct): - return AcquisitionProduct([self] + other.operands.sorted_params) return AcquisitionProduct([self, other]) def __setattr__(self, key, value): @@ -330,46 +455,32 @@ def __setattr__(self, key, value): self.highest_parent._needs_setup = True -class AcquisitionAggregation(Acquisition): - """ - Aggregates multiple acquisition functions, using a TensorFlow reduce operation. - """ +class AcquisitionWrapper(IAcquisition): - def __init__(self, operands, oper): + __metaclass__ = ABCMeta + + def __init__(self, wrap, *args, **kwargs): """ - :param operands: list of acquisition objects - :param oper: a tf.reduce operation (e.g., tf.reduce_sum) for aggregating the returned scores of each operand. + :param wrap: list of acquisition objects """ - super(AcquisitionAggregation, self).__init__() - assert (all([isinstance(x, Acquisition) for x in operands])) - self.operands = ParamList(operands) - self._oper = oper + self.wrapped = [] + super(AcquisitionWrapper, self).__init__() + assert all([isinstance(x, IAcquisition) for x in wrap]) + assert all(wrap[0].batch_size == oper.batch_size for oper in wrap) + self.wrapped = ParamList(wrap) + def _optimize_models(self): - for oper in self.operands: + for oper in self.wrapped: oper._optimize_models() - @Acquisition.models.getter - def models(self): - return [model for acq in self.operands for model in acq.models] - - def enable_scaling(self, domain): - for oper in self.operands: - oper.enable_scaling(domain) - - def set_data(self, X, Y): - offset = 0 - for operand in self.operands: - offset += operand.set_data(X, Y[:, offset:]) - return offset - def _setup_constraints(self): - for oper in self.operands: + for oper in self.wrapped: 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: + for oper in self.wrapped: oper._setup_objectives() def _setup(self): @@ -378,24 +489,73 @@ def _setup(self): # Then objectives as these might depend on the constraint acquisition self._setup_objectives() + @abstractmethod + def build_acquisition(self, *args): + pass + + @Acquisition.models.getter + def models(self): + return [model for acq in self.wrapped for model in acq.models] + + def enable_scaling(self, domain): + for oper in self.wrapped: + oper.enable_scaling(domain) + + def set_data(self, X, Y): + offset = 0 + for operand in self.wrapped: + offset += operand.set_data(X, Y[:, offset:]) + return offset + def constraint_indices(self): offset = [0] idx = [] - for operand in self.operands: + for operand in self.wrapped: idx.append(operand.constraint_indices()) offset.append(operand.data[1].shape[1]) return np.hstack([i + o for i, o in zip(idx, offset[:-1])]) def feasible_data_index(self): - return np.all(np.vstack(map(lambda o: o.feasible_data_index(), self.operands)), axis=0) + return np.all(np.vstack(map(lambda o: o.feasible_data_index(), self.wrapped)), axis=0) + + def __getitem__(self, item): + return self.wrapped[item] + + +class AcquisitionAggregation(AcquisitionWrapper, ParallelBatchAcquisition): + """ + Aggregates multiple acquisition functions, using a TensorFlow reduce operation. + """ + + def __init__(self, operands, oper): + """ + :param operands: list of acquisition objects + :param oper: a tf.reduce operation (e.g., tf.reduce_sum) for aggregating the returned scores of each operand. + """ + super(AcquisitionAggregation, self).__init__(operands) + self._oper = oper + + @property + def batch_size(self): + return self[0].batch_size + + @batch_size.setter + def batch_size(self, value): + for oper in self.operands: + oper.batch_size = value + + @property + def operands(self): + return self.wrapped + + @operands.setter + def operands(self, value): + self.wrapped = value def build_acquisition(self, *args): return self._oper(tf.concat(list(map(lambda operand: operand.build_acquisition(*args), self.operands)), 1), axis=1, keep_dims=True, name=self.__class__.__name__) - def __getitem__(self, item): - return self.operands[item] - class AcquisitionSum(AcquisitionAggregation): """ @@ -410,6 +570,9 @@ def __add__(self, other): else: return AcquisitionSum(self.operands.sorted_params + [other]) + def __radd__(self, other): + return AcquisitionSum([other] + self.operands.sorted_params) + class AcquisitionProduct(AcquisitionAggregation): """ @@ -424,6 +587,9 @@ def __mul__(self, other): else: return AcquisitionProduct(self.operands.sorted_params + [other]) + def __rmul__(self, other): + return AcquisitionProduct([other] + self.operands.sorted_params) + class MCMCAcquistion(AcquisitionSum): """ @@ -486,4 +652,18 @@ def _kill_autoflow(self): cause inconsistencies and errors. """ super(MCMCAcquistion, self)._kill_autoflow() - self._needs_new_copies = True \ No newline at end of file + self._needs_new_copies = True + + +class ParToSeqAcquisitionWrapper(AcquisitionWrapper, ParallelBatchAcquisition): + + def __init__(self, acquisition, batch_size): + assert isinstance(acquisition, Acquisition) + super(ParToSeqAcquisitionWrapper, self).__init__([acquisition], batch_size=batch_size) + + def build_acquisition(self, *args): + assert len(args) == self.batch_size + all_scores = self.wrapped.build_acquisition(tf.concat(args, 0)) + return tf.reduce_prod(tf.concat(tf.split(all_scores, axis=0, num_or_size_splits=self.batch_size), 1), axis=1) + + diff --git a/gpflowopt/bo.py b/gpflowopt/bo.py index 6fc4e0e..84f309c 100644 --- a/gpflowopt/bo.py +++ b/gpflowopt/bo.py @@ -13,14 +13,13 @@ # limitations under the License. from contextlib import contextmanager -import copy import numpy as np from scipy.optimize import OptimizeResult import tensorflow as tf from gpflow.gpr import GPR -from .acquisition import Acquisition, MCMCAcquistion +from .acquisition import IAcquisition, MCMCAcquistion from .design import Design, EmptyDesign from .objective import ObjectiveWrapper from .optim import Optimizer, SciPyOptimizer @@ -90,7 +89,7 @@ def __init__(self, domain, acquisition, optimizer=None, initial=None, scaling=Tr :class:`~.Acquisition` this allows several scenarios: do the optimization manually from the callback (optimize_restarts equals 0), or choose the starting point + some random restarts (optimize_restarts > 0). """ - assert isinstance(acquisition, Acquisition) + assert isinstance(acquisition, IAcquisition) assert hyper_draws is None or hyper_draws > 0 assert optimizer is None or isinstance(optimizer, Optimizer) assert initial is None or isinstance(initial, Design) diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index f406acd..a2762be 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -333,37 +333,37 @@ def test_multi_aggr(self): acq1, acq2, acq3, acq4 = acq joint = acq1 + acq2 + acq3 self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) - self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3]) + self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3]) joint = acq1 * acq2 * acq3 self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) - self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3]) + self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3]) first = acq2 + acq3 self.assertIsInstance(first, gpflowopt.acquisition.AcquisitionSum) - self.assertListEqual(first.operands.sorted_params, [acq2, acq3]) + self.assertListEqual(first.wrapped.sorted_params, [acq2, acq3]) joint = acq1 + first self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) - self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3]) + self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3]) first = acq2 * acq3 self.assertIsInstance(first, gpflowopt.acquisition.AcquisitionProduct) - self.assertListEqual(first.operands.sorted_params, [acq2, acq3]) + self.assertListEqual(first.wrapped.sorted_params, [acq2, acq3]) joint = acq1 * first self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) - self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3]) + self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3]) first = acq1 + acq2 second = acq3 + acq4 joint = first + second self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) - self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3, acq4]) + self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3, acq4]) first = acq1 * acq2 second = acq3 * acq4 joint = first * second self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) - self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3, acq4]) + self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3, acq4]) class TestRecompile(GPflowOptTestCase): @@ -377,7 +377,7 @@ def test_vgp(self): Y = np.sin(X[:,[0]]) m = gpflow.vgp.VGP(X, Y, gpflow.kernels.RBF(2), gpflow.likelihoods.Gaussian()) acq = gpflowopt.acquisition.ExpectedImprovement(m) - m._compile() + m.compile() self.assertFalse(m._needs_recompile) acq.evaluate(gpflowopt.design.RandomDesign(10, domain).generate()) self.assertTrue(hasattr(acq, '_evaluate_AF_storage')) From ec477e39896f3c4100fb11a24cb2d042d9dc895b Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Mon, 18 Sep 2017 00:20:38 +0200 Subject: [PATCH 09/15] Reworked some interfaces, started testing generating operators and restructured some tests --- gpflowopt/acquisition/__init__.py | 2 +- gpflowopt/acquisition/acquisition.py | 92 +++++++++----- testing/test_acquisition.py | 173 ++++++++++++++++++++------- testing/test_implementations.py | 27 +++++ 4 files changed, 221 insertions(+), 73 deletions(-) diff --git a/gpflowopt/acquisition/__init__.py b/gpflowopt/acquisition/__init__.py index a8f4dba..b9d2bb8 100644 --- a/gpflowopt/acquisition/__init__.py +++ b/gpflowopt/acquisition/__init__.py @@ -14,7 +14,7 @@ # Framework components and interfaces from .acquisition import (Acquisition, ParallelBatchAcquisition, AcquisitionAggregation, AcquisitionProduct, \ - AcquisitionSum, MCMCAcquistion, IAcquisition, setup_required) + AcquisitionSum, MCMCAcquistion, IAcquisition, ParToSeqAcquisitionWrapper, setup_required) # Single objective from .ei import ExpectedImprovement diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index 72b8698..7a0d8dd 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -55,7 +55,7 @@ def runnable(instance, *args, **kwargs): return runnable -class IAcquisition(Parameterized): # pragma: no cover +class IAcquisition(Parameterized): """ Interface for Acquisition functions mapping the belief represented by a Bayesian model into a score indicating how promising a point is for evaluation. @@ -159,6 +159,10 @@ def data(self): """ pass + @abstractproperty + def batch_size(self): + pass + @abstractmethod def set_data(self, X, Y): """ @@ -209,10 +213,21 @@ def feasible_data_index(self): @abstractmethod def evaluate(self, candidates): + """ + AutoFlow method to compute the acquisition scores for candidates, without returning the gradients. + + :return: acquisition scores, size N x 1 + """ pass @abstractmethod def evaluate_with_gradients(self, candidates): + """ + AutoFlow method to compute the acquisition scores for candidates, also returns the gradients. + + :return: acquisition scores, size N x 1 + the gradients of the acquisition scores, size N x D + """ pass @abstractmethod @@ -264,7 +279,7 @@ def __init__(self, models=[], optimize_restarts=5, batch_size=1): assert (optimize_restarts >= 0) self.optimize_restarts = optimize_restarts self._needs_setup = True - self.batch_size = batch_size + self._batch_size = batch_size def _optimize_models(self): if self.optimize_restarts == 0: @@ -327,6 +342,14 @@ def data(self): else: return self.models[0].X.value, np.hstack(map(lambda model: model.Y.value, self.models)) + @property + def batch_size(self): + return self._batch_size + + @batch_size.setter + def batch_size(self, value): + self._batch_size = value + def constraint_indices(self): return np.empty((0,), dtype=int) @@ -366,16 +389,27 @@ def __add__(self, other): if not isinstance(other, IAcquisition): return NotImplemented - if isinstance(other, Acquisition): - other = ParToSeqAcquisitionWrapper(other, self.batch_size) + if isinstance(other, AcquisitionSum): + return NotImplemented + + if self.batch_size < other.batch_size: + return ParToSeqAcquisitionWrapper(self, other.batch_size) + other + if self.batch_size > other.batch_size: + other = ParToSeqAcquisitionWrapper(other, self.batch_size) return AcquisitionSum([self, other]) def __mul__(self, other): if not isinstance(other, IAcquisition): return NotImplemented - if isinstance(other, Acquisition): + if isinstance(other, AcquisitionProduct): + return NotImplemented + + if self.batch_size < other.batch_size: + return ParToSeqAcquisitionWrapper(self, other.batch_size) * other + + if self.batch_size > other.batch_size: other = ParToSeqAcquisitionWrapper(other, self.batch_size) return AcquisitionProduct([self, other]) @@ -384,7 +418,7 @@ def __radd__(self, other): if not isinstance(other, IAcquisition): return NotImplemented - if isinstance(other, Acquisition): + if self.batch_size > other.batch_size: other = ParToSeqAcquisitionWrapper(other, self.batch_size) return AcquisitionSum([other, self]) @@ -393,13 +427,12 @@ def __rmul__(self, other): if not isinstance(other, IAcquisition): return NotImplemented - if isinstance(other, Acquisition): + if self.batch_size > other.batch_size: other = ParToSeqAcquisitionWrapper(other, self.batch_size) return AcquisitionProduct([other, self]) - class Acquisition(ParallelBatchAcquisition): """ Class to be implemented for standard single-point acquisition functions like standard EI. @@ -418,23 +451,12 @@ def build_acquisition(self, candidates): @setup_required @AutoFlow((float_type, [None, None])) def evaluate_with_gradients(self, Xcand): - """ - AutoFlow method to compute the acquisition scores for candidates, also returns the gradients. - - :return: acquisition scores, size N x 1 - the gradients of the acquisition scores, size N x D - """ 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): - """ - AutoFlow method to compute the acquisition scores for candidates, without returning the gradients. - - :return: acquisition scores, size N x 1 - """ return self.build_acquisition(Xcand) def __add__(self, other): @@ -464,12 +486,11 @@ def __init__(self, wrap, *args, **kwargs): :param wrap: list of acquisition objects """ self.wrapped = [] - super(AcquisitionWrapper, self).__init__() + super(AcquisitionWrapper, self).__init__(*args, **kwargs) assert all([isinstance(x, IAcquisition) for x in wrap]) assert all(wrap[0].batch_size == oper.batch_size for oper in wrap) self.wrapped = ParamList(wrap) - def _optimize_models(self): for oper in self.wrapped: oper._optimize_models() @@ -535,7 +556,7 @@ def __init__(self, operands, oper): super(AcquisitionAggregation, self).__init__(operands) self._oper = oper - @property + @ParallelBatchAcquisition.batch_size.getter def batch_size(self): return self[0].batch_size @@ -565,10 +586,13 @@ def __init__(self, operands): super(AcquisitionSum, self).__init__(operands, tf.reduce_sum) def __add__(self, other): - if isinstance(other, AcquisitionSum): - return AcquisitionSum(self.operands.sorted_params + other.operands.sorted_params) + if self.batch_size == other.batch_size: + if isinstance(other, AcquisitionSum): + return AcquisitionSum(self.operands.sorted_params + other.operands.sorted_params) + else: + return AcquisitionSum(self.operands.sorted_params + [other]) else: - return AcquisitionSum(self.operands.sorted_params + [other]) + super(AcquisitionSum, self).__add__(other) def __radd__(self, other): return AcquisitionSum([other] + self.operands.sorted_params) @@ -582,10 +606,13 @@ def __init__(self, operands): super(AcquisitionProduct, self).__init__(operands, tf.reduce_prod) def __mul__(self, other): - if isinstance(other, AcquisitionProduct): - return AcquisitionProduct(self.operands.sorted_params + other.operands.sorted_params) + if self.batch_size == other.batch_size: + if isinstance(other, AcquisitionProduct): + return AcquisitionProduct(self.operands.sorted_params + other.operands.sorted_params) + else: + return AcquisitionProduct(self.operands.sorted_params + [other]) else: - return AcquisitionProduct(self.operands.sorted_params + [other]) + super(AcquisitionSum, self).__add__(other) def __rmul__(self, other): return AcquisitionProduct([other] + self.operands.sorted_params) @@ -658,12 +685,13 @@ def _kill_autoflow(self): class ParToSeqAcquisitionWrapper(AcquisitionWrapper, ParallelBatchAcquisition): def __init__(self, acquisition, batch_size): - assert isinstance(acquisition, Acquisition) + assert isinstance(acquisition, IAcquisition) super(ParToSeqAcquisitionWrapper, self).__init__([acquisition], batch_size=batch_size) def build_acquisition(self, *args): - assert len(args) == self.batch_size - all_scores = self.wrapped.build_acquisition(tf.concat(args, 0)) - return tf.reduce_prod(tf.concat(tf.split(all_scores, axis=0, num_or_size_splits=self.batch_size), 1), axis=1) + splits = self.batch_size // self.wrapped[0].batch_size + assert len(args) == splits + all_scores = self.wrapped[0].build_acquisition(tf.concat(args, 0)) + return tf.reduce_prod(tf.concat(tf.split(all_scores, axis=0, num_or_size_splits=splits), 1), axis=1) diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index a2762be..11b5877 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -7,16 +7,25 @@ domain = np.sum([gpflowopt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)]) - +####### TESTING CLASSES ####### class SimpleAcquisition(gpflowopt.acquisition.Acquisition): def build_acquisition(self, Xcand): return -self.models[0].build_predict(Xcand)[0] +class SimpleAcquisitionConstrained(gpflowopt.acquisition.Acquisition): + + def constraint_indices(self): + return np.arange(self.data[1].shape[1]) + + def build_acquisition(self, Xcand): + return -self.models[0].build_predict(Xcand)[0] + + class SimpleParallelBatch(gpflowopt.acquisition.ParallelBatchAcquisition): - def __init__(self, model): - super(SimpleParallelBatch, self).__init__(model, batch_size=2) + def __init__(self, model, batch_size=2): + super(SimpleParallelBatch, self).__init__(model, batch_size=batch_size) self.n_args = 0 self.counter = 0 @@ -27,6 +36,7 @@ def _setup(self): def build_acquisition(self, *args): self.n_args = len(args) return -tf.add(self.models[0].build_predict(args[0])[0], self.models[0].build_predict(args[1]+1)[0]) +############################## class TestParallelBatchAcquisition(GPflowOptTestCase): @@ -174,10 +184,11 @@ def test_get_suggestion(self): ) -class TestAcquisitionAggregation(GPflowOptTestCase): +class TestAcquisitionAggregationImplementations(GPflowOptTestCase): @parameterized.expand(list(zip(aggregations))) def test_object_integrity(self, acquisition): + self.assertTrue(isinstance(acquisition, gpflowopt.acquisition.IAcquisition)) acquisition._kill_autoflow() with self.test_session(): for oper in acquisition.operands: @@ -236,15 +247,6 @@ 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): acquisition._kill_autoflow() @@ -287,32 +289,124 @@ def test_mcmc_acq(self): np.testing.assert_almost_equal(ei_mle, ei_mcmc, decimal=5) +class TestParallelToSequentialAdapter(GPflowOptTestCase): + + def test_wrapper(self): + acq = SimpleAcquisition(create_parabola_model(domain)) + adapter = gpflowopt.acquisition.ParToSeqAcquisitionWrapper(acq, 3) + self.assertEqual(acq.batch_size, 1) + self.assertEqual(adapter.batch_size, 3) + + a = np.array([[0.25, -0.25]]) + b = np.array([[0.67, 0]]) + c = np.array([[-0.1, 0.1]]) + + r1 = acq.evaluate(np.vstack((a, b, c))) + r2 = adapter.evaluate(np.hstack((a, b, c))) + self.assertTrue(np.allclose(np.prod(r1), r2)) + + class TestJointAcquisition(GPflowOptTestCase): - def test_constrained_ei(self): - with self.test_session(): - design = gpflowopt.design.LatinHyperCube(16, domain) - X = design.generate() - Yo = parabola2d(X) - Yc = -parabola2d(X) + 0.5 - m1 = gpflow.gpr.GPR(X, Yo, gpflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) - m2 = gpflow.gpr.GPR(X, Yc, gpflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) - ei = gpflowopt.acquisition.ExpectedImprovement(m1) - 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)) + def test_simple_generating_operators(self): + joint = SimpleAcquisition(create_parabola_model(domain)) + SimpleAcquisition(create_parabola_model(domain)) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + + joint = SimpleAcquisition(create_parabola_model(domain)) * SimpleAcquisition(create_parabola_model(domain)) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + + def test_batch_generating_operators(self): + joint = SimpleParallelBatch(create_parabola_model(domain), 2) + \ + SimpleParallelBatch(create_parabola_model(domain), 2) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + self.assertTrue(isinstance(joint.operands[0], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + + joint = SimpleParallelBatch(create_parabola_model(domain), 2) * \ + SimpleParallelBatch(create_parabola_model(domain), 2) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + self.assertTrue(isinstance(joint.operands[0], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + + joint = SimpleParallelBatch(create_parabola_model(domain), 6) + \ + SimpleParallelBatch(create_parabola_model(domain), 2) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + self.assertTrue(isinstance(joint.operands[0], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[1], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 6) + self.assertEqual(joint.operands[0].batch_size, 6) + self.assertEqual(joint.operands[1].batch_size, 6) + + joint = SimpleParallelBatch(create_parabola_model(domain), 2) + \ + SimpleParallelBatch(create_parabola_model(domain), 6) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 6) + self.assertEqual(joint.operands[0].batch_size, 6) + self.assertEqual(joint.operands[1].batch_size, 6) + + joint = SimpleParallelBatch(create_parabola_model(domain), 6) * \ + SimpleParallelBatch(create_parabola_model(domain), 2) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + self.assertTrue(isinstance(joint.operands[0], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[1], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 6) + self.assertEqual(joint.operands[0].batch_size, 6) + self.assertEqual(joint.operands[1].batch_size, 6) + + joint = SimpleParallelBatch(create_parabola_model(domain), 2) * \ + SimpleParallelBatch(create_parabola_model(domain), 6) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 6) + self.assertEqual(joint.operands[0].batch_size, 6) + self.assertEqual(joint.operands[1].batch_size, 6) + + def test_mixed_generating_operators(self): + joint = SimpleAcquisition(create_parabola_model(domain)) + SimpleParallelBatch(create_parabola_model(domain), 3) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionSum)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 3) + self.assertEqual(joint.operands[0].batch_size, 3) + self.assertEqual(joint.operands[1].batch_size, 3) + + joint = SimpleAcquisition(create_parabola_model(domain)) * SimpleParallelBatch(create_parabola_model(domain), 3) + self.assertTrue(isinstance(joint, gpflowopt.acquisition.AcquisitionProduct)) + self.assertTrue(isinstance(joint.operands[1], SimpleParallelBatch)) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.batch_size, 3) + self.assertEqual(joint.operands[0].batch_size, 3) + self.assertEqual(joint.operands[1].batch_size, 3) + + def test_incorrect_generating_operators(self): + with self.assertRaises(TypeError): + SimpleAcquisition(create_parabola_model(domain)) + 1 + + with self.assertRaises(TypeError): + 1 + SimpleAcquisition(create_parabola_model(domain)) + + with self.assertRaises(TypeError): + SimpleAcquisition(create_parabola_model(domain)) * 1 + + with self.assertRaises(TypeError): + 1 * SimpleAcquisition(create_parabola_model(domain)) + + with self.assertRaises(TypeError): + SimpleParallelBatch(create_parabola_model(domain), 2) + 1 + + with self.assertRaises(TypeError): + 1 + SimpleParallelBatch(create_parabola_model(domain), 2) + + with self.assertRaises(TypeError): + SimpleParallelBatch(create_parabola_model(domain), 2) * 1 - # 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") + with self.assertRaises(TypeError): + 1 * SimpleParallelBatch(create_parabola_model(domain), 2) - def test_hierarchy(self): + def test_simple_hierarchy(self): with self.test_session(): design = gpflowopt.design.LatinHyperCube(16, domain) X = design.generate() @@ -320,16 +414,15 @@ def test_hierarchy(self): 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) - + gpflowopt.acquisition.ExpectedImprovement(m2)) + joint = SimpleAcquisition(m1) * (SimpleAcquisitionConstrained(m3) + SimpleAcquisition(m2)) np.testing.assert_allclose(joint.objective_indices(), np.array([0, 2], dtype=int)) np.testing.assert_allclose(joint.constraint_indices(), np.array([1], dtype=int)) - def test_multi_aggr(self): + @parameterized.expand([([SimpleAcquisition(create_parabola_model(domain)) for i in range(4)],), + ([SimpleParallelBatch(create_parabola_model(domain)) for i in range(4)],)]) + def test_multi_aggr_simple(self, acq): with self.test_session(): - 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 index 3f3c3df..1850153 100644 --- a/testing/test_implementations.py +++ b/testing/test_implementations.py @@ -1,4 +1,5 @@ import gpflowopt +import gpflow import numpy as np from parameterized import parameterized from .utility import create_parabola_model, create_plane_model, create_vlmop2_model, parabola2d, load_data, GPflowOptTestCase @@ -186,3 +187,29 @@ def test_hvpoi_validity(self): with self.test_session(): scores = self.acquisition.evaluate(self.data['candidates']) np.testing.assert_almost_equal(scores, self.data['scores'], decimal=2) + + +class TestConstrainedExpectedImprovement(GPflowOptTestCase): + + def test_constrained_ei(self): + with self.test_session(): + design = gpflowopt.design.LatinHyperCube(16, domain) + X = design.generate() + Yo = parabola2d(X) + Yc = -parabola2d(X) + 0.5 + m1 = gpflow.gpr.GPR(X, Yo, gpflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m2 = gpflow.gpr.GPR(X, Yc, gpflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + ei = gpflowopt.acquisition.ExpectedImprovement(m1) + 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") \ No newline at end of file From 7c7649e57cb1b6873e2791305d5a4c6325d31b65 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Mon, 18 Sep 2017 00:34:47 +0200 Subject: [PATCH 10/15] Cleaned up file and added some documentation --- gpflowopt/acquisition/acquisition.py | 153 +++++++++++++++------------ 1 file changed, 87 insertions(+), 66 deletions(-) diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index 7a0d8dd..5ebc75f 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -255,8 +255,93 @@ def __mul__(self, other): pass +class AcquisitionWrapper(IAcquisition): + """ + Partial implementation of the Acquisition interface, useful for implementing acquisition functions which wrap + another acquisition function (and forward method calls). + + Typically used in combination with multiple inheritance: + + >>> class WrappedAcquisition(AcquisitionWrapper, Acquisition) + >>> def __init__(self, acquisition): + >>> super(WrappedAcquisition, self).__init__(acquisition, optimize_restarts=2) + >>> + >>> def build_acquisition(self, candidates): + >>> return tf.constant(1.0) + """ + + __metaclass__ = ABCMeta + + def __init__(self, wrap, *args, **kwargs): + """ + :param wrap: list of acquisition objects + """ + self.wrapped = [] + super(AcquisitionWrapper, self).__init__(*args, **kwargs) + assert all([isinstance(x, IAcquisition) for x in wrap]) + assert all(wrap[0].batch_size == oper.batch_size for oper in wrap) + self.wrapped = ParamList(np.atleast_1d(wrap).tolist()) + + def _optimize_models(self): + for oper in self.wrapped: + oper._optimize_models() + + def _setup_constraints(self): + for oper in self.wrapped: + if oper.constraint_indices().size > 0: # Small optimization, skip subtrees with objectives only + oper._setup_constraints() + + def _setup_objectives(self): + for oper in self.wrapped: + 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() + + @abstractmethod + def build_acquisition(self, *args): + pass + + @Acquisition.models.getter + def models(self): + return [model for acq in self.wrapped for model in acq.models] + + def enable_scaling(self, domain): + for oper in self.wrapped: + oper.enable_scaling(domain) + + def set_data(self, X, Y): + offset = 0 + for operand in self.wrapped: + offset += operand.set_data(X, Y[:, offset:]) + return offset + + def constraint_indices(self): + offset = [0] + idx = [] + for operand in self.wrapped: + idx.append(operand.constraint_indices()) + offset.append(operand.data[1].shape[1]) + return np.hstack([i + o for i, o in zip(idx, offset[:-1])]) + + def feasible_data_index(self): + return np.all(np.vstack(map(lambda o: o.feasible_data_index(), self.wrapped)), axis=0) + + def __getitem__(self, item): + return self.wrapped[item] + + class ParallelBatchAcquisition(IAcquisition): """ + Abstract Acquisition class for batch acquisition functions, optimizing the batch in parallel (as opposed to + sequential or greedy approaches). + + Subclasses deriving this class should accept multiple candidates tensors for build_acquisition. Each tensor + represents points x1, x2, ... xn in the batch. + ParallelBatchAcquisition objects 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 @@ -436,6 +521,8 @@ def __rmul__(self, other): class Acquisition(ParallelBatchAcquisition): """ Class to be implemented for standard single-point acquisition functions like standard EI. + + This abstract class inherits the lazy setup mechanism from its ancestor. """ __metaclass__ = ABCMeta @@ -477,72 +564,6 @@ def __setattr__(self, key, value): self.highest_parent._needs_setup = True -class AcquisitionWrapper(IAcquisition): - - __metaclass__ = ABCMeta - - def __init__(self, wrap, *args, **kwargs): - """ - :param wrap: list of acquisition objects - """ - self.wrapped = [] - super(AcquisitionWrapper, self).__init__(*args, **kwargs) - assert all([isinstance(x, IAcquisition) for x in wrap]) - assert all(wrap[0].batch_size == oper.batch_size for oper in wrap) - self.wrapped = ParamList(wrap) - - def _optimize_models(self): - for oper in self.wrapped: - oper._optimize_models() - - def _setup_constraints(self): - for oper in self.wrapped: - if oper.constraint_indices().size > 0: # Small optimization, skip subtrees with objectives only - oper._setup_constraints() - - def _setup_objectives(self): - for oper in self.wrapped: - 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() - - @abstractmethod - def build_acquisition(self, *args): - pass - - @Acquisition.models.getter - def models(self): - return [model for acq in self.wrapped for model in acq.models] - - def enable_scaling(self, domain): - for oper in self.wrapped: - oper.enable_scaling(domain) - - def set_data(self, X, Y): - offset = 0 - for operand in self.wrapped: - offset += operand.set_data(X, Y[:, offset:]) - return offset - - def constraint_indices(self): - offset = [0] - idx = [] - for operand in self.wrapped: - idx.append(operand.constraint_indices()) - offset.append(operand.data[1].shape[1]) - return np.hstack([i + o for i, o in zip(idx, offset[:-1])]) - - def feasible_data_index(self): - return np.all(np.vstack(map(lambda o: o.feasible_data_index(), self.wrapped)), axis=0) - - def __getitem__(self, item): - return self.wrapped[item] - - class AcquisitionAggregation(AcquisitionWrapper, ParallelBatchAcquisition): """ Aggregates multiple acquisition functions, using a TensorFlow reduce operation. From c9206cef6e851970bf8e095f8c787fd5b8b02468 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Mon, 18 Sep 2017 01:02:27 +0200 Subject: [PATCH 11/15] Restructuring exposed a bug --- gpflowopt/acquisition/acquisition.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index 5ebc75f..ce54797 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -305,7 +305,7 @@ def _setup(self): def build_acquisition(self, *args): pass - @Acquisition.models.getter + @property def models(self): return [model for acq in self.wrapped for model in acq.models] @@ -676,7 +676,7 @@ def _optimize_models(self): for idx, draw in enumerate(self.operands): draw.set_state(hypers[idx, :]) - @Acquisition.models.getter + @ParallelBatchAcquisition.models.getter def models(self): # Only return the models of the first operand, the copies remain hidden. return self.operands[0].models From 62ce738349114030f8aad7ea77d44c9ec8033a97 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Mon, 18 Sep 2017 09:58:40 +0200 Subject: [PATCH 12/15] Exluding abstract methods from the coverage report --- gpflowopt/acquisition/acquisition.py | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index ce54797..09c161e 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -55,7 +55,7 @@ def runnable(instance, *args, **kwargs): return runnable -class IAcquisition(Parameterized): +class IAcquisition(Parameterized): # pragma: no cover """ Interface for Acquisition functions mapping the belief represented by a Bayesian model into a score indicating how promising a point is for evaluation. @@ -302,7 +302,7 @@ def _setup(self): self._setup_objectives() @abstractmethod - def build_acquisition(self, *args): + def build_acquisition(self, *args): # pragma: no cover pass @property @@ -431,10 +431,6 @@ def data(self): def batch_size(self): return self._batch_size - @batch_size.setter - def batch_size(self, value): - self._batch_size = value - def constraint_indices(self): return np.empty((0,), dtype=int) @@ -532,7 +528,7 @@ def get_suggestion(self, optimizer): return result.x @abstractmethod - def build_acquisition(self, candidates): + def build_acquisition(self, candidates): # pragma: no cover pass @setup_required @@ -581,11 +577,6 @@ def __init__(self, operands, oper): def batch_size(self): return self[0].batch_size - @batch_size.setter - def batch_size(self, value): - for oper in self.operands: - oper.batch_size = value - @property def operands(self): return self.wrapped From e14f7f2a4fc256bade34d167342a2c4986094add Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Tue, 19 Sep 2017 00:21:02 +0200 Subject: [PATCH 13/15] Further testing of operators --- gpflowopt/acquisition/acquisition.py | 22 ++++++++--- testing/test_acquisition.py | 56 ++++++++++++++++++++++++---- 2 files changed, 65 insertions(+), 13 deletions(-) diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index 09c161e..b45bacd 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -463,7 +463,7 @@ def evaluate(self, Xcand): return self.build_acquisition(*tf.split(Xcand, num_or_size_splits=self.batch_size, axis=1)) @abstractmethod - def build_acquisition(self, *args): + def build_acquisition(self, *args): # pragma: no cover pass def __add__(self, other): @@ -499,6 +499,9 @@ def __radd__(self, other): if not isinstance(other, IAcquisition): return NotImplemented + if self.batch_size < other.batch_size: + return other + ParToSeqAcquisitionWrapper(self, other.batch_size) + if self.batch_size > other.batch_size: other = ParToSeqAcquisitionWrapper(other, self.batch_size) @@ -508,6 +511,9 @@ def __rmul__(self, other): if not isinstance(other, IAcquisition): return NotImplemented + if self.batch_size < other.batch_size: + return other * ParToSeqAcquisitionWrapper(self, other.batch_size) + if self.batch_size > other.batch_size: other = ParToSeqAcquisitionWrapper(other, self.batch_size) @@ -604,10 +610,13 @@ def __add__(self, other): else: return AcquisitionSum(self.operands.sorted_params + [other]) else: - super(AcquisitionSum, self).__add__(other) + return super(AcquisitionSum, self).__add__(other) def __radd__(self, other): - return AcquisitionSum([other] + self.operands.sorted_params) + if self.batch_size == other.batch_size: + return AcquisitionSum([other] + self.operands.sorted_params) + else: + return super(AcquisitionSum, self).__radd__(other) class AcquisitionProduct(AcquisitionAggregation): @@ -624,10 +633,13 @@ def __mul__(self, other): else: return AcquisitionProduct(self.operands.sorted_params + [other]) else: - super(AcquisitionSum, self).__add__(other) + return super(AcquisitionProduct, self).__mul__(other) def __rmul__(self, other): - return AcquisitionProduct([other] + self.operands.sorted_params) + if self.batch_size == other.batch_size: + return AcquisitionProduct([other] + self.operands.sorted_params) + else: + return super(AcquisitionProduct, self).__rmul__(other) class MCMCAcquistion(AcquisitionSum): diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index 11b5877..caf3120 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -422,41 +422,81 @@ def test_simple_hierarchy(self): @parameterized.expand([([SimpleAcquisition(create_parabola_model(domain)) for i in range(4)],), ([SimpleParallelBatch(create_parabola_model(domain)) for i in range(4)],)]) def test_multi_aggr_simple(self, acq): + # In these tests, the batch sizes align with self.test_session(): acq1, acq2, acq3, acq4 = acq joint = acq1 + acq2 + acq3 self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) - self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3]) + self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3]) joint = acq1 * acq2 * acq3 self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) - self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3]) + self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3]) first = acq2 + acq3 self.assertIsInstance(first, gpflowopt.acquisition.AcquisitionSum) - self.assertListEqual(first.wrapped.sorted_params, [acq2, acq3]) + self.assertListEqual(first.operands.sorted_params, [acq2, acq3]) joint = acq1 + first self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) - self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3]) + self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3]) first = acq2 * acq3 self.assertIsInstance(first, gpflowopt.acquisition.AcquisitionProduct) - self.assertListEqual(first.wrapped.sorted_params, [acq2, acq3]) + self.assertListEqual(first.operands.sorted_params, [acq2, acq3]) joint = acq1 * first self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) - self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3]) + self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3]) first = acq1 + acq2 second = acq3 + acq4 joint = first + second self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) - self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3, acq4]) + self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3, acq4]) first = acq1 * acq2 second = acq3 * acq4 joint = first * second self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) - self.assertListEqual(joint.wrapped.sorted_params, [acq1, acq2, acq3, acq4]) + self.assertListEqual(joint.operands.sorted_params, [acq1, acq2, acq3, acq4]) + + def test_multi_aggr_mixed(self): + # In these tests, batch sizes do not align + with self.test_session(): + acq1 = SimpleAcquisition(create_parabola_model(domain)) + acq2 = SimpleAcquisition(create_parabola_model(domain)) + acq3 = SimpleParallelBatch(create_parabola_model(domain)) + + first = acq1 + acq2 + joint = first + acq3 + self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) + self.assertTrue(len(joint.operands), 2) + self.assertEqual(joint.operands[1], acq3) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.operands[0].wrapped[0], first) + + joint = acq3 + first + self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionSum) + self.assertTrue(len(joint.operands), 2) + self.assertEqual(joint.operands[0], acq3) + self.assertTrue(isinstance(joint.operands[1], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.operands[1].wrapped[0], first) + + first = acq1 * acq2 + joint = first * acq3 + self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) + self.assertTrue(len(joint.operands), 2) + self.assertEqual(joint.operands[1], acq3) + self.assertTrue(isinstance(joint.operands[0], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.operands[0].wrapped[0], first) + + joint = acq3 * first + self.assertIsInstance(joint, gpflowopt.acquisition.AcquisitionProduct) + self.assertTrue(len(joint.operands), 2) + self.assertEqual(joint.operands[0], acq3) + self.assertTrue(isinstance(joint.operands[1], gpflowopt.acquisition.ParToSeqAcquisitionWrapper)) + self.assertEqual(joint.operands[1].wrapped[0], first) + + class TestRecompile(GPflowOptTestCase): From ffbe4990608cae865152206a1abbb5d6e0991c64 Mon Sep 17 00:00:00 2001 From: Joachim van der Herten Date: Thu, 9 Nov 2017 15:25:09 +0100 Subject: [PATCH 14/15] Fixed abstract classes in python 3 --- gpflowopt/acquisition/acquisition.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/gpflowopt/acquisition/acquisition.py b/gpflowopt/acquisition/acquisition.py index 3b027eb..d7f67fc 100644 --- a/gpflowopt/acquisition/acquisition.py +++ b/gpflowopt/acquisition/acquisition.py @@ -26,6 +26,7 @@ import copy from functools import wraps from abc import ABCMeta, abstractmethod, abstractproperty +import six float_type = settings.dtypes.float_type @@ -55,6 +56,7 @@ def runnable(instance, *args, **kwargs): return runnable +@six.add_metaclass(ABCMeta) class IAcquisition(Parameterized): # pragma: no cover """ Interface for Acquisition functions mapping the belief represented by a Bayesian model into a score indicating how @@ -72,8 +74,6 @@ class IAcquisition(Parameterized): # pragma: no cover for constrained optimization. The objects then form a tree hierarchy. """ - __metaclass__ = ABCMeta - @abstractmethod def _optimize_models(self): """ @@ -255,6 +255,7 @@ def __mul__(self, other): pass +@six.add_metaclass(ABCMeta) class AcquisitionWrapper(IAcquisition): """ Partial implementation of the Acquisition interface, useful for implementing acquisition functions which wrap @@ -270,8 +271,6 @@ class AcquisitionWrapper(IAcquisition): >>> return tf.constant(1.0) """ - __metaclass__ = ABCMeta - def __init__(self, wrap, *args, **kwargs): """ :param wrap: list of acquisition objects @@ -334,6 +333,7 @@ def __getitem__(self, item): return self.wrapped[item] +@six.add_metaclass(ABCMeta) class ParallelBatchAcquisition(IAcquisition): """ Abstract Acquisition class for batch acquisition functions, optimizing the batch in parallel (as opposed to @@ -349,8 +349,6 @@ class ParallelBatchAcquisition(IAcquisition): the objects handling objectives. """ - __metaclass__ = ABCMeta - def __init__(self, models=[], optimize_restarts=5, batch_size=1): """ :param models: list of GPflow models representing our beliefs about the problem @@ -520,6 +518,7 @@ def __rmul__(self, other): return AcquisitionProduct([other, self]) +@six.add_metaclass(ABCMeta) class Acquisition(ParallelBatchAcquisition): """ Class to be implemented for standard single-point acquisition functions like standard EI. @@ -527,8 +526,6 @@ class Acquisition(ParallelBatchAcquisition): This abstract class inherits the lazy setup mechanism from its ancestor. """ - __metaclass__ = ABCMeta - def get_suggestion(self, optimizer): result = optimizer.optimize(self._inverse_acquisition) return result.x From 4e749159c78596f90d3c134bb41e656847d87573 Mon Sep 17 00:00:00 2001 From: "nicolas.knudde@gmail.com" Date: Thu, 9 Nov 2017 17:31:34 +0100 Subject: [PATCH 15/15] Inital qEI implementation (doesn't work yet) --- gpflowopt/acquisition/__init__.py | 1 + gpflowopt/acquisition/qei.py | 85 +++++++++++++++++++++++++++++++ 2 files changed, 86 insertions(+) create mode 100644 gpflowopt/acquisition/qei.py diff --git a/gpflowopt/acquisition/__init__.py b/gpflowopt/acquisition/__init__.py index b9d2bb8..d34c132 100644 --- a/gpflowopt/acquisition/__init__.py +++ b/gpflowopt/acquisition/__init__.py @@ -25,6 +25,7 @@ from .hvpoi import HVProbabilityOfImprovement # Batch +from .qei import qExpectedImprovement # Black-box constraint from .pof import ProbabilityOfFeasibility diff --git a/gpflowopt/acquisition/qei.py b/gpflowopt/acquisition/qei.py new file mode 100644 index 0000000..a6560cb --- /dev/null +++ b/gpflowopt/acquisition/qei.py @@ -0,0 +1,85 @@ +# Copyright 2017 Joachim van der Herten +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from .acquisition import ParallelBatchAcquisition + +from gpflow.model import Model +from gpflow.param import DataHolder +from gpflow import settings + +import numpy as np +import tensorflow as tf +import tensorflow.contrib.distributions as ds + +stability = settings.numerics.jitter_level +float_type = settings.dtypes.float_type + + +class qExpectedImprovement(ParallelBatchAcquisition): + def __init__(self, model, batch_size=4): + """ + :param model: GPflow model (single output) representing our belief of the objective + """ + super(qExpectedImprovement, self).__init__(model, batch_size=batch_size) + self.fmin = DataHolder(np.zeros(1)) + self._setup() + + def _setup(self): + super(qExpectedImprovement, 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) + self.fmin.set_data(np.min(samples_mean, axis=0)) + + def build_acquisition(self, *args): + # Obtain predictive distributions for candidates + N, D = tf.shape(args[0])[0], tf.shape(args[0])[1] + q = self.batch_size + + Xcand = tf.transpose(tf.stack(args, axis=0), perm=[1, 0, 2]) + m, sig = tf.map_fn(lambda x: self.models[0].build_predict(x, full_cov=True), Xcand, + dtype=(float_type, float_type)) # N x q x 1, N x q x q + + eye = tf.tile(tf.expand_dims(tf.eye(q, dtype=float_type), 0), [q, 1, 1]) + A = eye + A = A - tf.transpose(eye, perm=[1, 2, 0]) + A = A - tf.transpose(eye, perm=[2, 0, 1]) # q x q x q (k x q x q) + + mk = tf.tensordot(A, m, [[2], [1]]) # N x q(k) x q Mean of Zk (k x q) + sigk = tf.tensordot(A, sig, [[2], [1]]) # N x q x q x q + sigk = tf.reduce_sum(tf.expand_dims(sigk, 3) * tf.expand_dims(tf.expand_dims(A, 0), 2), axis=-1)# N x q(k) x q x q + + a = tf.tile(tf.expand_dims(tf.eye(self.batch_size, self.batch_size, dtype=float_type), 0), [self.batch_size, 1, 1]) + A1 = tf.gather_nd(a, [[[i, j + (j >= i)] for j in range(self.batch_size)] for i in range(self.batch_size)]) + # q(i) x (q-1) x q + + bk = -self.fmin * tf.eye(q, dtype=float_type) # q(k) x q + + Sigk_t = tf.expand_dims(tf.transpose(tf.tensordot(sigk, A1, [[-2], [2]]), [0, 1, 3, 4, 2]), -2) # N x q(k) x q(i) x (q-1) x 1 x q + Sigk = tf.reduce_sum(tf.expand_dims(tf.expand_dims(tf.expand_dims(A1, 0), 0), -3)* Sigk_t, axis=-1) + #Sigk = tf.einsum('ijklm,lnk->ijlmn', Sigk_t, A1) # N x q(k) x q(i) x q-1 x q-1 + c = tf.tensordot(tf.expand_dims(bk, 0) - mk, A1, [[2], [2]]) + c = tf.tensordot(tf.expand_dims(bk, 0) - mk, A1, [[2], [2]]) # N x q(k) x q(i) x q-1 + F = tf.expand_dims(tf.expand_dims(tf.expand_dims(bk, 0) - mk, -1) / tf.matrix_diag_part(sigk), + -1) # N x q(k) x q(i) x 1 + F *= tf.transpose(tf.squeeze(tf.matrix_diag_part(tf.transpose(Sigk_t, [0, 1, 3, 4, 5, 2])), 4), [0, 1, 3, 2]) + c -= F + + MVN = ds.MultivariateNormalFullCovariance(loc=mk, covariance_matrix=sigk) + MVN2 = ds.MultivariateNormalFullCovariance(loc=tf.zeros(tf.shape(c), dtype=float_type), covariance_matrix=Sigk) + UVN = ds.MultivariateNormalDiag(loc=mk, scale_diag=tf.sqrt(tf.matrix_diag_part(sigk))) + t1 = tf.reduce_sum((self.fmin - m) * MVN.cdf(bk), axis=1) + sigkk = tf.transpose(tf.matrix_diag_part(tf.transpose(sigk, perm=[0, 3, 1, 2])), perm=[0, 2, 1]) + t2 = tf.reduce_sum(sigkk * UVN.pdf(-tf.expand_dims(bk, 0)) * MVN2.cdf(c), axis=[1, 2]) + return tf.add(t1, t2, name=self.__class__.__name__)