From c6f2082839d23074b8ae4a70d508b4e932199b8b Mon Sep 17 00:00:00 2001
From: Ricardo Andrade <acq11ra@sheffield.ac.uk>
Date: Wed, 6 Mar 2013 15:43:58 +0000
Subject: [PATCH 1/8] Sparse GP with EP is working now

---
 GPy/likelihoods/EP.py                   | 64 ++++++++++++-------------
 GPy/likelihoods/likelihood_functions.py |  4 +-
 GPy/models/sparse_GP.py                 | 36 ++++++++++----
 GPy/testing/unit_tests.py               | 20 +++++++-
 4 files changed, 78 insertions(+), 46 deletions(-)

diff --git a/GPy/likelihoods/EP.py b/GPy/likelihoods/EP.py
index efd887ae..cddc46ee 100644
--- a/GPy/likelihoods/EP.py
+++ b/GPy/likelihoods/EP.py
@@ -17,7 +17,7 @@ class EP(likelihood):
         self.epsilon = epsilon
         self.eta, self.delta = power_ep
         self.data = data
-        self.N = self.data.size
+        self.N, self.D = self.data.shape
         self.is_heteroscedastic = True
         self.Nparams = 0
 
@@ -29,7 +29,7 @@ class EP(likelihood):
         #initial values for the GP variables
         self.Y = np.zeros((self.N,1))
         self.covariance_matrix = np.eye(self.N)
-        self.precision = np.ones(self.N)
+        self.precision = np.ones(self.N)[:,None]
         self.Z = 0
         self.YYT = None
 
@@ -54,18 +54,14 @@ class EP(likelihood):
 
         self.Y =  mu_tilde[:,None]
         self.YYT = np.dot(self.Y,self.Y.T)
-        self.precision = self.tau_tilde
-        self.covariance_matrix = np.diag(1./self.precision)
+        self.covariance_matrix = np.diag(1./self.tau_tilde)
+        self.precision = self.tau_tilde[:,None]
 
     def fit_full(self,K):
         """
         The expectation-propagation algorithm.
         For nomenclature see Rasmussen & Williams 2006.
         """
-        #Prior distribution parameters: p(f|X) = N(f|0,K)
-
-        self.tau_tilde = np.zeros(self.N)
-        self.v_tilde = np.zeros(self.N)
         #Initial values - Posterior distribution parameters: q(f|X,Y) = N(f|mu,Sigma)
         mu = np.zeros(self.N)
         Sigma = K.copy()
@@ -124,13 +120,14 @@ class EP(likelihood):
 
         return self._compute_GP_variables()
 
-    def fit_DTC(self, Knn_diag, Kmn, Kmm):
+    #def fit_DTC(self, Knn_diag, Kmn, Kmm):
+    def fit_DTC(self, Kmm, Kmn):
         """
         The expectation-propagation algorithm with sparse pseudo-input.
         For nomenclature see ... 2013.
         """
 
-        #TODO: this doesn;t work with uncertain inputs!
+        #TODO: this doesn't work with uncertain inputs!
 
         """
         Prior approximation parameters:
@@ -158,12 +155,12 @@ class EP(likelihood):
         sigma_ = 1./tau_
         mu_ = v_/tau_
         """
-        tau_ = np.empty(self.N,dtype=float)
-        v_ = np.empty(self.N,dtype=float)
+        self.tau_ = np.empty(self.N,dtype=float)
+        self.v_ = np.empty(self.N,dtype=float)
 
         #Initial values - Marginal moments
         z = np.empty(self.N,dtype=float)
-        Z_hat = np.empty(self.N,dtype=float)
+        self.Z_hat = np.empty(self.N,dtype=float)
         phi = np.empty(self.N,dtype=float)
         mu_hat = np.empty(self.N,dtype=float)
         sigma2_hat = np.empty(self.N,dtype=float)
@@ -172,21 +169,21 @@ class EP(likelihood):
         epsilon_np1 = 1
         epsilon_np2 = 1
        	self.iterations = 0
-        np1 = [tau_tilde.copy()]
-        np2 = [v_tilde.copy()]
+        np1 = [self.tau_tilde.copy()]
+        np2 = [self.v_tilde.copy()]
         while epsilon_np1 > self.epsilon or epsilon_np2 > self.epsilon:
             update_order = np.random.permutation(self.N)
             for i in update_order:
                 #Cavity distribution parameters
-                tau_[i] = 1./Sigma_diag[i] - self.eta*tau_tilde[i]
-                v_[i] = mu[i]/Sigma_diag[i] - self.eta*v_tilde[i]
+                self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i]
+                self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i]
                 #Marginal moments
-                Z_hat[i], mu_hat[i], sigma2_hat[i] = self.likelihood_function.moments_match(self.data[i],tau_[i],v_[i])
+                self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.likelihood_function.moments_match(self.data[i],self.tau_[i],self.v_[i])
                 #Site parameters update
-                Delta_tau = delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
+                Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
                 Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
-                tau_tilde[i] = tau_tilde[i] + Delta_tau
-                v_tilde[i] = v_tilde[i] + Delta_v
+                self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau
+                self.v_tilde[i] = self.v_tilde[i] + Delta_v
                 #Posterior distribution parameters update
                 LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau
                 L = jitchol(LLT)
@@ -196,25 +193,26 @@ class EP(likelihood):
                 mu = mu + (Delta_v-Delta_tau*mu[i])*si
                 self.iterations += 1
             #Sigma recomputation with Cholesky decompositon
-            LLT0 = LLT0 + np.dot(Kmn*tau_tilde[None,:],Kmn.T)
+            LLT0 = LLT0 + np.dot(Kmn*self.tau_tilde[None,:],Kmn.T)
             L = jitchol(LLT)
             V,info = linalg.flapack.dtrtrs(L,Kmn,lower=1)
             V2,info = linalg.flapack.dtrtrs(L.T,V,lower=0)
             Sigma_diag = np.sum(V*V,-2)
-            Knmv_tilde = np.dot(Kmn,v_tilde)
+            Knmv_tilde = np.dot(Kmn,self.v_tilde)
             mu = np.dot(V2.T,Knmv_tilde)
-            epsilon_np1 = sum((tau_tilde-np1[-1])**2)/self.N
-            epsilon_np2 = sum((v_tilde-np2[-1])**2)/self.N
-            np1.append(tau_tilde.copy())
-            np2.append(v_tilde.copy())
+            epsilon_np1 = sum((self.tau_tilde-np1[-1])**2)/self.N
+            epsilon_np2 = sum((self.v_tilde-np2[-1])**2)/self.N
+            np1.append(self.tau_tilde.copy())
+            np2.append(self.v_tilde.copy())
 
         self._compute_GP_variables()
 
-    def fit_FITC(self, Knn_diag, Kmn):
+    def fit_FITC(self, Kmm, Kmn, Knn_diag):
         """
         The expectation-propagation algorithm with sparse pseudo-input.
         For nomenclature see Naish-Guzman and Holden, 2008.
         """
+        M = Kmm.shape[0]
 
         """
         Prior approximation parameters:
@@ -235,7 +233,7 @@ class EP(likelihood):
         mu = w + P*gamma
         """
         self.w = np.zeros(self.N)
-        self.gamma = np.zeros(self.M)
+        self.gamma = np.zeros(M)
         mu = np.zeros(self.N)
         P = P0.copy()
         R = R0.copy()
@@ -271,7 +269,7 @@ class EP(likelihood):
                 self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i]
                 self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i]
                 #Marginal moments
-                self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.likelihood_function.moments_match(data[i],self.tau_[i],self.v_[i])
+                self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.likelihood_function.moments_match(self.data[i],self.tau_[i],self.v_[i])
                 #Site parameters update
                 Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
                 Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
@@ -281,10 +279,10 @@ class EP(likelihood):
                 dtd1 = Delta_tau*Diag[i] + 1.
                 dii = Diag[i]
                 Diag[i] = dii - (Delta_tau * dii**2.)/dtd1
-                pi_ = P[i,:].reshape(1,self.M)
+                pi_ = P[i,:].reshape(1,M)
                 P[i,:] = pi_ - (Delta_tau*dii)/dtd1 * pi_
                 Rp_i = np.dot(R,pi_.T)
-                RTR = np.dot(R.T,np.dot(np.eye(self.M) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R))
+                RTR = np.dot(R.T,np.dot(np.eye(M) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R))
                 R = jitchol(RTR).T
                 self.w[i] = self.w[i] + (Delta_v - Delta_tau*self.w[i])*dii/dtd1
                 self.gamma = self.gamma + (Delta_v - Delta_tau*mu[i])*np.dot(RTR,P[i,:].T)
@@ -296,7 +294,7 @@ class EP(likelihood):
             Diag = Diag0/(1.+ Diag0 * self.tau_tilde)
             P = (Diag / Diag0)[:,None] * P0
             RPT0 = np.dot(R0,P0.T)
-            L = jitchol(np.eye(self.M) + np.dot(RPT0,(1./Diag0 - Diag/(Diag0**2))[:,None]*RPT0.T))
+            L = jitchol(np.eye(M) + np.dot(RPT0,(1./Diag0 - Diag/(Diag0**2))[:,None]*RPT0.T))
             R,info = linalg.flapack.dtrtrs(L,R0,lower=1)
             RPT = np.dot(R,P.T)
             Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1)
diff --git a/GPy/likelihoods/likelihood_functions.py b/GPy/likelihoods/likelihood_functions.py
index 23881899..3e2a0361 100644
--- a/GPy/likelihoods/likelihood_functions.py
+++ b/GPy/likelihoods/likelihood_functions.py
@@ -37,8 +37,8 @@ class probit(likelihood_function):
         :param tau_i: precision of the cavity distribution (float)
         :param v_i: mean/variance of the cavity distribution (float)
         """
-        # TODO: some version of assert np.sum(np.abs(Y)-1) == 0, "Output values must be either -1 or 1"
-        if data_i == 0: data_i = -1 #NOTE Binary classification works better classes {-1,1}, 1D-plotting works better with classes {0,1}.
+        if data_i == 0: data_i = -1 #NOTE Binary classification algorithm works better with classes {-1,1}, 1D-plotting works better with classes {0,1}.
+        # TODO: some version of assert
         z = data_i*v_i/np.sqrt(tau_i**2 + tau_i)
         Z_hat = stats.norm.cdf(z)
         phi = stats.norm.pdf(z)
diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 6932154d..12cc1769 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -82,9 +82,9 @@ class sparse_GP(GP):
         if self.likelihood.is_heteroscedastic:
             assert self.likelihood.D == 1 #TODO: what is the likelihood is heterscedatic and there are multiple independent outputs?
             if self.has_uncertain_inputs:
-                self.psi2_beta_scaled = (self.psi2*(self.likelihood.precision.reshape(self.N,1,1)/sf2)).sum(0)
+                self.psi2_beta_scaled = (self.psi2*(self.likelihood.precision.flatten().reshape(self.N,1,1)/sf2)).sum(0)
             else:
-                tmp = self.psi1.T*(np.sqrt(self.likelihood.precision.reshape(1,self.N))/sf)
+                tmp = self.psi1*(np.sqrt(self.likelihood.precision.flatten().reshape(1,self.N))/sf)
                 self.psi2_beta_scaled = np.dot(tmp,tmp.T)
         else:
             if self.has_uncertain_inputs:
@@ -107,14 +107,18 @@ class sparse_GP(GP):
         self.E = mdot(self.C, self.psi1VVpsi1/sf2, self.C.T)
 
         # Compute dL_dpsi # FIXME: this is untested for the het. case
-        self.dL_dpsi0 = - 0.5 * self.D * self.likelihood.precision * np.ones(self.N)
+        self.dL_dpsi0 = - 0.5 * self.D * (self.likelihood.precision * np.ones([self.N,1])).flatten()
         self.dL_dpsi1 = mdot(self.V, self.psi1V.T,self.C).T
         if self.likelihood.is_heteroscedastic:
-            self.dL_dpsi2 = 0.5 * self.likelihood.precision[:,None,None] * self.D * self.Kmmi[None,:,:] # dB
-            self.dL_dpsi2 += - 0.5 * self.likelihood.precision[:,None,None]/sf2 * self.D * self.C[None,:,:] # dC
-            self.dL_dpsi2 += - 0.5 * self.likelihood.precision[:,None,None]* self.E[None,:,:] # dD
-            if not self.has_uncertain_inputs:
-                raise NotImplementedError, "TODO: recaste derivatibes in psi2 back into psi1"
+            if self.has_uncertain_inputs:
+                self.dL_dpsi2 = 0.5 * self.likelihood.precision[:,None,None] * self.D * self.Kmmi[None,:,:] # dB
+                self.dL_dpsi2 += - 0.5 * self.likelihood.precision[:,None,None]/sf2 * self.D * self.C[None,:,:] # dC
+                self.dL_dpsi2 += - 0.5 * self.likelihood.precision[:,None,None]* self.E[None,:,:] # dD
+            else:
+                self.dL_dpsi1 += mdot(self.Kmmi,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)) #dB
+                self.dL_dpsi1 += -mdot(self.C,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)/sf2) #dC
+                self.dL_dpsi1 += -mdot(self.E,self.psi1*self.likelihood.precision.flatten().reshape(1,self.N)) #dD
+                self.dL_dpsi2 = None
 
         else:
             self.dL_dpsi2 = 0.5 * self.likelihood.precision * self.D * self.Kmmi # dB
@@ -166,14 +170,28 @@ class sparse_GP(GP):
     def _get_param_names(self):
         return sum([['iip_%i_%i'%(i,j) for j in range(self.Z.shape[1])] for i in range(self.Z.shape[0])],[]) + GP._get_param_names(self)
 
+    def update_likelihood_approximation(self):
+        """
+        Approximates a non-gaussian likelihood using Expectation Propagation
+
+        For a Gaussian (or direct: TODO) likelihood, no iteration is required:
+        this function does nothing
+        """
+        if self.has_uncertain_inputs:
+            raise NotImplementedError, "EP approximation not implemented for uncertain inputs"
+        else:
+            self.likelihood.fit_DTC(self.Kmm,self.psi1)
+            self._set_params(self._get_params()) # update the GP
+
     def log_likelihood(self):
         """ Compute the (lower bound on the) log marginal likelihood """
         sf2 = self.scale_factor**2
         if self.likelihood.is_heteroscedastic:
             A = -0.5*self.N*self.D*np.log(2.*np.pi) +0.5*np.sum(np.log(self.likelihood.precision)) -0.5*np.sum(self.V*self.likelihood.Y)
+            B = -0.5*self.D*(np.sum(self.likelihood.precision.flatten()*self.psi0) - np.trace(self.A)*sf2)
         else:
             A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.likelihood.precision)) -0.5*self.likelihood.precision*self.likelihood.trYYT
-        B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2)
+            B = -0.5*self.D*(np.sum(self.likelihood.precision*self.psi0) - np.trace(self.A)*sf2)
         C = -0.5*self.D * (self.B_logdet + self.M*np.log(sf2))
         D = +0.5*np.sum(self.psi1VVpsi1 * self.C)
         return A+B+C+D
diff --git a/GPy/testing/unit_tests.py b/GPy/testing/unit_tests.py
index 61fb15bb..90037dcb 100644
--- a/GPy/testing/unit_tests.py
+++ b/GPy/testing/unit_tests.py
@@ -157,13 +157,29 @@ class GradientTests(unittest.TestCase):
     def test_GP_EP_probit(self):
         N = 20
         X = np.hstack([np.random.normal(5,2,N/2),np.random.normal(10,2,N/2)])[:,None]
-        Y = np.hstack([np.ones(N/2),np.repeat(-1,N/2)])[:,None]
+        Y = np.hstack([np.ones(N/2),np.zeros(N/2)])[:,None]
         kernel = GPy.kern.rbf(1)
         distribution = GPy.likelihoods.likelihood_functions.probit()
         likelihood = GPy.likelihoods.EP(Y, distribution)
         m = GPy.models.GP(X, likelihood, kernel)
         m.ensure_default_constraints()
-        self.assertTrue(m.EPEM)
+        m.update_likelihood_approximation()
+        self.assertTrue(m.checkgrad())
+        #self.assertTrue(m.EPEM)
+
+    def test_sparse_EP_DTC_probit(self):
+        N = 20
+        X = np.hstack([np.random.normal(5,2,N/2),np.random.normal(10,2,N/2)])[:,None]
+        Y = np.hstack([np.ones(N/2),np.zeros(N/2)])[:,None]
+        Z = np.linspace(0,15,4)[:,None]
+        kernel = GPy.kern.rbf(1)
+        distribution = GPy.likelihoods.likelihood_functions.probit()
+        likelihood = GPy.likelihoods.EP(Y, distribution)
+        m = GPy.models.sparse_GP(X, likelihood, kernel,Z)
+        m.ensure_default_constraints()
+        m.update_likelihood_approximation()
+        self.assertTrue(m.checkgrad())
+
 
     @unittest.skip("FITC will be broken for a while")
     def test_generalized_FITC(self):

From 24d705417418d152bc98e102b77d3afa7e79e694 Mon Sep 17 00:00:00 2001
From: Ricardo Andrade <acq11ra@sheffield.ac.uk>
Date: Fri, 8 Mar 2013 11:46:17 +0000
Subject: [PATCH 2/8] some small changes.

---
 GPy/models/__init__.py         |   1 +
 GPy/models/generalized_FITC.py | 162 +++++++++++++++++++++++++++++++++
 GPy/models/sparse_GP.py        |   7 +-
 3 files changed, 167 insertions(+), 3 deletions(-)
 create mode 100644 GPy/models/generalized_FITC.py

diff --git a/GPy/models/__init__.py b/GPy/models/__init__.py
index c099d0d5..61591320 100644
--- a/GPy/models/__init__.py
+++ b/GPy/models/__init__.py
@@ -11,3 +11,4 @@ from warped_GP import warpedGP
 from sparse_GPLVM import sparse_GPLVM
 from uncollapsed_sparse_GP import uncollapsed_sparse_GP
 from BGPLVM import Bayesian_GPLVM
+from generalized_FITC import generalized_FITC
diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py
new file mode 100644
index 00000000..7e0c656e
--- /dev/null
+++ b/GPy/models/generalized_FITC.py
@@ -0,0 +1,162 @@
+# 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
+from ..util.linalg import mdot, jitchol, chol_inv, pdinv
+from ..util.plot import gpplot
+from scipy import linalg
+from .. import kern
+from sparse_GP import sparse_GP
+
+"""
+import numpy as np
+import pylab as pb
+from scipy import stats, linalg
+from .. import kern
+from ..core import model
+from ..util.linalg import pdinv,mdot
+from ..util.plot import gpplot
+#from ..inference.Expectation_Propagation import FITC
+from ..likelihoods.EP import FITC
+from ..likelihoods import likelihood,probit
+"""
+
+class generalized_FITC(sparse_GP):
+    def __init__(self, X, likelihood, kernel, Z, X_uncertainty=None, Xslices=None,Zslices=None, normalize_X=False):
+    #def __init__(self, X, likelihood, kernel=None, inducing=10, epsilon_ep=1e-3, powerep=[1.,1.]):
+        """
+        Naish-Guzman, A. and Holden, S. (2008) implemantation of EP with FITC.
+
+        :param X: input observations
+        :param likelihood: Output's likelihood (likelihood class)
+        :param kernel: a GPy kernel
+        :param Z:  Either an array specifying the inducing points location or a scalar defining their number.
+        """
+
+        if type(Z) == int:
+            self.M = Z
+            self.Z = (np.random.random_sample(self.D*self.M)*(self.X.max()-self.X.min())+self.X.min()).reshape(self.M,-1)
+        elif type(Z) == np.ndarray:
+            self.Z = Z
+            self.M = self.Z.shape[0]
+
+        self._precision = likelihood.precision
+
+        sparse_GP.__init__(self, X, likelihood, kernel=kernel, Z=self.Z, X_uncertainty=None, Xslices=None,Zslices=None, normalize_X=False)
+        self.scale_factor = 100.
+
+    def update_likelihood_approximation(self):
+        """
+        Approximates a non-gaussian likelihood using Expectation Propagation
+
+        For a Gaussian (or direct: TODO) likelihood, no iteration is required:
+        this function does nothing
+        """
+        if self.has_uncertain_inputs:
+            raise NotImplementedError, "FITC approximation not implemented for uncertain inputs"
+        else:
+            self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
+            self._precision = self.likelihood.precision # Save the true precision
+            self.likelihood.precision = self.likelihood.precision/(1. + self.likelihood.precision*self.Diag0[:,None]) # Add the diagonal element of the FITC approximation
+            self._set_params(self._get_params()) # update the GP
+
+    def _set_params(self, p):
+        self.Z = p[:self.M*self.Q].reshape(self.M, self.Q)
+        self.kern._set_params(p[self.Z.size:self.Z.size+self.kern.Nparam])
+        self.likelihood._set_params(p[self.Z.size+self.kern.Nparam:])
+        self._compute_kernel_matrices()
+        self._computations()
+        self._FITC_computations()
+
+    def _FITC_computations(self):
+        """
+        FITC approximation doesn't have the correction term in the log-likelihood bound,
+        but adds a diagonal term to the covariance matrix.
+        This function:
+            - computes the diagonal term
+            - eliminates the extra terms computed in the sparse_GP approximation
+            - computes the likelihood gradients wrt the true precision.
+        """
+        # Compute FITC's diagonal term of the covariance
+        sf = self.scale_factor
+        sf2 = sf**2
+        self.Qnn = mdot(self.psi1.T,self.Kmmi,self.psi1)
+        self.Diag0 = self.psi0 - np.diag(self.Qnn)
+
+        self.Diag = self.Diag0/(1.+ self.Diag0 * self._precision.flatten())
+        self.P = (self.Diag / self.Diag0)[:,None] * self.psi1.T
+        self.RPT0 = np.dot(self.Lmi,self.psi1)
+        self.L = np.linalg.cholesky(np.eye(self.M) + np.dot(self.RPT0,(1./self.Diag0 - self.Diag/(self.Diag0**2))[:,None]*self.RPT0.T))
+        self.R,info = linalg.flapack.dtrtrs(self.L,self.Lmi,lower=1)
+        self.RPT = np.dot(self.R,self.P.T)
+        self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT)
+        self.w = self.Diag * self.likelihood.v_tilde
+        self.gamma = np.dot(self.R.T, np.dot(self.RPT,self.likelihood.v_tilde))
+        self.mu = self.w + np.dot(self.P,self.gamma)
+        self.mu_tilde = (self.likelihood.v_tilde/self.likelihood.tau_tilde)[:,None]
+
+        # Remove extra term from dL_dpsi
+        self.dL_dpsi0 = np.zeros(self.N)
+        # Remove extra term from dL_dKmm
+        self.dL_dKmm = +0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi)*sf2 # dB
+        #the partial derivative vector for the likelihood with the true precision
+        if self.likelihood.Nparams ==0:
+            #save computation here
+            self.partial_for_likelihood = None
+        elif self.likelihood.is_heteroscedastic:
+            raise NotImplementedError, "heteroscedatic derivates not implemented"
+        else:
+            beta = self.likelihood._precision # NOTE the true precison is now '_precison' not 'precision'
+            dbeta =   0.5 * self.N*self.D/beta - 0.5 * np.sum(np.square(self.likelihood.Y))
+            #dbeta += - 0.5 * self.D * (self.psi0.sum() - np.trace(self.A)/beta*sf2)
+            dbeta += - 0.5 * self.D * np.sum(self.Bi*self.A)/beta
+            dbeta += np.sum((self.C - 0.5 * mdot(self.C,self.psi2_beta_scaled,self.C) ) * self.psi1VVpsi1 )/beta
+            self.partial_for_likelihood = -dbeta*self.likelihood.precision**2
+
+
+
+
+    def _raw_predict(self, Xnew, slices, full_cov=True):
+        """
+        Make a prediction for the vsGP model
+
+        Arguments
+        ---------
+        X : Input prediction data - Nx1 numpy array (floats)
+        """
+        Kx = self.kern.K(self.Z, Xnew)
+        #K_x = self.kernel.K(self.Z,X)
+        if full_cov:
+            Kxx = self.kern.K(Xnew)
+        else:
+            Kxx = self.kern.K(Xnew)#FIXME
+            #raise NotImplementedError
+            #Kxx = self.kern.Kdiag(Xnew)
+
+        # q(u|f) = N(u| R0i*mu_u*f, R0i*C*R0i.T)
+
+        # Ci = I + (RPT0)Di(RPT0).T
+        # C = I - [RPT0] * (D+[RPT0].T*[RPT0])^-1*[RPT0].T
+        #   = I - [RPT0] * (D + self.Qnn)^-1 * [RPT0].T
+        #   = I - [RPT0] * (U*U.T)^-1 * [RPT0].T
+        #   = I - V.T * V
+        U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn)
+        V,info = linalg.flapack.dtrtrs(U,self.RPT0.T,lower=1)
+        C = np.eye(self.M) - np.dot(V.T,V)
+        mu_u = np.dot(C,self.RPT0)*(1./self.Diag0[None,:])
+        #self.C = C
+        #self.RPT0 = np.dot(self.R0,self.Knm.T) P0.T
+        #self.mu_u = mu_u
+        #self.U = U
+        # q(u|y) = N(u| R0i*mu_H,R0i*Sigma_H*R0i.T)
+        mu_H = np.dot(mu_u,self.mu)
+        self.mu_H = mu_H
+        Sigma_H = C + np.dot(mu_u,np.dot(self.Sigma,mu_u.T))
+        # q(f_star|y) = N(f_star|mu_star,sigma2_star)
+        KR0T = np.dot(Kx.T,self.Lmi.T)
+        mu_star = np.dot(KR0T,mu_H)
+        sigma2_star = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T))
+        vdiag = np.diag(sigma2_star)
+        return mu_star[:,None],vdiag[:,None]
diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 12cc1769..54eebd2f 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -72,7 +72,7 @@ class sparse_GP(GP):
             self.psi2 = None
 
     def _computations(self):
-        # TODO find routine to multiply triangular matrices
+        #TODO: find routine to multiply triangular matrices
         #TODO: slices for psi statistics (easy enough)
 
         sf = self.scale_factor
@@ -106,7 +106,7 @@ class sparse_GP(GP):
         self.C = mdot(self.Lmi.T, self.Bi, self.Lmi)
         self.E = mdot(self.C, self.psi1VVpsi1/sf2, self.C.T)
 
-        # Compute dL_dpsi # FIXME: this is untested for the het. case
+        # Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertin inputs case
         self.dL_dpsi0 = - 0.5 * self.D * (self.likelihood.precision * np.ones([self.N,1])).flatten()
         self.dL_dpsi1 = mdot(self.V, self.psi1V.T,self.C).T
         if self.likelihood.is_heteroscedastic:
@@ -180,7 +180,8 @@ class sparse_GP(GP):
         if self.has_uncertain_inputs:
             raise NotImplementedError, "EP approximation not implemented for uncertain inputs"
         else:
-            self.likelihood.fit_DTC(self.Kmm,self.psi1)
+            #self.likelihood.fit_DTC(self.Kmm,self.psi1)
+            self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
             self._set_params(self._get_params()) # update the GP
 
     def log_likelihood(self):

From 0ade786385a97f6c1fda221080430c7676609a0c Mon Sep 17 00:00:00 2001
From: Ricardo Andrade <acq11ra@sheffield.ac.uk>
Date: Mon, 11 Mar 2013 11:39:48 +0000
Subject: [PATCH 3/8] Plot function moved to GP model

---
 GPy/models/sparse_GP.py | 15 ++-------------
 1 file changed, 2 insertions(+), 13 deletions(-)

diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 54eebd2f..ff00faea 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -180,8 +180,8 @@ class sparse_GP(GP):
         if self.has_uncertain_inputs:
             raise NotImplementedError, "EP approximation not implemented for uncertain inputs"
         else:
-            #self.likelihood.fit_DTC(self.Kmm,self.psi1)
-            self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
+            self.likelihood.fit_DTC(self.Kmm,self.psi1)
+            #self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
             self._set_params(self._get_params()) # update the GP
 
     def log_likelihood(self):
@@ -240,14 +240,3 @@ class sparse_GP(GP):
             var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.C/self.scale_factor**2, Kx),0)
 
         return mu,var[:,None]
-
-    def plot(self, *args, **kwargs):
-        """
-        Plot the fitted model: just call the GP plot function and then add inducing inputs
-        """
-        GP.plot(self,*args,**kwargs)
-        if self.Q==1:
-            if self.has_uncertain_inputs:
-                pb.errorbar(self.X[:,0], pb.ylim()[0]+np.zeros(self.N), xerr=2*np.sqrt(self.X_uncertainty.flatten()))
-        if self.Q==2:
-            pb.plot(self.Z[:,0],self.Z[:,1],'wo')

From 5dbc5bdb6e283d5705c4c376aa662119bff1e808 Mon Sep 17 00:00:00 2001
From: Ricardo Andrade <acq11ra@sheffield.ac.uk>
Date: Mon, 11 Mar 2013 11:40:17 +0000
Subject: [PATCH 4/8] Plotting functions for sparse_GP added

---
 GPy/models/GP.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/GPy/models/GP.py b/GPy/models/GP.py
index 08ac1bb1..11cd174d 100644
--- a/GPy/models/GP.py
+++ b/GPy/models/GP.py
@@ -269,6 +269,8 @@ class GP(model):
             if hasattr(self,'Z'):
                 Zu = self.Z*self._Xstd + self._Xmean
                 pb.plot(Zu,Zu*0+pb.ylim()[0],'r|',mew=1.5,markersize=12)
+                if self.has_uncertain_inputs:
+                    pb.errorbar(self.X[:,0], pb.ylim()[0]+np.zeros(self.N), xerr=2*np.sqrt(self.X_uncertainty.flatten()))
 
         elif self.X.shape[1]==2: #FIXME
             resolution = resolution or 50
@@ -281,5 +283,8 @@ class GP(model):
             pb.scatter(self.X[:,0], self.X[:,1], 40, Yf, cmap=pb.cm.jet,vmin=m.min(),vmax=m.max(), linewidth=0.)
             pb.xlim(xmin[0],xmax[0])
             pb.ylim(xmin[1],xmax[1])
+            if hasattr(self,'Z'):
+                pb.plot(self.Z[:,0],self.Z[:,1],'wo')
+
         else:
             raise NotImplementedError, "Cannot define a frame with more than two input dimensions"

From 1ddc05925178857bd3e3354c1d2b99a65abdde33 Mon Sep 17 00:00:00 2001
From: Ricardo Andrade <acq11ra@sheffield.ac.uk>
Date: Mon, 11 Mar 2013 11:41:22 +0000
Subject: [PATCH 5/8] Test for EP_DTC added

---
 GPy/testing/unit_tests.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/GPy/testing/unit_tests.py b/GPy/testing/unit_tests.py
index 90037dcb..5ec1d766 100644
--- a/GPy/testing/unit_tests.py
+++ b/GPy/testing/unit_tests.py
@@ -180,7 +180,6 @@ class GradientTests(unittest.TestCase):
         m.update_likelihood_approximation()
         self.assertTrue(m.checkgrad())
 
-
     @unittest.skip("FITC will be broken for a while")
     def test_generalized_FITC(self):
         N = 20

From addb5da4e439f76b13a47a53a30ba35b6fe51bee Mon Sep 17 00:00:00 2001
From: Ricardo Andrade <acq11ra@sheffield.ac.uk>
Date: Mon, 11 Mar 2013 11:41:46 +0000
Subject: [PATCH 6/8] Irrelevant changes

---
 GPy/examples/classification.py |  22 +++----
 GPy/examples/sparse_ep_fix.py  | 113 +++++++++++++++++++++------------
 2 files changed, 82 insertions(+), 53 deletions(-)

diff --git a/GPy/examples/classification.py b/GPy/examples/classification.py
index 592299d8..031cc915 100644
--- a/GPy/examples/classification.py
+++ b/GPy/examples/classification.py
@@ -11,7 +11,7 @@ import GPy
 
 default_seed=10000
 
-def crescent_data(model_type='Full', inducing=10, seed=default_seed): #FIXME
+def crescent_data(seed=default_seed): #FIXME
     """Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
 
     :param model_type: type of model to fit ['Full', 'FITC', 'DTC'].
@@ -31,11 +31,8 @@ def crescent_data(model_type='Full', inducing=10, seed=default_seed): #FIXME
     likelihood = GPy.likelihoods.EP(data['Y'],distribution)
 
 
-    if model_type=='Full':
-        m = GPy.models.GP(data['X'],likelihood,kernel)
-    else:
-        # create sparse GP EP model
-        m = GPy.models.sparse_GP_EP(data['X'],likelihood=likelihood,inducing=inducing,ep_proxy=model_type)
+    m = GPy.models.GP(data['X'],likelihood,kernel)
+    m.ensure_default_constraints()
 
     m.update_likelihood_approximation()
     print(m)
@@ -94,16 +91,13 @@ def toy_linear_1d_classification(seed=default_seed):
 
     # Model definition
     m = GPy.models.GP(data['X'],likelihood=likelihood,kernel=kernel)
+    m.ensure_default_constraints()
 
     # Optimize
-    """
-    EPEM runs a loop that consists of two steps:
-    1) EP likelihood approximation:
-        m.update_likelihood_approximation()
-    2) Parameters optimization:
-        m.optimize()
-    """
-    m.EPEM()
+    m.update_likelihood_approximation()
+    # Parameters optimization:
+    m.optimize()
+    #m.EPEM() #FIXME
 
     # Plot
     pb.subplot(211)
diff --git a/GPy/examples/sparse_ep_fix.py b/GPy/examples/sparse_ep_fix.py
index defcb4eb..acbd506c 100644
--- a/GPy/examples/sparse_ep_fix.py
+++ b/GPy/examples/sparse_ep_fix.py
@@ -10,51 +10,86 @@ import pylab as pb
 import numpy as np
 import GPy
 np.random.seed(2)
-pb.ion()
 N = 500
 M = 5
 
-pb.close('all')
-######################################
-## 1 dimensional example
+default_seed=10000
 
-# sample inputs and outputs
-X = np.random.uniform(-3.,3.,(N,1))
-#Y = np.sin(X)+np.random.randn(N,1)*0.05
-F = np.sin(X)+np.random.randn(N,1)*0.05
-Y = np.ones([F.shape[0],1])
-Y[F<0] = -1
-likelihood = GPy.inference.likelihoods.probit(Y)
+def crescent_data(inducing=10, seed=default_seed):
+    """Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
 
-# construct kernel
-rbf =  GPy.kern.rbf(1)
-noise = GPy.kern.white(1)
-kernel = rbf + noise
+    :param model_type: type of model to fit ['Full', 'FITC', 'DTC'].
+    :param seed : seed value for data generation.
+    :type seed: int
+    :param inducing : number of inducing variables (only used for 'FITC' or 'DTC').
+    :type inducing: int
+    """
 
-# create simple GP model
-#m = GPy.models.sparse_GP(X,Y=None, kernel=kernel, M=M,likelihood= likelihood)
+    data = GPy.util.datasets.crescent_data(seed=seed)
 
-# contrain all parameters to be positive
-#m.constrain_fixed('prec',100.)
-m = GPy.models.sparse_GP(X, Y, kernel, M=M)
-m.ensure_default_constraints()
-#if not isinstance(m.likelihood,GPy.inference.likelihoods.gaussian):
-#    m.approximate_likelihood()
-print m.checkgrad()
-m.optimize('tnc', messages = 1)
-m.plot(samples=3)
-print m
+    # Kernel object
+    kernel = GPy.kern.rbf(data['X'].shape[1]) + GPy.kern.white(data['X'].shape[1])
 
-n = GPy.models.sparse_GP(X,Y=None, kernel=kernel, M=M,likelihood= likelihood)
-n.ensure_default_constraints()
-if not isinstance(n.likelihood,GPy.inference.likelihoods.gaussian):
-    n.approximate_likelihood()
-print n.checkgrad()
-pb.figure()
-n.plot()
+    # Likelihood object
+    distribution = GPy.likelihoods.likelihood_functions.probit()
+    likelihood = GPy.likelihoods.EP(data['Y'],distribution)
+
+    sample = np.random.randint(0,data['X'].shape[0],inducing)
+    Z = data['X'][sample,:]
+    #Z = (np.random.random_sample(2*inducing)*(data['X'].max()-data['X'].min())+data['X'].min()).reshape(inducing,-1)
+
+    # create sparse GP EP model
+    m = GPy.models.sparse_GP(data['X'],likelihood=likelihood,kernel=kernel,Z=Z)
+    m.ensure_default_constraints()
+
+    m.update_likelihood_approximation()
+    print(m)
+
+    # optimize
+    m.optimize()
+    print(m)
+
+    # plot
+    m.plot()
+    return m
+
+
+def toy_linear_1d_classification(seed=default_seed):
+    """
+    Simple 1D classification example
+    :param seed : seed value for data generation (default is 4).
+    :type seed: int
+    """
+
+    data = GPy.util.datasets.toy_linear_1d_classification(seed=seed)
+    Y = data['Y'][:, 0:1]
+    Y[Y == -1] = 0
+
+    # Kernel object
+    kernel = GPy.kern.rbf(1)
+
+    # Likelihood object
+    distribution = GPy.likelihoods.likelihood_functions.probit()
+    likelihood = GPy.likelihoods.EP(Y,distribution)
+
+    Z = np.random.uniform(data['X'].min(),data['X'].max(),(10,1))
+
+    # Model definition
+    m = GPy.models.sparse_GP(data['X'],likelihood=likelihood,kernel=kernel,Z=Z)
+
+    m.ensure_default_constraints()
+    # Optimize
+    m.update_likelihood_approximation()
+    # Parameters optimization:
+    m.optimize()
+    #m.EPEM() #FIXME
+
+    # Plot
+    pb.subplot(211)
+    m.plot_f()
+    pb.subplot(212)
+    m.plot()
+    print(m)
+
+    return m
 
-"""
-m = GPy.models.sparse_GP_regression(X, Y, kernel, M=M)
-m.ensure_default_constraints()
-print m.checkgrad()
-"""

From 9126c1086601bf128a8bf4842d45ce9829b2ce50 Mon Sep 17 00:00:00 2001
From: Ricardo Andrade <acq11ra@sheffield.ac.uk>
Date: Mon, 11 Mar 2013 11:42:30 +0000
Subject: [PATCH 7/8] Removed generalized_FITC.py

---
 GPy/models/__init__.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/GPy/models/__init__.py b/GPy/models/__init__.py
index 61591320..c099d0d5 100644
--- a/GPy/models/__init__.py
+++ b/GPy/models/__init__.py
@@ -11,4 +11,3 @@ from warped_GP import warpedGP
 from sparse_GPLVM import sparse_GPLVM
 from uncollapsed_sparse_GP import uncollapsed_sparse_GP
 from BGPLVM import Bayesian_GPLVM
-from generalized_FITC import generalized_FITC

From 3f88381dc9ed2cdce3e46bb31a81c1892edc552c Mon Sep 17 00:00:00 2001
From: Ricardo Andrade <acq11ra@sheffield.ac.uk>
Date: Mon, 11 Mar 2013 11:43:04 +0000
Subject: [PATCH 8/8] generalized_FITC removed

---
 GPy/models/generalized_FITC.py | 162 ---------------------------------
 1 file changed, 162 deletions(-)
 delete mode 100644 GPy/models/generalized_FITC.py

diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py
deleted file mode 100644
index 7e0c656e..00000000
--- a/GPy/models/generalized_FITC.py
+++ /dev/null
@@ -1,162 +0,0 @@
-# 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
-from ..util.linalg import mdot, jitchol, chol_inv, pdinv
-from ..util.plot import gpplot
-from scipy import linalg
-from .. import kern
-from sparse_GP import sparse_GP
-
-"""
-import numpy as np
-import pylab as pb
-from scipy import stats, linalg
-from .. import kern
-from ..core import model
-from ..util.linalg import pdinv,mdot
-from ..util.plot import gpplot
-#from ..inference.Expectation_Propagation import FITC
-from ..likelihoods.EP import FITC
-from ..likelihoods import likelihood,probit
-"""
-
-class generalized_FITC(sparse_GP):
-    def __init__(self, X, likelihood, kernel, Z, X_uncertainty=None, Xslices=None,Zslices=None, normalize_X=False):
-    #def __init__(self, X, likelihood, kernel=None, inducing=10, epsilon_ep=1e-3, powerep=[1.,1.]):
-        """
-        Naish-Guzman, A. and Holden, S. (2008) implemantation of EP with FITC.
-
-        :param X: input observations
-        :param likelihood: Output's likelihood (likelihood class)
-        :param kernel: a GPy kernel
-        :param Z:  Either an array specifying the inducing points location or a scalar defining their number.
-        """
-
-        if type(Z) == int:
-            self.M = Z
-            self.Z = (np.random.random_sample(self.D*self.M)*(self.X.max()-self.X.min())+self.X.min()).reshape(self.M,-1)
-        elif type(Z) == np.ndarray:
-            self.Z = Z
-            self.M = self.Z.shape[0]
-
-        self._precision = likelihood.precision
-
-        sparse_GP.__init__(self, X, likelihood, kernel=kernel, Z=self.Z, X_uncertainty=None, Xslices=None,Zslices=None, normalize_X=False)
-        self.scale_factor = 100.
-
-    def update_likelihood_approximation(self):
-        """
-        Approximates a non-gaussian likelihood using Expectation Propagation
-
-        For a Gaussian (or direct: TODO) likelihood, no iteration is required:
-        this function does nothing
-        """
-        if self.has_uncertain_inputs:
-            raise NotImplementedError, "FITC approximation not implemented for uncertain inputs"
-        else:
-            self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
-            self._precision = self.likelihood.precision # Save the true precision
-            self.likelihood.precision = self.likelihood.precision/(1. + self.likelihood.precision*self.Diag0[:,None]) # Add the diagonal element of the FITC approximation
-            self._set_params(self._get_params()) # update the GP
-
-    def _set_params(self, p):
-        self.Z = p[:self.M*self.Q].reshape(self.M, self.Q)
-        self.kern._set_params(p[self.Z.size:self.Z.size+self.kern.Nparam])
-        self.likelihood._set_params(p[self.Z.size+self.kern.Nparam:])
-        self._compute_kernel_matrices()
-        self._computations()
-        self._FITC_computations()
-
-    def _FITC_computations(self):
-        """
-        FITC approximation doesn't have the correction term in the log-likelihood bound,
-        but adds a diagonal term to the covariance matrix.
-        This function:
-            - computes the diagonal term
-            - eliminates the extra terms computed in the sparse_GP approximation
-            - computes the likelihood gradients wrt the true precision.
-        """
-        # Compute FITC's diagonal term of the covariance
-        sf = self.scale_factor
-        sf2 = sf**2
-        self.Qnn = mdot(self.psi1.T,self.Kmmi,self.psi1)
-        self.Diag0 = self.psi0 - np.diag(self.Qnn)
-
-        self.Diag = self.Diag0/(1.+ self.Diag0 * self._precision.flatten())
-        self.P = (self.Diag / self.Diag0)[:,None] * self.psi1.T
-        self.RPT0 = np.dot(self.Lmi,self.psi1)
-        self.L = np.linalg.cholesky(np.eye(self.M) + np.dot(self.RPT0,(1./self.Diag0 - self.Diag/(self.Diag0**2))[:,None]*self.RPT0.T))
-        self.R,info = linalg.flapack.dtrtrs(self.L,self.Lmi,lower=1)
-        self.RPT = np.dot(self.R,self.P.T)
-        self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT)
-        self.w = self.Diag * self.likelihood.v_tilde
-        self.gamma = np.dot(self.R.T, np.dot(self.RPT,self.likelihood.v_tilde))
-        self.mu = self.w + np.dot(self.P,self.gamma)
-        self.mu_tilde = (self.likelihood.v_tilde/self.likelihood.tau_tilde)[:,None]
-
-        # Remove extra term from dL_dpsi
-        self.dL_dpsi0 = np.zeros(self.N)
-        # Remove extra term from dL_dKmm
-        self.dL_dKmm = +0.5 * self.D * mdot(self.Lmi.T, self.A, self.Lmi)*sf2 # dB
-        #the partial derivative vector for the likelihood with the true precision
-        if self.likelihood.Nparams ==0:
-            #save computation here
-            self.partial_for_likelihood = None
-        elif self.likelihood.is_heteroscedastic:
-            raise NotImplementedError, "heteroscedatic derivates not implemented"
-        else:
-            beta = self.likelihood._precision # NOTE the true precison is now '_precison' not 'precision'
-            dbeta =   0.5 * self.N*self.D/beta - 0.5 * np.sum(np.square(self.likelihood.Y))
-            #dbeta += - 0.5 * self.D * (self.psi0.sum() - np.trace(self.A)/beta*sf2)
-            dbeta += - 0.5 * self.D * np.sum(self.Bi*self.A)/beta
-            dbeta += np.sum((self.C - 0.5 * mdot(self.C,self.psi2_beta_scaled,self.C) ) * self.psi1VVpsi1 )/beta
-            self.partial_for_likelihood = -dbeta*self.likelihood.precision**2
-
-
-
-
-    def _raw_predict(self, Xnew, slices, full_cov=True):
-        """
-        Make a prediction for the vsGP model
-
-        Arguments
-        ---------
-        X : Input prediction data - Nx1 numpy array (floats)
-        """
-        Kx = self.kern.K(self.Z, Xnew)
-        #K_x = self.kernel.K(self.Z,X)
-        if full_cov:
-            Kxx = self.kern.K(Xnew)
-        else:
-            Kxx = self.kern.K(Xnew)#FIXME
-            #raise NotImplementedError
-            #Kxx = self.kern.Kdiag(Xnew)
-
-        # q(u|f) = N(u| R0i*mu_u*f, R0i*C*R0i.T)
-
-        # Ci = I + (RPT0)Di(RPT0).T
-        # C = I - [RPT0] * (D+[RPT0].T*[RPT0])^-1*[RPT0].T
-        #   = I - [RPT0] * (D + self.Qnn)^-1 * [RPT0].T
-        #   = I - [RPT0] * (U*U.T)^-1 * [RPT0].T
-        #   = I - V.T * V
-        U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn)
-        V,info = linalg.flapack.dtrtrs(U,self.RPT0.T,lower=1)
-        C = np.eye(self.M) - np.dot(V.T,V)
-        mu_u = np.dot(C,self.RPT0)*(1./self.Diag0[None,:])
-        #self.C = C
-        #self.RPT0 = np.dot(self.R0,self.Knm.T) P0.T
-        #self.mu_u = mu_u
-        #self.U = U
-        # q(u|y) = N(u| R0i*mu_H,R0i*Sigma_H*R0i.T)
-        mu_H = np.dot(mu_u,self.mu)
-        self.mu_H = mu_H
-        Sigma_H = C + np.dot(mu_u,np.dot(self.Sigma,mu_u.T))
-        # q(f_star|y) = N(f_star|mu_star,sigma2_star)
-        KR0T = np.dot(Kx.T,self.Lmi.T)
-        mu_star = np.dot(KR0T,mu_H)
-        sigma2_star = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T))
-        vdiag = np.diag(sigma2_star)
-        return mu_star[:,None],vdiag[:,None]