From 7d9352c7330d9c826c21c9e8f8cb4aee930037b5 Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Sun, 28 Apr 2013 21:37:36 +0100
Subject: [PATCH] weaved coregionalise. much performance gained

---
 GPy/kern/coregionalise.py   | 62 ++++++++++++++++++++++++++++++++++---
 GPy/kern/kern.py            |  1 +
 GPy/kern/prod.py            |  9 ++++--
 GPy/kern/prod_orthogonal.py |  9 ++++--
 4 files changed, 70 insertions(+), 11 deletions(-)

diff --git a/GPy/kern/coregionalise.py b/GPy/kern/coregionalise.py
index a76bb31e..a4d22c2d 100644
--- a/GPy/kern/coregionalise.py
+++ b/GPy/kern/coregionalise.py
@@ -5,10 +5,11 @@ from kernpart import kernpart
 import numpy as np
 from GPy.util.linalg import mdot, pdinv
 import pdb
+from scipy import weave
 
 class coregionalise(kernpart):
     """
-    Kernel for Intrisec Corregionalization Models
+    Kernel for Intrinsic Corregionalization Models
     """
     def __init__(self,Nout,R=1, W=None, kappa=None):
         self.D = 1
@@ -42,19 +43,70 @@ class coregionalise(kernpart):
 
     def K(self,index,index2,target):
         index = np.asarray(index,dtype=np.int)
+
+        #here's the old code (numpy)
+        #if index2 is None:
+            #index2 = index
+        #else:
+            #index2 = np.asarray(index2,dtype=np.int)
+        #false_target = target.copy()
+        #ii,jj = np.meshgrid(index,index2)
+        #ii,jj = ii.T, jj.T
+        #false_target += self.B[ii,jj]
+
         if index2 is None:
-            index2 = index
+            code="""
+            for(int i=0;i<N; i++){
+              target[i+i*N] += B[index[i]+Nout*index[i]];
+              for(int j=0; j<i; j++){
+                  target[j+i*N] += B[index[i]+Nout*index[j]];
+                  target[i+j*N] += target[j+i*N];
+                }
+              }
+            """
+            N,B,Nout = index.size, self.B, self.Nout
+            weave.inline(code,['target','index','N','B','Nout'])
         else:
             index2 = np.asarray(index2,dtype=np.int)
-        ii,jj = np.meshgrid(index,index2)
-        ii,jj = ii.T, jj.T
-        target += self.B[ii,jj]
+            code="""
+            for(int i=0;i<M; i++){
+              for(int j=0; j<N; j++){
+                  target[i+j*M] += B[Nout*index[j]+index2[i]];
+                }
+              }
+            """
+            N,M,B,Nout = index.size,index2.size, self.B, self.Nout
+            weave.inline(code,['target','index','index2','N','M','B','Nout'])
+
 
     def Kdiag(self,index,target):
         target += np.diag(self.B)[np.asarray(index,dtype=np.int).flatten()]
 
     def dK_dtheta(self,dL_dK,index,index2,target):
         index = np.asarray(index,dtype=np.int)
+        dL_dK_small = np.zeros_like(self.B)
+        if index2 is None:
+            index2 = index
+        else:
+            index2 = np.asarray(index2,dtype=np.int)
+
+        code="""
+        for(int i=0; i<M; i++){
+          for(int j=0; j<N; j++){
+            dL_dK_small[index[j] + Nout*index2[i]] += dL_dK[i+j*M];
+          }
+        }
+        """
+        N, M, Nout = index.size, index2.size, self.Nout
+        weave.inline(code, ['N','M','Nout','dL_dK','dL_dK_small','index','index2'])
+
+        dkappa = np.diag(dL_dK_small)
+        dL_dK_small += dL_dK_small.T
+        dW = (self.W[:,None,:]*dL_dK_small[:,:,None]).sum(0)
+
+        target += np.hstack([dW.flatten(),dkappa])
+
+    def dK_dtheta_old(self,dL_dK,index,index2,target):
         if index2 is None:
             index2 = index
         else:
diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py
index 4547fadc..67333765 100644
--- a/GPy/kern/kern.py
+++ b/GPy/kern/kern.py
@@ -9,6 +9,7 @@ from kernpart import kernpart
 import itertools
 from prod_orthogonal import prod_orthogonal
 from prod import prod
+from ..util.linalg import symmetrify
 
 class kern(parameterised):
     def __init__(self, D, parts=[], input_slices=None):
diff --git a/GPy/kern/prod.py b/GPy/kern/prod.py
index 6a59c220..c16d6034 100644
--- a/GPy/kern/prod.py
+++ b/GPy/kern/prod.py
@@ -40,9 +40,12 @@ class prod(kernpart):
 
     def K(self,X,X2,target):
         """Compute the covariance matrix between X and X2."""
-        if X2 is None: X2 = X
-        target1 = np.zeros((X.shape[0],X2.shape[0]))
-        target2 = np.zeros((X.shape[0],X2.shape[0]))
+        if X2 is None:
+            target1 = np.zeros((X.shape[0],X2.shape[0]))
+            target2 = np.zeros((X.shape[0],X2.shape[0]))
+        else:
+            target1 = np.zeros((X.shape[0],X.shape[0]))
+            target2 = np.zeros((X.shape[0],X.shape[0]))
         self.k1.K(X,X2,target1)
         self.k2.K(X,X2,target2)
         target += target1 * target2
diff --git a/GPy/kern/prod_orthogonal.py b/GPy/kern/prod_orthogonal.py
index fc349da8..6ba9965f 100644
--- a/GPy/kern/prod_orthogonal.py
+++ b/GPy/kern/prod_orthogonal.py
@@ -39,11 +39,14 @@ class prod_orthogonal(kernpart):
 
     def K(self,X,X2,target):
         """Compute the covariance matrix between X and X2."""
-        if X2 is None: X2 = X
         target1 = np.zeros_like(target)
         target2 = np.zeros_like(target)
-        self.k1.K(X[:,:self.k1.D],X2[:,:self.k1.D],target1)
-        self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2)
+        if X2 is None:
+            self.k1.K(X[:,:self.k1.D],None,target1)
+            self.k2.K(X[:,self.k1.D:],None,target2)
+        else:
+            self.k1.K(X[:,:self.k1.D],X2[:,:self.k1.D],target1)
+            self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2)
         target += target1 * target2
 
     def dK_dtheta(self,dL_dK,X,X2,target):