diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py index 87b291a7..fa4ec62f 100644 --- a/GPy/core/parameterization/variational.py +++ b/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) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 0e440ab7..dfd922f0 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/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]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 - ZAinner = self._ZAinner(variational_posterior, Z) - return np.dot(ZAinner, ZA.T) + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + gamma = variational_posterior.binary_prob + 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 - 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) + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + gamma = variational_posterior.binary_prob + mu = variational_posterior.mean + S = variational_posterior.variance + mu2S = np.square(mu)+S + + _dpsi2_dvariance, _, _, _, _ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma) + grad = np.einsum('n,nq,nq->q',dL_dpsi0,gamma,mu2S) + np.einsum('nm,nq,mq,nq->q',dL_dpsi1,gamma,Z,mu) +\ + np.einsum('nmo,nmoq->q',dL_dpsi2,_dpsi2_dvariance) + 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 - grad = self.gradients_X(dL_dpsi1.T, Z, variational_posterior.mean) - #psi2 - self._weave_dpsi2_dZ(dL_dpsi2, Z, variational_posterior, grad) - return grad + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + gamma = variational_posterior.binary_prob + mu = variational_posterior.mean + S = variational_posterior.variance + _, _, _, _, _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) - # 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 + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + gamma = variational_posterior.binary_prob + mu = variational_posterior.mean + S = variational_posterior.variance + mu2S = np.square(mu)+S + _, _dpsi2_dgamma, _dpsi2_dmu, _dpsi2_dS, _ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma) + + 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) + 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 # diff --git a/GPy/kern/_src/rbf_psi_comp/__init__.py b/GPy/kern/_src/psi_comp/__init__.py similarity index 100% rename from GPy/kern/_src/rbf_psi_comp/__init__.py rename to GPy/kern/_src/psi_comp/__init__.py diff --git a/GPy/kern/_src/psi_comp/linear_psi_comp.py b/GPy/kern/_src/psi_comp/linear_psi_comp.py new file mode 100644 index 00000000..22147366 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py b/GPy/kern/_src/psi_comp/ssrbf_psi_comp.py similarity index 100% rename from GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py rename to GPy/kern/_src/psi_comp/ssrbf_psi_comp.py diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 498ab0ac..cd6c41e9 100644 --- a/GPy/kern/_src/rbf.py +++ b/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): """ diff --git a/GPy/kern/_src/ssrbf.py b/GPy/kern/_src/ssrbf.py index 391ef1c7..c566c414 100644 --- a/GPy/kern/_src/ssrbf.py +++ b/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): """ diff --git a/GPy/kern/_src/sympykern.py b/GPy/kern/_src/sympykern.py index 920f47f3..9878ec68 100644 --- a/GPy/kern/_src/sympykern.py +++ b/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 diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index 701a5a2f..aff55533 100644 --- a/GPy/likelihoods/likelihood.py +++ b/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. diff --git a/GPy/plotting/matplot_dep/variational_plots.py b/GPy/plotting/matplot_dep/variational_plots.py index 72b857a6..cf00d8a2 100644 --- a/GPy/plotting/matplot_dep/variational_plots.py +++ b/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