From d36ba011ef33102d9c9daf304c0d3565409a8fd7 Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Fri, 19 Jul 2013 10:27:58 +0100
Subject: [PATCH] changes to psi2 in linear

---
 GPy/kern/parts/linear.py | 45 ++++++++++++++++++++++------------------
 1 file changed, 25 insertions(+), 20 deletions(-)

diff --git a/GPy/kern/parts/linear.py b/GPy/kern/parts/linear.py
index 04bd78a4..82a8c263 100644
--- a/GPy/kern/parts/linear.py
+++ b/GPy/kern/parts/linear.py
@@ -140,30 +140,26 @@ class Linear(Kernpart):
     def dpsi1_dZ(self, dL_dpsi1, Z, mu, S, target):
         self.dK_dX(dL_dpsi1.T, Z, mu, target)
 
-    def psi2(self, Z, mu, S, target):
-        """
-        returns N,num_inducing,num_inducing matrix
-        """
+    def psi2_old(self, Z, mu, S, target):
         self._psi_computations(Z, mu, S)
-#         psi2_old = self.ZZ * np.square(self.variances) * self.mu2_S[:, None, None, :]
-#         target += psi2.sum(-1)
-        # slow way of doing it, but right
-#         psi2_real = rm np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
-#         for n in range(mu.shape[0]):
-#             for m_prime in range(Z.shape[0]):
-#                 for m in range(Z.shape[0]):
-#                     tmp = self._Z[m:m + 1] * self.variances
-#                     tmp = np.dot(tmp, (tdot(self._mu[n:n + 1].T) + np.diag(S[n])))
-#                     psi2_real[n, m, m_prime] = np.dot(tmp, (
-#                             self._Z[m_prime:m_prime + 1] * self.variances).T)
-#         mu2_S = (self._mu[:, None, :] * self._mu[:, :, None])
-#         mu2_S[:, np.arange(self.input_dim), np.arange(self.input_dim)] += self._S
-#         psi2 = (self.ZA[None, :, None, :] * mu2_S[:, None]).sum(-1)
-#         psi2 = (psi2[:, :, None] * self.ZA[None, None]).sum(-1)
-#         psi2_tensor = np.tensordot(self.ZZ[None, :, :, :] * np.square(self.variances), self.mu2_S[:, None, None, :], ((3), (3))).squeeze().T
         target += self._psi2
 
+    def psi2(self,Z,mu,S,target):
+        tmp = np.zeros((mu.shape[0], Z.shape[0]))
+        self.K(mu,Z,tmp)
+        target += tmp[:,:,None]*tmp[:,None,:] + np.sum(S[:,None,None,:]*self.variances**2*Z[None,:,None,:]*Z[None,None,:,:],-1)
+
     def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, target):
+        tmp = np.zeros((mu.shape[0], Z.shape[0]))
+        self.K(mu,Z,tmp)
+        self.dK_dtheta(2.*np.sum(dL_dpsi2*tmp[:,None,:],2),mu,Z,target)
+        result= 2.*(dL_dpsi2[:,:,:,None]*S[:,None,None,:]*self.variances*Z[None,:,None,:]*Z[None,None,:,:]).sum(0).sum(0).sum(0)
+        if self.ARD:
+            target += result.sum(0).sum(0).sum(0)
+        else:
+            target += result.sum()
+
+    def dpsi2_dtheta_old(self, dL_dpsi2, Z, mu, S, target):
         self._psi_computations(Z, mu, S)
         tmp = dL_dpsi2[:, :, :, None] * (self.ZAinner[:, :, None, :] * (2 * Z)[None, None, :, :])
         if self.ARD:
@@ -172,6 +168,15 @@ class Linear(Kernpart):
             target += tmp.sum()
 
     def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
+        tmp = np.zeros((mu.shape[0], Z.shape[0]))
+        self.K(mu,Z,tmp)
+        self.dK_dX(2.*np.sum(dL_dpsi2*tmp[:,None,:],2),mu,Z,target_mu)
+
+        Zs = Z*self.variances
+        Zs_sq = Zs[:,None,:]*Zs[None,:,:]
+        target_S += (dL_dpsi2[:,:,:,None]*Zs_sq[None,:,:,:]).sum(1).sum(1)
+
+    def dpsi2_dmuS_old(self, dL_dpsi2, Z, mu, S, target_mu, target_S):
         """Think N,num_inducing,num_inducing,input_dim """
         self._psi_computations(Z, mu, S)
         AZZA = self.ZA.T[:, None, :, None] * self.ZA[None, :, None, :]