initial cython commit

GPy/kern/_src/stationary.py

@@ -11,6 +11,7 @@ import numpy as np
 from scipy import integrate
 from ...util.config import config # for assesing whether to use weave
 from ...util.caching import Cache_this
+import stationary_cython
 
 try:
     from scipy import weave
@@ -245,6 +246,20 @@ class Stationary(Kern):
 
         return ret
 
+    def gradients_X_cython(self, dL_dK, X, X2=None):
+        invdist = self._inv_dist(X, X2)
+        dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
+        tmp = invdist*dL_dr
+        if X2 is None:
+            tmp = tmp + tmp.T
+            X2 = X
+        grad = np.zeros_like(X)
+        if hasattr(X, 'values'):X = X.values #remove the GPy wrapping to make passing into weave safe
+        if hasattr(X2, 'values'):X2 = X2.values
+        stationary_cython.grad_X(X.shape[0], X.shape[1], X2.shape[0], X, X2, tmp, grad)
+        return grad/self.lengthscale**2
+
+
     def gradients_X_weave(self, dL_dK, X, X2=None):
         invdist = self._inv_dist(X, X2)
         dL_dr = self.dK_dr_via_X(X, X2) * dL_dK

GPy/kern/_src/stationary_cython.pyx

@@ -0,0 +1,20 @@
+#cython: boundscheck=False
+#cython: wraparound=False
+import numpy as np
+cimport numpy as np
+
+ctypedef np.float64_t DTYPE_t
+ 
+cdef extern from "stationary_utils.h":
+    void _grad_X "_grad_X" (int N, int D, int M, double* X, double* X2, double* tmp, double* grad)
+ 
+def grad_X(int N, int D, int M,
+        np.ndarray[DTYPE_t, ndim=2] _X,
+        np.ndarray[DTYPE_t, ndim=2] _X2,
+        np.ndarray[DTYPE_t, ndim=2] _tmp,
+        np.ndarray[DTYPE_t, ndim=2] _grad):
+    cdef double *X = <double*> _X.data
+    cdef double *X2 = <double*> _X2.data
+    cdef double *tmp = <double*> _tmp.data
+    cdef double *grad = <double*> _grad.data
+    _grad_X(N, D, M, X, X2, tmp, grad) # return nothing, work in place.

GPy/kern/_src/stationary_utils.c

@@ -0,0 +1,20 @@
+void _grad_X(int N, int D, int M, double* X, double* X2, double* tmp, double* grad){
+int n,m,d;
+double retnd;
+#pragma omp parallel for private(n,d, retnd, m)
+for(d=0;d<D;d++){
+  for(n=0;n<N;n++){
+    retnd = 0.0;
+    for(m=0;m<M;m++){
+      retnd += tmp[n*M+m]*(X[n*D+d]-X2[m*D+d]);
+    }
+    grad[n*D+d] = retnd;
+  }
+}
+} //grad_X
+
+//#weave_options = {'headers'           : ['<omp.h>'],
+         //'extra_compile_args': ['-fopenmp -O3'], # -march=native'],
+         //'extra_link_args'   : ['-lgomp']}
+
+

GPy/kern/_src/stationary_utils.h

@@ -0,0 +1,2 @@
+#include <omp.h>
+void _grad_X(int N, int D, int M, double*X, double* X2, double* tmp, double* grad);