From a6725b55e11e68ba5d97ac37bd2a4e34864031f6 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 16 Dec 2013 11:37:42 +0000 Subject: [PATCH] pca adjustements to lvm models --- GPy/models_modules/bayesian_gplvm.py | 92 +++++++++++++++++++++++----- GPy/models_modules/gplvm.py | 21 +++---- GPy/models_modules/mrd.py | 14 ++--- 3 files changed, 94 insertions(+), 33 deletions(-) diff --git a/GPy/models_modules/bayesian_gplvm.py b/GPy/models_modules/bayesian_gplvm.py index 90e54111..57e50955 100644 --- a/GPy/models_modules/bayesian_gplvm.py +++ b/GPy/models_modules/bayesian_gplvm.py @@ -2,17 +2,18 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import numpy as np +import itertools +from matplotlib import pyplot + from ..core.sparse_gp import SparseGP from ..likelihoods import Gaussian from .. import kern -import itertools -from matplotlib.colors import colorConverter -from GPy.inference.optimization import SCG -from GPy.util import plot_latent, linalg -from .gplvm import GPLVM -from GPy.util.plot_latent import most_significant_input_dimensions -from matplotlib import pyplot -from GPy.core.model import Model +from ..inference.optimization import SCG +from ..util import plot_latent, linalg +from .gplvm import GPLVM, initialise_latent +from ..util.plot_latent import most_significant_input_dimensions +from ..core.model import Model +from ..util.subarray_and_sorting import common_subarrays class BayesianGPLVM(SparseGP, GPLVM): """ @@ -34,7 +35,7 @@ class BayesianGPLVM(SparseGP, GPLVM): likelihood = likelihood_or_Y if X == None: - X = self.initialise_latent(init, input_dim, likelihood.Y) + X = initialise_latent(init, input_dim, likelihood.Y) self.init = init if X_variance is None: @@ -308,14 +309,36 @@ class BayesianGPLVMWithMissingData(Model): :type init: 'PCA' | 'random' """ def __init__(self, likelihood_or_Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10, - Z=None, kernel=None, missing=np.nan, **kwargs): + Z=None, kernel=None, **kwargs): + #======================================================================= + # Filter Y, such that same missing data is at same positions. + # If full rows are missing, delete them entirely! if type(likelihood_or_Y) is np.ndarray: - likelihood = Gaussian(likelihood_or_Y) + Y = likelihood_or_Y + likelihood = Gaussian + params = 1. + normalize=None else: - likelihood = likelihood_or_Y + Y = likelihood_or_Y.Y + likelihood = likelihood_or_Y.__class__ + params = likelihood_or_Y._get_params() + if isinstance(likelihood_or_Y, Gaussian): + normalize = True + scale = likelihood_or_Y._scale + offset = likelihood_or_Y._offset + # Get common subrows + filter_ = np.isnan(Y) + self.subarray_indices = common_subarrays(filter_,axis=1) + likelihoods = [likelihood(Y[~np.array(v,dtype=bool),:][:,ind]) for v,ind in self.subarray_indices.iteritems()] + for l in likelihoods: + l._set_params(params) + if normalize: # get normalization in common + l._scale = scale + l._offset = offset + #======================================================================= if X == None: - X = self.initialise_latent(init, input_dim, likelihood.Y) + X = initialise_latent(init, input_dim, Y[:,np.any(np.isnan(Y),1)]) self.init = init if X_variance is None: @@ -328,13 +351,52 @@ class BayesianGPLVMWithMissingData(Model): if kernel is None: kernel = kern.rbf(input_dim) # + kern.white(input_dim) - SparseGP.__init__(self, X, likelihood, kernel, Z=Z, X_variance=X_variance, **kwargs) + self.submodels = [BayesianGPLVM(l, input_dim, X, X_variance, init, num_inducing, Z, kernel) for l in likelihoods] + self.gref = self.submodels[0] + #:type self.gref: BayesianGPLVM self.ensure_default_constraints() + def log_likelihood(self): + ll = -self.gref.KL_divergence() + for g in self.submodels: + ll += SparseGP.log_likelihood(g) + return ll + + def _log_likelihood_gradients(self): + dLdmu, dLdS = reduce(lambda a, b: [a[0] + b[0], a[1] + b[1]], (g.dL_dmuS() for g in self.bgplvms)) + dKLmu, dKLdS = self.gref.dKL_dmuS() + dLdmu -= dKLmu + dLdS -= dKLdS + dLdmuS = np.hstack((dLdmu.flatten(), dLdS.flatten())).flatten() + dldzt1 = reduce(lambda a, b: a + b, (SparseGP._log_likelihood_gradients(g)[:self.gref.num_inducing*self.gref.input_dim] for g in self.submodels)) + + return np.hstack((dLdmuS, + dldzt1, + np.hstack([np.hstack([g.dL_dtheta(), + g.likelihood._gradients(\ + partial=g.partial_for_likelihood)]) \ + for g in self.submodels]))) + + def getstate(self): + return Model.getstate(self)+[self.submodels,self.subarray_indices] + + def setstate(self, state): + self.subarray_indices = state.pop() + self.submodels = state.pop() + self.gref = self.submodels[0] + Model.setstate(self, state) + self._set_params(self._get_params()) + def _get_param_names(self): X_names = sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], []) S_names = sum([['X_variance_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], []) - return (X_names + S_names + SparseGP._get_param_names(self)) + return (X_names + S_names + SparseGP._get_param_names(self.gref)) + + def _get_params(self): + return self.gref._get_params() + def _set_params(self, x): + [g._set_params(x) for g in self.submodels] + pass diff --git a/GPy/models_modules/gplvm.py b/GPy/models_modules/gplvm.py index f27f861c..541b3176 100644 --- a/GPy/models_modules/gplvm.py +++ b/GPy/models_modules/gplvm.py @@ -10,6 +10,13 @@ from ..core import GP from ..likelihoods import Gaussian from .. import util +def initialise_latent(init, input_dim, Y): + Xr = np.random.randn(Y.shape[0], input_dim) + if init == 'pca': + from ..util.linalg import pca + PC = pca(Y, input_dim)[0] + Xr[:PC.shape[0], :PC.shape[1]] = PC + return Xr class GPLVM(GP): """ @@ -20,12 +27,12 @@ class GPLVM(GP): :param input_dim: latent dimensionality :type input_dim: int :param init: initialisation method for the latent space - :type init: 'PCA'|'random' + :type init: 'pca'|'random' """ - def __init__(self, Y, input_dim, init='PCA', X=None, kernel=None, normalize_Y=False): + def __init__(self, Y, input_dim, init='pca', X=None, kernel=None, normalize_Y=False): if X is None: - X = self.initialise_latent(init, input_dim, Y) + X = initialise_latent(init, input_dim, Y) if kernel is None: kernel = kern.rbf(input_dim, ARD=input_dim > 1) + kern.bias(input_dim, np.exp(-2)) likelihood = Gaussian(Y, normalize=normalize_Y, variance=np.exp(-2.)) @@ -33,14 +40,6 @@ class GPLVM(GP): self.set_prior('.*X', priors.Gaussian(0, 1)) self.ensure_default_constraints() - def initialise_latent(self, init, input_dim, Y): - Xr = np.random.randn(Y.shape[0], input_dim) - if init == 'PCA': - from ..util.linalg import PCA - PC = PCA(Y, input_dim)[0] - Xr[:PC.shape[0], :PC.shape[1]] = PC - return Xr - def _get_param_names(self): return sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], []) + GP._get_param_names(self) diff --git a/GPy/models_modules/mrd.py b/GPy/models_modules/mrd.py index b9c99a64..2376993d 100644 --- a/GPy/models_modules/mrd.py +++ b/GPy/models_modules/mrd.py @@ -5,7 +5,7 @@ Created on 10 Apr 2013 ''' from GPy.core import Model from GPy.core import SparseGP -from GPy.util.linalg import PCA +from GPy.util.linalg import pca import numpy import itertools import pylab @@ -26,8 +26,8 @@ class MRD(Model): :type input_dim: int :param initx: initialisation method for the latent space : - * 'concat' - PCA on concatenation of all datasets - * 'single' - Concatenation of PCA on datasets, respectively + * 'concat' - pca on concatenation of all datasets + * 'single' - Concatenation of pca on datasets, respectively * 'random' - Random draw from a normal :type initx: ['concat'|'single'|'random'] @@ -42,7 +42,7 @@ class MRD(Model): """ def __init__(self, likelihood_or_Y_list, input_dim, num_inducing=10, names=None, - kernels=None, initx='PCA', + kernels=None, initx='pca', initz='permute', _debug=False, **kw): if names is None: self.names = ["{}".format(i) for i in range(len(likelihood_or_Y_list))] @@ -237,7 +237,7 @@ class MRD(Model): partial=g.partial_for_likelihood)]) \ for g in self.bgplvms]))) - def _init_X(self, init='PCA', likelihood_list=None): + def _init_X(self, init='pca', likelihood_list=None): if likelihood_list is None: likelihood_list = self.likelihood_list Ylist = [] @@ -248,11 +248,11 @@ class MRD(Model): Ylist.append(likelihood_or_Y.Y) del likelihood_list if init in "PCA_concat": - X = PCA(numpy.hstack(Ylist), self.input_dim)[0] + X = pca(numpy.hstack(Ylist), self.input_dim)[0] elif init in "PCA_single": X = numpy.zeros((Ylist[0].shape[0], self.input_dim)) for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.input_dim), len(Ylist)), Ylist): - X[:, qs] = PCA(Y, len(qs))[0] + X[:, qs] = pca(Y, len(qs))[0] else: # init == 'random': X = numpy.random.randn(Ylist[0].shape[0], self.input_dim) self.X = X