Skip to content

Optimization tmp fix #5

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
javdrher merged 5 commits into from
May 5, 2017
Merged

Optimization tmp fix #5

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
70e4d8a
Adjusted model initial hyperparameters in the tests. Also added funct…
javdrher May 5, 2017
58a41d3
Fixed lengthscales for python 2.7
javdrher May 5, 2017
95d5abe
Prior should solve the optimization instabilities
javdrher May 5, 2017
dd5aea4
Implemented a restart strategy for optimization of the hypers
javdrher May 5, 2017
d69c046
Fixed initialization error of Gamma prior for python 2.7
javdrher May 5, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 19 additions & 6 deletions GPflowOpt/acquisition.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
71 changes: 41 additions & 30 deletions testing/test_acquisition.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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):
Expand All @@ -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
Expand All @@ -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):
Expand All @@ -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")
Expand All @@ -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):
Expand All @@ -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):
Expand All @@ -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)])
Expand All @@ -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))
Expand All @@ -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))
3 changes: 2 additions & 1 deletion testing/test_optimizers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down