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

GPy/core/parameterization/variational.py

@@ -135,4 +135,4 @@ class SpikeAndSlabPosterior(VariationalPosterior):
         import sys
         assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
         from ...plotting.matplot_dep import variational_plots
-        return variational_plots.plot(self,*args)
+        return variational_plots.plot_SpikeSlab(self,*args)

GPy/examples/dimensionality_reduction.py

@@ -515,3 +515,28 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose
         lvm_visualizer.close()
 
     return m
+
+def ssgplvm_simulation_linear():
+    import numpy as np
+    import GPy
+    N, D, Q = 1000, 20, 5
+    pi = 0.2
+    
+    def sample_X(Q, pi):
+        x = np.empty(Q)
+        dies = np.random.rand(Q)
+        for q in xrange(Q):
+            if dies[q]<pi:
+                x[q] = np.random.randn()
+            else:
+                x[q] = 0.
+        return x
+    
+    Y = np.empty((N,D))
+    X = np.empty((N,Q))
+    # Generate data from random sampled weight matrices
+    for n in xrange(N):
+        X[n] = sample_X(Q,pi)
+        w = np.random.randn(D,Q)
+        Y[n] = np.dot(w,X[n])
+    

GPy/examples/regression.py

@@ -284,7 +284,7 @@ def toy_poisson_rbf_1d_laplace(optimize=True, plot=True):
 
     kern = GPy.kern.RBF(1)
     poisson_lik = GPy.likelihoods.Poisson()
-    laplace_inf = GPy.inference.latent_function_inference.LaplaceInference()
+    laplace_inf = GPy.inference.latent_function_inference.Laplace()
 
     # create simple GP Model
     m = GPy.core.GP(X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf)

GPy/inference/latent_function_inference/var_dtc.py

@@ -80,7 +80,7 @@ class VarDTC(object):
         
         Kmm = kern.K(Z) +np.eye(Z.shape[0]) * self.const_jitter
 
-        Lm = jitchol(Kmm)
+        Lm = jitchol(Kmm+np.eye(Z.shape[0])*self.const_jitter)
 
         # The rather complex computations of A
         if uncertain_inputs:

GPy/kern/_src/linear.py

@@ -6,10 +6,12 @@ import numpy as np
 from scipy import weave
 from kern import Kern
 from ...util.linalg import tdot
-from ...util.misc import fast_array_equal, param_to_array
+from ...util.misc import param_to_array
 from ...core.parameterization import Param
 from ...core.parameterization.transformations import Logexp
 from ...util.caching import Cache_this
+from ...core.parameterization import variational
+from psi_comp import linear_psi_comp
 
 class Linear(Kern):
     """
@@ -104,49 +106,113 @@ class Linear(Kern):
     #---------------------------------------#
 
     def psi0(self, Z, variational_posterior):
-        return np.sum(self.variances * self._mu2S(variational_posterior), 1)
+        if isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
+            gamma = variational_posterior.binary_prob
+            mu = variational_posterior.mean
+            S = variational_posterior.variance
+
+            return np.einsum('q,nq,nq->n',self.variances,gamma,np.square(mu)+S)
+#            return (self.variances*gamma*(np.square(mu)+S)).sum(axis=1)
+        else:
+            return np.sum(self.variances * self._mu2S(variational_posterior), 1)
 
     def psi1(self, Z, variational_posterior):
-        return self.K(variational_posterior.mean, Z) #the variance, it does nothing
+        if isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
+            gamma = variational_posterior.binary_prob
+            mu = variational_posterior.mean
+            return np.einsum('nq,q,mq,nq->nm',gamma,self.variances,Z,mu)
+#            return (self.variances*gamma*mu).sum(axis=1)       
+        else:
+            return self.K(variational_posterior.mean, Z) #the variance, it does nothing
 
     @Cache_this(limit=1)
     def psi2(self, Z, variational_posterior):
-        ZA = Z * self.variances
+        if isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
-        ZAinner = self._ZAinner(variational_posterior, Z)
+            gamma = variational_posterior.binary_prob
-        return np.dot(ZAinner, ZA.T)
+            mu = variational_posterior.mean
+            S = variational_posterior.variance
+            mu2 = np.square(mu)
+            variances2 = np.square(self.variances)
+            tmp = np.einsum('nq,q,mq,nq->nm',gamma,self.variances,Z,mu)
+            return np.einsum('nq,q,mq,oq,nq->nmo',gamma,variances2,Z,Z,mu2+S)+\
+                np.einsum('nm,no->nmo',tmp,tmp) - np.einsum('nq,q,mq,oq,nq->nmo',np.square(gamma),variances2,Z,Z,mu2)
+        else:
+            ZA = Z * self.variances
+            ZAinner = self._ZAinner(variational_posterior, Z)
+            return np.dot(ZAinner, ZA.T)
 
     def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
-        #psi1
+        if isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
-        self.update_gradients_full(dL_dpsi1, variational_posterior.mean, Z)
+            gamma = variational_posterior.binary_prob
-        # psi0:
+            mu = variational_posterior.mean
-        tmp = dL_dpsi0[:, None] * self._mu2S(variational_posterior)
+            S = variational_posterior.variance
-        if self.ARD: self.variances.gradient += tmp.sum(0)
+            mu2S = np.square(mu)+S
-        else: self.variances.gradient += tmp.sum()
+            
-        #psi2
+            _dpsi2_dvariance, _, _, _, _ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma)
-        if self.ARD:
+            grad = np.einsum('n,nq,nq->q',dL_dpsi0,gamma,mu2S) + np.einsum('nm,nq,mq,nq->q',dL_dpsi1,gamma,Z,mu) +\
-            tmp = dL_dpsi2[:, :, :, None] * (self._ZAinner(variational_posterior, Z)[:, :, None, :] * Z[None, None, :, :])
+                 np.einsum('nmo,nmoq->q',dL_dpsi2,_dpsi2_dvariance)
-            self.variances.gradient += 2.*tmp.sum(0).sum(0).sum(0)
+            if self.ARD:
+                self.variances.gradient = grad
+            else:
+                self.variances.gradient = grad.sum()
         else:
-            self.variances.gradient += 2.*np.sum(dL_dpsi2 * self.psi2(Z, variational_posterior))/self.variances
+            #psi1
+            self.update_gradients_full(dL_dpsi1, variational_posterior.mean, Z)
+            # psi0:
+            tmp = dL_dpsi0[:, None] * self._mu2S(variational_posterior)
+            if self.ARD: self.variances.gradient += tmp.sum(0)
+            else: self.variances.gradient += tmp.sum()
+            #psi2
+            if self.ARD:
+                tmp = dL_dpsi2[:, :, :, None] * (self._ZAinner(variational_posterior, Z)[:, :, None, :] * Z[None, None, :, :])
+                self.variances.gradient += 2.*tmp.sum(0).sum(0).sum(0)
+            else:
+                self.variances.gradient += 2.*np.sum(dL_dpsi2 * self.psi2(Z, variational_posterior))/self.variances
 
     def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
-        #psi1
+        if isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
-        grad = self.gradients_X(dL_dpsi1.T, Z, variational_posterior.mean)
+            gamma = variational_posterior.binary_prob
-        #psi2
+            mu = variational_posterior.mean
-        self._weave_dpsi2_dZ(dL_dpsi2, Z, variational_posterior, grad)
+            S = variational_posterior.variance
-        return grad
+            _, _, _, _, _dpsi2_dZ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma)
+            
+            grad =  np.einsum('nm,nq,q,nq->mq',dL_dpsi1,gamma, self.variances,mu) +\
+                 np.einsum('nmo,noq->mq',dL_dpsi2,_dpsi2_dZ)
+            
+            return grad
+        else:
+            #psi1
+            grad = self.gradients_X(dL_dpsi1.T, Z, variational_posterior.mean)
+            #psi2
+            self._weave_dpsi2_dZ(dL_dpsi2, Z, variational_posterior, grad)
+            return grad
 
     def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
-        grad_mu, grad_S = np.zeros(variational_posterior.mean.shape), np.zeros(variational_posterior.mean.shape)
+        if isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
-        # psi0
+            gamma = variational_posterior.binary_prob
-        grad_mu += dL_dpsi0[:, None] * (2.0 * variational_posterior.mean * self.variances)
+            mu = variational_posterior.mean
-        grad_S += dL_dpsi0[:, None] * self.variances
+            S = variational_posterior.variance
-        # psi1
+            mu2S = np.square(mu)+S            
-        grad_mu += (dL_dpsi1[:, :, None] * (Z * self.variances)).sum(1)
+            _, _dpsi2_dgamma, _dpsi2_dmu, _dpsi2_dS, _ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma)
-        # psi2
+            
-        self._weave_dpsi2_dmuS(dL_dpsi2, Z, variational_posterior, grad_mu, grad_S)
+            grad_gamma = np.einsum('n,q,nq->nq',dL_dpsi0,self.variances,mu2S) + np.einsum('nm,q,mq,nq->nq',dL_dpsi1,self.variances,Z,mu) +\
-
+                 np.einsum('nmo,nmoq->nq',dL_dpsi2,_dpsi2_dgamma)
-        return grad_mu, grad_S
+            grad_mu = np.einsum('n,nq,q,nq->nq',dL_dpsi0,gamma,2.*self.variances,mu) + np.einsum('nm,nq,q,mq->nq',dL_dpsi1,gamma,self.variances,Z) +\
+                 np.einsum('nmo,nmoq->nq',dL_dpsi2,_dpsi2_dmu)
+            grad_S = np.einsum('n,nq,q->nq',dL_dpsi0,gamma,self.variances) + np.einsum('nmo,nmoq->nq',dL_dpsi2,_dpsi2_dS)
+            
+            return grad_mu, grad_S, grad_gamma
+        else:
+            grad_mu, grad_S = np.zeros(variational_posterior.mean.shape), np.zeros(variational_posterior.mean.shape)
+            # psi0
+            grad_mu += dL_dpsi0[:, None] * (2.0 * variational_posterior.mean * self.variances)
+            grad_S += dL_dpsi0[:, None] * self.variances
+            # psi1
+            grad_mu += (dL_dpsi1[:, :, None] * (Z * self.variances)).sum(1)
+            # psi2
+            self._weave_dpsi2_dmuS(dL_dpsi2, Z, variational_posterior, grad_mu, grad_S)
+    
+            return grad_mu, grad_S
 
     #--------------------------------------------------#
     #            Helpers for psi statistics            #

GPy/kern/_src/psi_comp/linear_psi_comp.py

@@ -0,0 +1,51 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+"""
+The package for the Psi statistics computation of the linear kernel for SSGPLVM
+"""
+
+import numpy as np
+from GPy.util.caching import Cache_this
+
+#@Cache_this(limit=1)
+def _psi2computations(variance, Z, mu, S, gamma):
+    """
+    Z - MxQ
+    mu - NxQ
+    S - NxQ
+    gamma - NxQ
+    """
+    # here are the "statistics" for psi1 and psi2
+    # Produced intermediate results:
+    # _psi2                NxMxM
+    # _psi2_dvariance      NxMxMxQ
+    # _psi2_dZ             NxMxQ
+    # _psi2_dgamma         NxMxMxQ
+    # _psi2_dmu            NxMxMxQ
+    # _psi2_dS             NxMxMxQ
+    
+    mu2 = np.square(mu)
+    gamma2 = np.square(gamma)
+    variance2 = np.square(variance)
+    mu2S = mu2+S # NxQ
+    common_sum = np.einsum('nq,q,mq,nq->nm',gamma,variance,Z,mu) # NxM
+            
+    _dpsi2_dvariance = np.einsum('nq,q,mq,oq->nmoq',2.*(gamma*mu2S-gamma2*mu2),variance,Z,Z)+\
+        np.einsum('nq,mq,nq,no->nmoq',gamma,Z,mu,common_sum)+\
+        np.einsum('nq,oq,nq,nm->nmoq',gamma,Z,mu,common_sum)
+        
+    _dpsi2_dgamma = np.einsum('q,mq,oq,nq->nmoq',variance2,Z,Z,(mu2S-2.*gamma*mu2))+\
+        np.einsum('q,mq,nq,no->nmoq',variance,Z,mu,common_sum)+\
+        np.einsum('q,oq,nq,nm->nmoq',variance,Z,mu,common_sum)
+    
+    _dpsi2_dmu = np.einsum('q,mq,oq,nq,nq->nmoq',variance2,Z,Z,mu,2.*(gamma-gamma2))+\
+        np.einsum('nq,q,mq,no->nmoq',gamma,variance,Z,common_sum)+\
+        np.einsum('nq,q,oq,nm->nmoq',gamma,variance,Z,common_sum)
+        
+    _dpsi2_dS = np.einsum('nq,q,mq,oq->nmoq',gamma,variance2,Z,Z)
+    
+    _dpsi2_dZ = 2.*(np.einsum('nq,q,mq,nq->nmq',gamma,variance2,Z,mu2S)+np.einsum('nq,q,nq,nm->nmq',gamma,variance,mu,common_sum)
+                    -np.einsum('nq,q,mq,nq->nmq',gamma2,variance2,Z,mu2))
+
+    return _dpsi2_dvariance, _dpsi2_dgamma, _dpsi2_dmu, _dpsi2_dS, _dpsi2_dZ

GPy/kern/_src/rbf.py

@@ -8,7 +8,7 @@ from ...util.misc import param_to_array
 from stationary import Stationary
 from GPy.util.caching import Cache_this
 from ...core.parameterization import variational
-from rbf_psi_comp import ssrbf_psi_comp
+from psi_comp import ssrbf_psi_comp
 
 class RBF(Stationary):
     """

GPy/kern/_src/ssrbf.py

@@ -7,7 +7,7 @@ import numpy as np
 from ...util.linalg import tdot
 from ...util.config import *
 from stationary import Stationary
-from rbf_psi_comp import ssrbf_psi_comp
+from psi_comp import ssrbf_psi_comp
 
 class SSRBF(Stationary):
     """

GPy/kern/_src/sympykern.py

@@ -1,11 +1,10 @@
 # Check Matthew Rocklin's blog post.
-try: 
+try:
     import sympy as sp
     sympy_available=True
     from sympy.utilities.lambdify import lambdify
 except ImportError:
     sympy_available=False
-    exit()
 
 import numpy as np
 from kern import Kern
@@ -36,7 +35,7 @@ class Sympykern(Kern):
         super(Sympykern, self).__init__(input_dim, name)
 
         self._sp_k = k
-        
+
         # pull the variable names out of the symbolic covariance function.
         sp_vars = [e for e in k.atoms() if e.is_Symbol]
         self._sp_x= sorted([e for e in sp_vars if e.name[0:2]=='x_'],key=lambda x:int(x.name[2:]))
@@ -51,7 +50,7 @@ class Sympykern(Kern):
         self._sp_kdiag = k
         for x, z in zip(self._sp_x, self._sp_z):
             self._sp_kdiag = self._sp_kdiag.subs(z, x)
-            
+
         # If it is a multi-output covariance, add an input for indexing the outputs.
         self._real_input_dim = x_dim
         # Check input dim is number of xs + 1 if output_dim is >1
@@ -73,7 +72,7 @@ class Sympykern(Kern):
 
             # Extract names of shared parameters (those without a subscript)
             self._sp_theta = [theta for theta in thetas if theta not in self._sp_theta_i and theta not in self._sp_theta_j]
-            
+
             self.num_split_params = len(self._sp_theta_i)
             self._split_theta_names = ["%s"%theta.name[:-2] for theta in self._sp_theta_i]
             # Add split parameters to the model.
@@ -82,11 +81,11 @@ class Sympykern(Kern):
                 setattr(self, theta, Param(theta, np.ones(self.output_dim), None))
                 self.add_parameter(getattr(self, theta))
 
-            
+
             self.num_shared_params = len(self._sp_theta)
             for theta_i, theta_j in zip(self._sp_theta_i, self._sp_theta_j):
                 self._sp_kdiag = self._sp_kdiag.subs(theta_j, theta_i)
-            
+
         else:
             self.num_split_params = 0
             self._split_theta_names = []
@@ -107,10 +106,10 @@ class Sympykern(Kern):
         derivative_arguments = self._sp_x + self._sp_theta
         if self.output_dim > 1:
             derivative_arguments += self._sp_theta_i
-        
+
         self.derivatives = {theta.name : sp.diff(self._sp_k,theta).simplify() for theta in derivative_arguments}
         self.diag_derivatives = {theta.name : sp.diff(self._sp_kdiag,theta).simplify() for theta in derivative_arguments}
-        
+
         # This gives the parameters for the arg list.
         self.arg_list = self._sp_x + self._sp_z + self._sp_theta
         self.diag_arg_list = self._sp_x + self._sp_theta
@@ -137,7 +136,7 @@ class Sympykern(Kern):
         for key in self.derivatives.keys():
             setattr(self, '_Kdiag_diff_' + key, lambdify(self.diag_arg_list, self.diag_derivatives[key], 'numpy'))
 
-    def K(self,X,X2=None):        
+    def K(self,X,X2=None):
         self._K_computations(X, X2)
         return self._K_function(**self._arguments)
 
@@ -145,11 +144,11 @@ class Sympykern(Kern):
     def Kdiag(self,X):
         self._K_computations(X)
         return self._Kdiag_function(**self._diag_arguments)
-    
+
     def _param_grad_helper(self,partial,X,Z,target):
         pass
 
-               
+
     def gradients_X(self, dL_dK, X, X2=None):
         #if self._X is None or X.base is not self._X.base or X2 is not None:
         self._K_computations(X, X2)
@@ -168,7 +167,7 @@ class Sympykern(Kern):
             gf = getattr(self, '_Kdiag_diff_' + x.name)
             dX[:, i] = gf(**self._diag_arguments)*dL_dK
         return dX
-    
+
     def update_gradients_full(self, dL_dK, X, X2=None):
         # Need to extract parameters to local variables first
         self._K_computations(X, X2)
@@ -193,7 +192,7 @@ class Sympykern(Kern):
                     gradient += np.asarray([A[np.where(self._output_ind2==i)].T.sum()
                                  for i in np.arange(self.output_dim)])
                 setattr(parameter, 'gradient', gradient)
-                
+
 
     def update_gradients_diag(self, dL_dKdiag, X):
         self._K_computations(X)
@@ -209,7 +208,7 @@ class Sympykern(Kern):
                 setattr(parameter, 'gradient',
                         np.asarray([a[np.where(self._output_ind==i)].sum()
                          for i in np.arange(self.output_dim)]))
-    
+
     def _K_computations(self, X, X2=None):
         """Set up argument lists for the derivatives."""
         # Could check if this needs doing or not, there could

GPy/likelihoods/likelihood.py

@@ -358,7 +358,7 @@ class Likelihood(Parameterized):
 
         return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta
 
-    def predictive_values(self, mu, var, full_cov=False, sampling=False, num_samples=10000):
+    def predictive_values(self, mu, var, full_cov=False, sampling=True, num_samples=10000):
         """
         Compute  mean, variance and conficence interval (percentiles 5 and 95) of the  prediction.
 

GPy/plotting/matplot_dep/variational_plots.py

@@ -44,3 +44,48 @@ def plot(parameterized, fignum=None, ax=None, colors=None):
     pb.draw()
     fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
     return fig
+
+def plot_SpikeSlab(parameterized, fignum=None, ax=None, colors=None):
+    """
+    Plot latent space X in 1D:
+
+        - if fig is given, create input_dim subplots in fig and plot in these
+        - if ax is given plot input_dim 1D latent space plots of X into each `axis`
+        - if neither fig nor ax is given create a figure with fignum and plot in there
+
+    colors:
+        colors of different latent space dimensions input_dim
+
+    """
+    if ax is None:
+        fig = pb.figure(num=fignum, figsize=(8, min(12, (2 * parameterized.mean.shape[1]))))
+    if colors is None:
+        colors = pb.gca()._get_lines.color_cycle
+        pb.clf()
+    else:
+        colors = iter(colors)
+    plots = []
+    means, variances, gamma = param_to_array(parameterized.mean, parameterized.variance, parameterized.binary_prob)
+    x = np.arange(means.shape[0])
+    for i in range(means.shape[1]):
+        # mean and variance plot
+        a = fig.add_subplot(means.shape[1]*2, 1, 2*i + 1)
+        a.plot(means, c='k', alpha=.3)
+        plots.extend(a.plot(x, means.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
+        a.fill_between(x,
+                        means.T[i] - 2 * np.sqrt(variances.T[i]),
+                        means.T[i] + 2 * np.sqrt(variances.T[i]),
+                        facecolor=plots[-1].get_color(),
+                        alpha=.3)
+        a.legend(borderaxespad=0.)
+        a.set_xlim(x.min(), x.max())
+        if i < means.shape[1] - 1:
+            a.set_xticklabels('')
+        # binary prob plot
+        a = fig.add_subplot(means.shape[1]*2, 1, 2*i + 2)
+        a.bar(x,gamma[:,i],bottom=0.,linewidth=0,align='center')
+        a.set_xlim(x.min(), x.max())
+        a.set_ylim([0.,1.])
+    pb.draw()
+    fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
+    return fig