Merge branch 'devel' of github.com:SheffieldML/GPy into devel

GPy/core/gp.py

@@ -108,9 +108,15 @@ class GP(Model):
 
         # The predictive variable to be used to predict using the posterior object's
         # woodbury_vector and woodbury_inv is defined as predictive_variable
+        # as long as the posterior has the right woodbury entries.
+        # It is the input variable used for the covariance between
+        # X_star and the posterior of the GP.
         # This is usually just a link to self.X (full GP) or self.Z (sparse GP).
         # Make sure to name this variable and the predict functions will "just work"
-        # as long as the posterior has the right woodbury entries.
+        # In maths the predictive variable is:
+        #         K_{xx} - K_{xp}W_{pp}^{-1}K_{px}
+        #         W_{pp} := \texttt{Woodbury inv}
+        #         p := _predictive_variable
         self._predictive_variable = self.X
 
 
@@ -213,7 +219,7 @@ class GP(Model):
             Kxx = kern.K(Xnew)
             if self.posterior.woodbury_inv.ndim == 2:
                 var = Kxx - np.dot(Kx.T, np.dot(self.posterior.woodbury_inv, Kx))
-            elif self.posterior.woodbury_inv.ndim == 3:
+            elif self.posterior.woodbury_inv.ndim == 3: # Missing data
                 var = np.empty((Kxx.shape[0],Kxx.shape[1],self.posterior.woodbury_inv.shape[2]))
                 from ..util.linalg import mdot
                 for i in range(var.shape[2]):
@@ -223,7 +229,7 @@ class GP(Model):
             Kxx = kern.Kdiag(Xnew)
             if self.posterior.woodbury_inv.ndim == 2:
                 var = (Kxx - np.sum(np.dot(self.posterior.woodbury_inv.T, Kx) * Kx, 0))[:,None]
-            elif self.posterior.woodbury_inv.ndim == 3:
+            elif self.posterior.woodbury_inv.ndim == 3: # Missing data
                 var = np.empty((Kxx.shape[0],self.posterior.woodbury_inv.shape[2]))
                 for i in range(var.shape[1]):
                     var[:, i] = (Kxx - (np.sum(np.dot(self.posterior.woodbury_inv[:, :, i].T, Kx) * Kx, 0)))
@@ -364,11 +370,15 @@ class GP(Model):
                 var_jac = dK2_dXdX - np.einsum('qim,miq->iq', dK_dXnew_full.T.dot(wi), dK_dXnew_full)
             return var_jac
 
-        if self.posterior.woodbury_inv.ndim == 3:
-            var_jac = []
-            for d in range(self.posterior.woodbury_inv.shape[2]):
-                var_jac.append(compute_cov_inner(self.posterior.woodbury_inv[:, :, d]))
-            var_jac = np.concatenate(var_jac)
+        if self.posterior.woodbury_inv.ndim == 3: # Missing data:
+            if full_cov:
+                var_jac = np.empty((Xnew.shape[0],Xnew.shape[0],Xnew.shape[1],self.output_dim))
+                for d in range(self.posterior.woodbury_inv.shape[2]):
+                    var_jac[:, :, :, d] = compute_cov_inner(self.posterior.woodbury_inv[:, :, d])
+            else:
+                var_jac = np.empty((Xnew.shape[0],Xnew.shape[1],self.output_dim))
+                for d in range(self.posterior.woodbury_inv.shape[2]):
+                    var_jac[:, :, d] = compute_cov_inner(self.posterior.woodbury_inv[:, :, d])
         else:
             var_jac = compute_cov_inner(self.posterior.woodbury_inv)
         return mean_jac, var_jac
@@ -391,10 +401,11 @@ class GP(Model):
 
         mu_jac, var_jac = self.predict_jacobian(Xnew, kern, full_cov=False)
         mumuT = np.einsum('iqd,ipd->iqp', mu_jac, mu_jac)
+        Sigma = np.zeros(mumuT.shape)
         if var_jac.ndim == 3:
-            Sigma = np.einsum('iqd,ipd->iqp', var_jac, var_jac)
+            Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = var_jac.sum(-1)
         else:
-            Sigma = self.output_dim*np.einsum('iq,ip->iqp', var_jac, var_jac)
+            Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = self.output_dim*var_jac
         G = 0.
         if mean:
             G += mumuT
@@ -412,8 +423,13 @@ class GP(Model):
         """
         G = self.predict_wishard_embedding(Xnew, kern, mean, covariance)
         from ..util.linalg import jitchol
-        return np.array([np.sqrt(np.exp(2*np.sum(np.log(np.diag(jitchol(G[n, :, :])))))) for n in range(Xnew.shape[0])])
-        #return np.array([np.sqrt(np.linalg.det(G[n, :, :])) for n in range(Xnew.shape[0])])
+        mag = np.empty(Xnew.shape[0])
+        for n in range(Xnew.shape[0]):
+            try:
+                mag[n] = np.sqrt(np.exp(2*np.sum(np.log(np.diag(jitchol(G[n, :, :]))))))
+            except:
+                mag[n] = np.sqrt(np.linalg.det(G[n, :, :]))
+        return mag
 
     def posterior_samples_f(self,X,size=10, full_cov=True):
         """

GPy/likelihoods/link_functions.py

@@ -141,7 +141,7 @@ class Log_ex_1(GPTransformation):
 
     """
     def transf(self,f):
-        return scipy.log1p(safe_exp(f))
+        return scipy.special.log1p(safe_exp(f))
 
     def dtransf_df(self,f):
         ef = safe_exp(f)

GPy/models/gplvm.py

@@ -36,8 +36,10 @@ class GPLVM(GP):
         likelihood = Gaussian()
 
         super(GPLVM, self).__init__(X, Y, kernel, likelihood, name='GPLVM')
+
         self.X = Param('latent_mean', X)
         self.link_parameter(self.X, index=0)
+        self._predictive_variable = self.X
 
     def parameters_changed(self):
         super(GPLVM, self).parameters_changed()

GPy/plotting/matplot_dep/dim_reduction_plots.py

@@ -119,7 +119,7 @@ def plot_latent(model, labels=None, which_indices=None,
         Xtest_full[:, [input_1, input_2]] = x
         _, var = model.predict(Xtest_full, **predict_kwargs)
         var = var[:, :1]
-        return np.log(var)
+        return 2*np.sqrt(var)
 
     #Create an IMshow controller that can re-plot the latent space shading at a good resolution
     if plot_limits is None:

GPy/testing/link_function_tests.py

@@ -1,5 +1,5 @@
 import numpy as np
-import scipy as sp
+import scipy
 from scipy.special import cbrt
 from GPy.models import GradientChecker
 _lim_val = np.finfo(np.float64).max
@@ -92,18 +92,18 @@ class LinkFunctionTests(np.testing.TestCase):
         link = Log_ex_1()
         lim_of_inf = _lim_val_exp
 
-        np.testing.assert_almost_equal(np.log1p(np.exp(self.mid_f)), link.transf(self.mid_f))
-        assert np.isinf(np.log1p(np.exp(np.log(self.f_upper_lim))))
+        np.testing.assert_almost_equal(scipy.special.log1p(np.exp(self.mid_f)), link.transf(self.mid_f))
+        assert np.isinf(scipy.special.log1p(np.exp(np.log(self.f_upper_lim))))
         #Check the clipping works
         np.testing.assert_almost_equal(link.transf(self.f_lower_lim), 0, decimal=5)
         #Need to look at most significant figures here rather than the decimals
-        np.testing.assert_approx_equal(link.transf(self.f_upper_lim), np.log1p(_lim_val), significant=5)
+        np.testing.assert_approx_equal(link.transf(self.f_upper_lim), scipy.special.log1p(_lim_val), significant=5)
         self.check_overflow(link, lim_of_inf)
 
         #Check that it would otherwise fail
         beyond_lim_of_inf = lim_of_inf + 10.0
         old_err_state = np.seterr(over='ignore')
-        self.assertTrue(np.isinf(np.log1p(np.exp(beyond_lim_of_inf))))
+        self.assertTrue(np.isinf(scipy.special.log1p(np.exp(beyond_lim_of_inf))))
         np.seterr(**old_err_state)
 
 

GPy/testing/run_coverage.sh

@@ -1 +1 @@
-nosetests . --with-coverage --logging-level=INFO --cover-html --cover-html-dir=coverage --cover-package=GPy --cover-erase
+nosetests . --with-coverage --logging-level=INFO --cover-html --cover-html-dir=coverage --cover-package=GPy --cover-erase --cover-omit=GPy.examples