adding choleskies cythonized

GPy/kern/_src/coregionalize.py

@@ -6,6 +6,7 @@ import numpy as np
 from ...core.parameterization import Param
 from ...core.parameterization.transformations import Logexp
 from ...util.config import config # for assesing whether to use weave
+import coregionalize_cython
 
 try:
     from scipy import weave
@@ -169,9 +170,9 @@ class Coregionalize(Kern):
                 dL_dK_small[j,i] = tmp1[:,index2==j].sum()
         return dL_dK_small
 
-    def gradient_reduce_cython(self, dL_dK, index, index2):
+    def _gradient_reduce_cython(self, dL_dK, index, index2):
         index, index2 = index[:,0], index2[:,0]
-        return coregionalize_cython.gradient_reduce(self.output_dim, dL_dK, index, index2
+        return coregionalize_cython.gradient_reduce(self.B.shape[0], dL_dK, index, index2)
 
 
     def update_gradients_diag(self, dL_dKdiag, X):

GPy/kern/_src/coregionalize_cython.pyx

@@ -1,40 +1,33 @@
+#cython: boundscheck=False
+#cython: wraparound=False
 import cython
 import numpy as np
 cimport numpy as np
 
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nonecheck(False)
 def K_symmetric(np.ndarray[double, ndim=2] B, np.ndarray[int, ndim=1] X):
-    N = X.size
+    cdef int N = X.size
-    K = np.zeros((N, N))
+    cdef np.ndarray[np.double_t, ndim=2] K = np.empty((N, N))
     for n in range(N):
         for m in range(N):
             K[n,m] = B[X[n],X[m]]
     return K
 
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nonecheck(False)
 def K_asymmetric(np.ndarray[double, ndim=2] B, np.ndarray[int, ndim=1] X, np.ndarray[int, ndim=1] X2):
-    N = X.size
+    cdef int N = X.size
-    M = X2.size
+    cdef int M = X2.size
-    K = np.zeros((N, M))
+    cdef np.ndarray[np.double_t, ndim=2] K = np.empty((N, M))
     for n in range(N):
         for m in range(M):
             K[n,m] = B[X[n],X2[m]]
     return K
 
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nonecheck(False)
 def gradient_reduce(int D, np.ndarray[double, ndim=2] dL_dK, np.ndarray[int, ndim=1] index, np.ndarray[int, ndim=1] index2):
-        dL_dK_small = np.zeros((D, D))
+        cdef np.ndarray[np.double_t, ndim=2] dL_dK_small = np.zeros((D, D))
-        N = index.size
+        cdef int N = index.size
-        M = index2.size
+        cdef int M = index2.size
         for i in range(M):
             for j in range(N):
-                dL_dK_small[index[j] + D*index2[i]] += dL_dK[i+j*M];
+                dL_dK_small[index[j],index2[i]] += dL_dK[i,j];
         return dL_dK_small
 
 

GPy/util/choleskies_cython.pyx

@@ -0,0 +1,42 @@
+# Copyright James Hensman and Alan Saul 2015
+
+import numpy as np
+cimport numpy as np
+
+from . import linalg
+
+def flat_to_triang(np.ndarray[double, ndim=3] flat, int M):
+    """take a matrix N x D and return a M X M x D array where
+
+    N = M(M+1)/2
+
+    the lower triangluar portion of the d'th slice of the result is filled by the d'th column of flat.
+    This is the weave implementation
+    """
+    cdef int N = flat.shape[0]
+    cdef int D = flat.shape[1]
+    cdef int count = 0
+    cdef np.ndarray[double, ndim=2] ret = np.empty((M, M, D))
+    for d in range(D):
+        count = 0
+        for m in range(M):
+            for mm in range(m):
+                ret[m, mm, d] = flat[count,d]
+                count += 1
+    return ret
+
+def triang_to_flat(np.ndarray[double, ndim=3] L):
+    cdef int M = L.shape[0]
+    cdef int D = L.shape[2]
+    cdef int N = M*(M+1)/2
+    cdef int count = 0
+    cdef np.ndarray[double, ndim=2] flat = np.empty((N, D))
+    for d in range(D):
+        count = 0
+        for m in range(M):
+            for mm in range(m):
+                flat[count,d] = L[m, mm, d]
+                count += 1
+    return flat
+
+

setup.py

@@ -18,8 +18,12 @@ ext_mods = [Extension(name='GPy.kern._src.stationary_cython',
                       include_dirs=[np.get_include()],
                       extra_compile_args=compile_flags,
                       extra_link_args = ['-lgomp']),
+            Extension(name='GPy.util.choleskies_cython',
+                      sources=['GPy/util/choleskies_cython.c'],
+                      include_dirs=[np.get_include()],
+                      extra_compile_args=compile_flags),
             Extension(name='GPy.kern._src.coregionalize_cython',
-                      sources=['GPy/kern/_src/coregionalize_cython.c','GPy/kern/_src/coregionalize_cython.c'],
+                      sources=['GPy/kern/_src/coregionalize_cython.c'],
                       include_dirs=[np.get_include()],
                       extra_compile_args=compile_flags)]