pca adjustements to lvm models

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 ..inference.optimization import SCG
-from matplotlib.colors import colorConverter
+from ..util import plot_latent, linalg
-from GPy.inference.optimization import SCG
+from .gplvm import GPLVM, initialise_latent
-from GPy.util import plot_latent, linalg
+from ..util.plot_latent import most_significant_input_dimensions
-from .gplvm import GPLVM
+from ..core.model import Model
-from GPy.util.plot_latent import most_significant_input_dimensions
+from ..util.subarray_and_sorting import common_subarrays
-from matplotlib import pyplot
-from GPy.core.model import Model
 
 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
 

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)
 

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
+        * 'concat' - pca on concatenation of all datasets
-        * 'single' - Concatenation of PCA on datasets, respectively
+        * '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