Adds a couple of tests to make sure optimizer acq_max is not completely off. And finishes a couple of things in the migration to the new GP objects.

Author: fmfn

bayes_opt/helpers.py

@@ -86,27 +86,19 @@ def utility(self, x, gp, y_max):
 
     @staticmethod
     def _ucb(x, gp, kappa):
-        mean, var = gp.predict(x, return_cov=True)
-        return mean + kappa * np.sqrt(var)
+        mean, std = gp.predict(x, return_std=True)
+        return mean + kappa * std
 
     @staticmethod
     def _ei(x, gp, y_max, xi):
-        mean, var = gp.predict(x, return_cov=True)
-
-        # Avoid points with zero variance
-        var = np.maximum(var, 1e-9 + 0 * var)
-
-        z = (mean - y_max - xi)/np.sqrt(var)
-        return (mean - y_max - xi) * norm.cdf(z) + np.sqrt(var) * norm.pdf(z)
+        mean, std = gp.predict(x, return_std=True)
+        z = (mean - y_max - xi)/std
+        return (mean - y_max - xi) * norm.cdf(z) + std * norm.pdf(z)
 
     @staticmethod
     def _poi(x, gp, y_max, xi):
-        mean, var = gp.predict(x, return_cov=True)
-
-        # Avoid points with zero variance
-        var = np.maximum(var, 1e-9 + 0 * var)
-
-        z = (mean - y_max - xi)/np.sqrt(var)
+        mean, std = gp.predict(x, return_std=True)
+        z = (mean - y_max - xi)/std
         return norm.cdf(z)
 
 

setup.py

@@ -2,13 +2,13 @@
 
 setup(
     name='bayesian-optimization',
-    version='0.1.0',
+    version='0.3.0',
     url='https://github.com/fmfn/BayesianOptimization',
     packages=find_packages(),
     description='Bayesian Optimization package',
     install_requires=[
         "numpy >= 1.9.0",
         "scipy >= 0.14.0",
-        "scikit-learn >= 0.16.1",
+        "scikit-learn >= 0.18.0",
     ],
 )

tests/tests_helper_functions.py

@@ -1,28 +1,82 @@
 from __future__ import print_function
 from __future__ import division
 import unittest
+import numpy as np
 from sklearn.gaussian_process import GaussianProcessRegressor
 from sklearn.gaussian_process.kernels import Matern
+import sys
+sys.path.append("./")
+from bayes_opt.helpers import UtilityFunction, acq_max
+
+
+def get_globals():
+    X = np.array([
+        [0.00, 0.00],
+        [0.99, 0.99],
+        [0.00, 0.99],
+        [0.99, 0.00],
+        [0.50, 0.50],
+        [0.25, 0.50],
+        [0.50, 0.25],
+        [0.75, 0.50],
+        [0.50, 0.75],
+    ])
+
+    def get_y(X):
+        return -(X[:, 0] - 0.3) ** 2 - 0.5 * (X[:, 1] - 0.6)**2 + 2
+    y = get_y(X)
+
+    mesh = np.dstack(
+        np.meshgrid(np.arange(0, 1, 0.01), np.arange(0, 1, 0.01))
+    ).reshape(-1, 2)
+
+    GP = GaussianProcessRegressor(
+        kernel=Matern(),
+        n_restarts_optimizer=25,
+    )
+    GP.fit(X, y)
+
+    return {'x': X, 'y': y, 'gp': GP, 'mesh': mesh}
+
+
+def brute_force_maximum(MESH, GP, kind='ucb', kappa=1.0, xi=1e-6):
+    uf = UtilityFunction(kind=kind, kappa=kappa, xi=xi)
+
+    mesh_vals = uf.utility(MESH, GP, 2)
+    max_val = mesh_vals.max()
+    max_arg_val = MESH[np.argmax(mesh_vals)]
+
+    return max_val, max_arg_val
+
+
+GLOB = get_globals()
+X, Y, GP, MESH = GLOB['x'], GLOB['y'], GLOB['gp'], GLOB['mesh']
 
 
 class TestMaximizationOfAcquisitionFunction(unittest.TestCase):
 
-    def setUp(self, msg=""):
-        self.gp = GaussianProcessRegressor(
-            kernel=Matern(),
-            n_restarts_optimizer=25,
+    def setUp(self, kind='ucb', kappa=1.0, xi=1e-6):
+        self.util = UtilityFunction(kind=kind, kappa=kappa, xi=xi)
+        self.episilon = 1e-2
+        self.y_max = 2.0
+
+    def test_acq_max_function_with_ucb_algo(self):
+        self.setUp(kind='ucb', kappa=1.0, xi=1.0)
+        max_arg = acq_max(
+            self.util.utility, GP, self.y_max, bounds=np.array([[0, 1], [0, 1]])
         )
+        _, brute_max_arg = brute_force_maximum(MESH, GP)
 
-        X = np.array([
-            [1, 2],
-            [2, 3],
-            [0, -1],
-        ])
+        self.assertTrue( all(abs(brute_max_arg - max_arg) < self.episilon))
 
-        y = X[:, 0]**2 - 2 * X[:, 1]
+    def test_ei_max_function_with_ucb_algo(self):
+        self.setUp(kind='ei', kappa=1.0, xi=1e-6)
+        max_arg = acq_max(
+            self.util.utility, GP, self.y_max, bounds=np.array([[0, 1], [0, 1]])
+        )
+        _, brute_max_arg = brute_force_maximum(MESH, GP, kind='ei')
 
-    def test_something(self):
-        pass
+        self.assertTrue( all(abs(brute_max_arg - max_arg) < self.episilon))
 
 
 if __name__ == '__main__':