diff --git a/GPflowOpt/acquisition.py b/GPflowOpt/acquisition.py index c24e553..adcd209 100644 --- a/GPflowOpt/acquisition.py +++ b/GPflowOpt/acquisition.py @@ -37,17 +37,30 @@ class Acquisition(Parameterized): optimization) """ - def __init__(self, models=[]): + def __init__(self, models=[], optimize_restarts=5): super(Acquisition, self).__init__() self.models = ParamList(np.atleast_1d(models).tolist()) + self._default_params = map(lambda m: m.get_free_state(), self.models) + + assert(optimize_restarts >= 0) + self._optimize_restarts = optimize_restarts self._optimize_all() def _optimize_all(self): - for model in self.models: - # If likelihood variance is close to zero, updating data may result in non-invertible K - # Increase likelihood variance a bit. - model.likelihood.variance = 4.0 - model.optimize() + for model, hypers in zip(self.models, self._default_params): + runs = [] + # Start from supplied hyperparameters + model.set_state(hypers) + for i in range(self._optimize_restarts): + if i > 0: + model.randomize() + try: + result = model.optimize() + runs.append(result) + except tf.errors.InvalidArgumentError: + print("Warning - optimization restart {0}/{1} failed".format(i + 1, self._optimize_restarts)) + best_idx = np.argmin(map(lambda r: r.fun, runs)) + model.set_state(runs[best_idx].x) def _build_acquisition_wrapper(self, Xcand, gradients=True): acq = self.build_acquisition(Xcand) diff --git a/testing/test_acquisition.py b/testing/test_acquisition.py index 651fee2..e6afb6d 100644 --- a/testing/test_acquisition.py +++ b/testing/test_acquisition.py @@ -22,6 +22,21 @@ class _TestAcquisition(object): def domain(self): return np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)]) + def create_parabola_model(self, design=None): + if design is None: + design = GPflowOpt.design.FactorialDesign(4, self.domain) + X, Y = design.generate(), parabola2d(design.generate()) + m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.var(axis=0))) + m.kern.variance.prior = GPflow.priors.Gamma(3.0, 1.0 / 3.0) + return m + + def create_plane_model(self, design=None): + if design is None: + design = GPflowOpt.design.FactorialDesign(5, self.domain) + X, Y = design.generate(), plane(design.generate()) + m = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + return m + def setUp(self): self.acquisition = None @@ -73,9 +88,7 @@ def test_data_update(self): class TestExpectedImprovement(_TestAcquisition, unittest.TestCase): def setUp(self): super(TestExpectedImprovement, self).setUp() - design = GPflowOpt.design.FactorialDesign(4, self.domain) - X, Y = design.generate(), parabola2d(design.generate()) - self.model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) + self.model = self.create_parabola_model() self.acquisition = GPflowOpt.acquisition.ExpectedImprovement(self.model) def test_object_integrity(self): @@ -93,7 +106,7 @@ def test_setup(self): p = np.array([[0.0, 0.0]]) self.acquisition.set_data(np.vstack((self.model.X.value, p)), np.vstack((self.model.Y.value, parabola2d(p)))) - self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-2), msg="fmin not updated") + self.assertTrue(np.allclose(self.acquisition.fmin.value, 0, atol=1e-1), msg="fmin not updated") def test_EI_validity(self): Xcenter = np.random.rand(20, 2) * 0.25 - 0.125 @@ -109,9 +122,7 @@ def test_EI_validity(self): class TestProbabilityOfFeasibility(_TestAcquisition, unittest.TestCase): def setUp(self): super(TestProbabilityOfFeasibility, self).setUp() - design = GPflowOpt.design.FactorialDesign(5, self.domain) - X, Y = design.generate(), plane(design.generate()) - self.model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) + self.model = self.create_plane_model() self.acquisition = GPflowOpt.acquisition.ProbabilityOfFeasibility(self.model) def test_object_integrity(self): @@ -128,7 +139,6 @@ def test_PoF_validity(self): class _TestAcquisitionBinaryOperator(_TestAcquisition): - def test_object_integrity(self): self.assertTrue(isinstance(self.acquisition.lhs, GPflowOpt.acquisition.Acquisition), msg="Left hand operand should be an acquisition object") @@ -150,22 +160,23 @@ def test_data(self): class TestAcquisitionSum(_TestAcquisitionBinaryOperator, unittest.TestCase): def setUp(self): super(TestAcquisitionSum, self).setUp() - design = GPflowOpt.design.FactorialDesign(4, self.domain) - X, Y = design.generate(), parabola2d(design.generate()) - self.model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) - self.acquisition = GPflowOpt.acquisition.AcquisitionSum(GPflowOpt.acquisition.ExpectedImprovement(self.model), - GPflowOpt.acquisition.ExpectedImprovement(self.model)) + self.models = [self.create_parabola_model(), self.create_parabola_model()] + self.acquisition = GPflowOpt.acquisition.AcquisitionSum( + GPflowOpt.acquisition.ExpectedImprovement(self.models[0]), + GPflowOpt.acquisition.ExpectedImprovement(self.models[1])) def test_sum_validity(self): - single_ei = GPflowOpt.acquisition.ExpectedImprovement(self.model) design = GPflowOpt.design.FactorialDesign(4, self.domain) + m = self.create_parabola_model(design) + single_ei = GPflowOpt.acquisition.ExpectedImprovement(m) p1 = self.acquisition.evaluate(design.generate()) p2 = single_ei.evaluate(design.generate()) np.testing.assert_allclose(p2, p1 / 2, rtol=1e-3, err_msg="The sum of 2 EI should be the double of only EI") def test_generating_operator(self): - joint = GPflowOpt.acquisition.ExpectedImprovement(self.model) + \ - GPflowOpt.acquisition.ExpectedImprovement(self.model) + design = GPflowOpt.design.FactorialDesign(4, self.domain) + joint = GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model(design)) + \ + GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model(design)) self.assertTrue(isinstance(joint, GPflowOpt.acquisition.AcquisitionSum)) def test_indices(self): @@ -178,24 +189,23 @@ def test_indices(self): class TestAcquisitionProduct(_TestAcquisitionBinaryOperator, unittest.TestCase): def setUp(self): super(TestAcquisitionProduct, self).setUp() - design = GPflowOpt.design.FactorialDesign(4, self.domain) - X, Y = design.generate(), parabola2d(design.generate()) - self.model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) + self.models = [self.create_parabola_model(), self.create_parabola_model()] self.acquisition = GPflowOpt.acquisition.AcquisitionProduct( - GPflowOpt.acquisition.ExpectedImprovement(self.model), - GPflowOpt.acquisition.ExpectedImprovement(self.model)) + GPflowOpt.acquisition.ExpectedImprovement(self.models[0]), + GPflowOpt.acquisition.ExpectedImprovement(self.models[1])) def test_product_validity(self): - single_ei = GPflowOpt.acquisition.ExpectedImprovement(self.model) design = GPflowOpt.design.FactorialDesign(4, self.domain) + m = self.create_parabola_model(design) + single_ei = GPflowOpt.acquisition.ExpectedImprovement(m) p1 = self.acquisition.evaluate(design.generate()) p2 = single_ei.evaluate(design.generate()) np.testing.assert_allclose(p2, np.sqrt(p1), rtol=1e-3, err_msg="The product of 2 EI should be the square of one EI") def test_generating_operator(self): - joint = GPflowOpt.acquisition.ExpectedImprovement(self.model) * \ - GPflowOpt.acquisition.ExpectedImprovement(self.model) + joint = GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model()) * \ + GPflowOpt.acquisition.ExpectedImprovement(self.create_parabola_model()) self.assertTrue(isinstance(joint, GPflowOpt.acquisition.AcquisitionProduct)) def test_indices(self): @@ -206,7 +216,6 @@ def test_indices(self): class TestJointAcquisition(unittest.TestCase): - @property def domain(self): return np.sum([GPflowOpt.domain.ContinuousParameter("x{0}".format(i), -1, 1) for i in range(1, 3)]) @@ -216,8 +225,8 @@ def test_constrained_EI(self): X = design.generate() Yo = parabola2d(X) Yc = plane(X) - m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True)) - m2 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True)) + m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m2 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0) / 2)) joint = GPflowOpt.acquisition.ExpectedImprovement(m1) * GPflowOpt.acquisition.ProbabilityOfFeasibility(m2) np.testing.assert_allclose(joint.objective_indices(), np.array([0], dtype=int)) @@ -228,10 +237,12 @@ def test_hierarchy(self): X = design.generate() Yo = parabola2d(X) Yc = plane(X) - m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True)) - m2 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True)) + m1 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m2 = GPflow.gpr.GPR(X, Yo, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0))) + m3 = GPflow.gpr.GPR(X, Yc, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.std(axis=0) / 2)) joint = GPflowOpt.acquisition.ExpectedImprovement(m1) * \ - (GPflowOpt.acquisition.ProbabilityOfFeasibility(m2) + GPflowOpt.acquisition.ExpectedImprovement(m1)) + (GPflowOpt.acquisition.ProbabilityOfFeasibility(m3) + + GPflowOpt.acquisition.ExpectedImprovement(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)) diff --git a/testing/test_optimizers.py b/testing/test_optimizers.py index 0b76eba..5f6cdf5 100644 --- a/testing/test_optimizers.py +++ b/testing/test_optimizers.py @@ -134,7 +134,8 @@ def setUp(self): super(TestBayesianOptimizer, self).setUp() design = GPflowOpt.design.FactorialDesign(4, self.domain) X, Y = design.generate(), parabola2d(design.generate())[0] - model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True)) + model = GPflow.gpr.GPR(X, Y, GPflow.kernels.RBF(2, ARD=True, lengthscales=X.var(axis=0))) + model.kern.variance.prior = GPflow.priors.Gamma(3, 1.0 / 3.0) acquisition = GPflowOpt.acquisition.ExpectedImprovement(model) self.optimizer = GPflowOpt.BayesianOptimizer(self.domain, acquisition)