From fe405271e5a498e118148629b399c1a11646f83a Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Wed, 13 Aug 2014 10:36:54 +0100
Subject: [PATCH] gradients of predictions for Trevor

---
 GPy/core/gp.py              | 22 ++++++++++++++++++++++
 GPy/kern/_src/stationary.py | 18 +++++++++---------
 GPy/models/gplvm.py         |  6 +++---
 3 files changed, 34 insertions(+), 12 deletions(-)

diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index e0f5755c..ba7c38f2 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -148,6 +148,28 @@ class GP(Model):
         m, v = self._raw_predict(X,  full_cov=False)
         return self.likelihood.predictive_quantiles(m, v, quantiles, Y_metadata)
 
+    def predictive_gradients(self, Xnew):
+        """
+        Compute the derivatives of the latent function with respect to X*
+
+        Given a set of points at which to predict X* (size [N*,Q]), compute the
+        derivatives of the mean and variance. Resulting arrays are sized:
+         dmu_dX* -- [N*, Q ,D], where D is the number of output in this GP (usually one).
+         dv_dX*  -- [N*, Q],    (since all outputs have the same variance)
+
+        """
+        dmu_dX = np.empty((Xnew.shape[0],Xnew.shape[1],self.output_dim))
+        for i in range(self.output_dim):
+            dmu_dX[:,:,i] = self.kern.gradients_X(self.posterior.woodbury_vector[:,i:i+1].T, Xnew, self.X)
+
+        # gradients wrt the diagonal part k_{xx}
+        dv_dX = self.kern.gradients_X(np.eye(Xnew.shape[0]), Xnew)
+        #grads wrt 'Schur' part K_{xf}K_{ff}^{-1}K_{fx}
+        alpha = -2.*np.dot(self.kern.K(Xnew, self.X),self.posterior.woodbury_inv)
+        dv_dX += self.kern.gradients_X(alpha, Xnew, self.X)
+        return dmu_dX, dv_dX
+
+
     def posterior_samples_f(self,X,size=10, full_cov=True):
         """
         Samples the posterior GP at the points X.
diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py
index 67014ab5..c0553238 100644
--- a/GPy/kern/_src/stationary.py
+++ b/GPy/kern/_src/stationary.py
@@ -160,20 +160,20 @@ class Stationary(Kern):
         """
         invdist = self._inv_dist(X, X2)
         dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
-        #The high-memory numpy way:
-        #d =  X[:, None, :] - X2[None, :, :]
-        #ret = np.sum((invdist*dL_dr)[:,:,None]*d,1)/self.lengthscale**2
-        #if X2 is None:
-            #ret *= 2.
-
-        #the lower memory way with a loop
         tmp = invdist*dL_dr
-        ret = np.empty(X.shape, dtype=np.float64)
         if X2 is None:
             tmp = tmp + tmp.T
             X2 = X
-        [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)]
+
+        #The high-memory numpy way:
+        #d =  X[:, None, :] - X2[None, :, :]
+        #ret = np.sum(tmp[:,:,None]*d,1)/self.lengthscale**2
+
+        #the lower memory way with a loop
+        ret = np.empty(X.shape, dtype=np.float64)
+        [np.sum(tmp*(X[:,q][:,None]-X2[:,q][None,:]), axis=1, out=ret[:,q]) for q in xrange(self.input_dim)]
         ret /= self.lengthscale**2
+
         return ret
 
     def gradients_X_diag(self, dL_dKdiag, X):
diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py
index 542dcd31..8f5432ba 100644
--- a/GPy/models/gplvm.py
+++ b/GPy/models/gplvm.py
@@ -45,10 +45,10 @@ class GPLVM(GP):
         self.X.gradient = self.kern.gradients_X(self.grad_dict['dL_dK'], self.X, None)
 
     def jacobian(self,X):
-        target = np.zeros((X.shape[0],X.shape[1],self.output_dim))
+        J = np.zeros((X.shape[0],X.shape[1],self.output_dim))
         for i in range(self.output_dim):
-            target[:,:,i]=self.kern.gradients_X(np.dot(self.Ki,self.likelihood.Y[:,i])[None, :],X,self.X)
-        return target
+            J[:,:,i] = self.kern.gradients_X(self.posterior.woodbury_vector[:,i:i+1], X, self.X)
+        return J
 
     def magnification(self,X):
         target=np.zeros(X.shape[0])