Merge pull request #121 from akangasr/gpy_fix

Improve consistency of GPyModel fitting, add tests

elfi/bo/gpy_model.py

@@ -61,7 +61,7 @@ def __init__(self, input_dim=1, bounds=None, kernel=None,
             self.bounds = [(0,1)] * self.input_dim
         if len(self.bounds) != self.input_dim:
             raise ValueError("Number of bounds should match input dimension. " +
-                    "Expected {}. Received {}.".format(self.input_dim, len(self_bounds)))
+                    "Expected {}. Received {}.".format(self.input_dim, len(self.bounds)))
         self.noise_var = noise_var
         self.optimizer = optimizer
         self.max_opt_iters = max_opt_iters
@@ -128,29 +128,33 @@ def set_kernel(self, kernel=None, kernel_class=None, kernel_var=None,
         ----------
         see constructor
         """
-        if kernel is not None:
+        self.kernel = kernel
+        self.kernel_class = kernel_class
+        self.kernel_var = kernel_var
+        self.kernel_scale = kernel_scale
+        if self.kernel is not None:
             # explicit kernel supplied
             if kernel.input_dim != self.input_dim:
                 raise ValueError("Kernel input_dim must match model input_dim.")
-            self.kernel = kernel
-        else:
-            self.kernel_class = kernel_class or self.kernel_class
-            self.kernel_var = kernel_var or self.kernel_var
-            self.kernel_scale = kernel_scale or self.kernel_scale
-            if isinstance(self.kernel_class, str):
-                self.kernel_class = getattr(GPy.kern, self.kernel_class)
-            self.kernel = self.kernel_class(input_dim=self.input_dim,
-                                            variance=self.kernel_var,
-                                            lengthscale=self.kernel_scale)
         if self.gp is not None:
             # re-fit gp with new kernel
             self._fit_gp(self.gp.X, self.gp.Y)
 
+    def get_kernel(self):
+        if self.kernel is not None:
+            return self.kernel.copy()
+        else:
+            if isinstance(self.kernel_class, str):
+                self.kernel_class = getattr(GPy.kern, self.kernel_class)
+            return self.kernel_class(input_dim=self.input_dim,
+                                     variance=self.kernel_var,
+                                     lengthscale=self.kernel_scale)
+
     def _fit_gp(self, X, Y):
         """Constructs the gp model.
         """
         self.gp = GPy.models.GPRegression(X=X, Y=Y,
-                                          kernel=self.kernel,
+                                          kernel=self.get_kernel(),
                                           noise_var=self.noise_var)
 
         # FIXME: move to initialization
@@ -247,11 +251,14 @@ def n_observations(self):
     def copy(self):
         model = GPyModel(input_dim=self.input_dim,
                          bounds=self.bounds[:],
-                         kernel=self.kernel.copy(),
+                         kernel=self.kernel,
+                         kernel_class=self.kernel_class,
+                         kernel_var=self.kernel_var,
+                         kernel_scale=self.kernel_scale,
                          noise_var=self.noise_var,
                          optimizer=self.optimizer,
                          max_opt_iters=self.max_opt_iters)
         if self.gp is not None:
-            model._fit_gp(self.gp.X[:], self.gp.Y[:])
+            model.update(self.gp.X[:], self.gp.Y[:])
         return model
 

tests/unit/bo/test_gpy_model.py

@@ -52,6 +52,42 @@ def test_four_2d_samples(self):
         assert abs(pred1[0] + pred2[0]) < 1e-3
         np.testing.assert_allclose(pred1[1:2], pred2[1:2], atol=1e-3)
 
+    def test_copy_of_gp_gives_same_results(self):
+        n = 10
+        X = np.random.uniform(size=(n, 2))
+        Y = np.random.uniform(size=(n, 1))
+        gp1 = GPyModel(input_dim=2, bounds=((0,1),(0,1)))
+        gp1.update(X, Y)
+        gp2 = gp1.copy()
+        loc = np.random.uniform(size=(n, 2))
+        for i in range(n):
+            m, s, s2 = gp1.evaluate(loc[i][None,:])
+            m_, s_, s2_ = gp2.evaluate(loc[i][None,:])
+            assert np.abs(m - m_) < 1e-5
+            assert np.abs(s - s_) < 1e-5
+            assert np.abs(s2 - s2_) < 1e-5
+
+    def test_update_order_is_irrelevant_for_end_result(self):
+        n = 10
+        X = np.random.uniform(size=(n, 2))
+        Y = np.random.uniform(size=(n, 1))
+        order1 = np.random.permutation(n)
+        order2 = np.random.permutation(n)
+        gp1 = GPyModel(input_dim=2, bounds=((0,1),(0,1)))
+        gp2 = GPyModel(input_dim=2, bounds=((0,1),(0,1)))
+        for i in order1:
+            gp1.update(X[i][None,:], Y[i][None,:])
+        for i in order2:
+            gp2.update(X[i][None,:], Y[i][None,:])
+        loc = np.random.uniform(size=(n, 2))
+        for i in range(n):
+            m, s, s2 = gp1.evaluate(loc[i][None,:])
+            m_, s_, s2_ = gp2.evaluate(loc[i][None,:])
+            assert np.abs(m - m_) < 1e-5
+            assert np.abs(s - s_) < 1e-5
+            assert np.abs(s2 - s2_) < 1e-5
+
+
     # FIXME
     # def test_change_kernel(self):
     #     bounds = ((0, 1), )