From e99993b8ce0cd0b1ad3084a01dc593926a38222b Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 11 Apr 2013 15:11:41 +0100 Subject: [PATCH] GPLVM readded --- GPy/models/GPLVM.py | 125 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 GPy/models/GPLVM.py diff --git a/GPy/models/GPLVM.py b/GPy/models/GPLVM.py new file mode 100644 index 00000000..32594594 --- /dev/null +++ b/GPy/models/GPLVM.py @@ -0,0 +1,125 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + + +import numpy as np +import pylab as pb +import sys, pdb +from .. import kern +from ..core import model +from ..util.linalg import pdinv, PCA +from GP import GP +from ..likelihoods import Gaussian +from .. import util + +class GPLVM(GP): + """ + Gaussian Process Latent Variable Model + + :param Y: observed data + :type Y: np.ndarray + :param Q: latent dimensionality + :type Q: int + :param init: initialisation method for the latent space + :type init: 'PCA'|'random' + + """ + def __init__(self, Y, Q, init='PCA', X = None, kernel=None, **kwargs): + if X is None: + X = self.initialise_latent(init, Q, Y) + if kernel is None: + kernel = kern.rbf(Q) + kern.bias(Q) + likelihood = Gaussian(Y) + GP.__init__(self, X, likelihood, kernel, **kwargs) + + def initialise_latent(self, init, Q, Y): + if init == 'PCA': + return PCA(Y, Q)[0] + else: + return np.random.randn(Y.shape[0], Q) + + def _get_param_names(self): + return sum([['X_%i_%i'%(n,q) for q in range(self.Q)] for n in range(self.N)],[]) + GP._get_param_names(self) + + def _get_params(self): + return np.hstack((self.X.flatten(), GP._get_params(self))) + + def _set_params(self,x): + self.X = x[:self.X.size].reshape(self.N,self.Q).copy() + GP._set_params(self, x[self.X.size:]) + + def _log_likelihood_gradients(self): + dL_dX = 2.*self.kern.dK_dX(self.dL_dK,self.X) + + return np.hstack((dL_dX.flatten(),GP._log_likelihood_gradients(self))) + + def plot(self): + assert self.likelihood.Y.shape[1]==2 + pb.scatter(self.likelihood.Y[:,0],self.likelihood.Y[:,1],40,self.X[:,0].copy(),linewidth=0,cmap=pb.cm.jet) + Xnew = np.linspace(self.X.min(),self.X.max(),200)[:,None] + mu, var, upper, lower = self.predict(Xnew) + pb.plot(mu[:,0], mu[:,1],'k',linewidth=1.5) + + def plot_latent(self,labels=None, which_indices=None, resolution=50): + """ + :param labels: a np.array of size self.N containing labels for the points (can be number, strings, etc) + :param resolution: the resolution of the grid on which to evaluate the predictive variance + """ + + util.plot.Tango.reset() + + if labels is None: + labels = np.ones(self.N) + if which_indices is None: + if self.Q==1: + input_1 = 0 + input_2 = None + if self.Q==2: + input_1, input_2 = 0,1 + else: + #try to find a linear of RBF kern in the kernel + k = [p for p in self.kern.parts if p.name in ['rbf','linear']] + if (not len(k)==1) or (not k[0].ARD): + raise ValueError, "cannot Atomatically determine which dimensions to plot, please pass 'which_indices'" + k = k[0] + if k.name=='rbf': + input_1, input_2 = np.argsort(k.lengthscale)[:2] + elif k.name=='linear': + input_1, input_2 = np.argsort(k.variances)[::-1][:2] + + #first, plot the output variance as a function of the latent space + Xtest, xx,yy,xmin,xmax = util.plot.x_frame2D(self.X[:,[input_1, input_2]],resolution=resolution) + Xtest_full = np.zeros((Xtest.shape[0], self.X.shape[1])) + Xtest_full[:, :2] = Xtest + mu, var, low, up = self.predict(Xtest_full) + var = var[:, :2] + pb.imshow(var.reshape(resolution,resolution).T[::-1,:],extent=[xmin[0],xmax[0],xmin[1],xmax[1]],cmap=pb.cm.binary,interpolation='bilinear') + + + for i,ul in enumerate(np.unique(labels)): + if type(ul) is np.string_: + this_label = ul + elif type(ul) is np.int64: + this_label = 'class %i'%ul + else: + this_label = 'class %i'%i + + index = np.nonzero(labels==ul)[0] + if self.Q==1: + x = self.X[index,input_1] + y = np.zeros(index.size) + else: + x = self.X[index,input_1] + y = self.X[index,input_2] + pb.plot(x,y,marker='o',color=util.plot.Tango.nextMedium(),mew=0,label=this_label,linewidth=0) + + pb.xlabel('latent dimension %i'%input_1) + pb.ylabel('latent dimension %i'%input_2) + + if not np.all(labels==1.): + pb.legend(loc=0,numpoints=1) + + pb.xlim(xmin[0],xmax[0]) + pb.ylim(xmin[1],xmax[1]) + + return input_1, input_2