tests passing (cython)

GPy/kern/_src/stationary.py

@@ -226,9 +226,8 @@ class Stationary(Kern):
         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
+        X, X2 = np.ascontiguousarray(X), np.ascontiguousarray(X2)
+        grad = np.zeros(X.shape)
         stationary_cython.grad_X(X.shape[0], X.shape[1], X2.shape[0], X, X2, tmp, grad)
         return grad/self.lengthscale**2
 

GPy/kern/_src/stationary_utils.c

@@ -1,7 +1,7 @@
 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)
+//#pragma omp parallel for private(n,d, retnd, m)
 for(d=0;d<D;d++){
   for(n=0;n<N;n++){
     retnd = 0.0;

GPy/util/choleskies_cython.pyx

@@ -11,7 +11,6 @@ def flat_to_triang(np.ndarray[double, ndim=2] flat, int M):
     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]

GPy/util/linalg.py

@@ -15,11 +15,7 @@ import warnings
 import os
 from .config import config
 import logging
-
-try:
-    from scipy import weave
-except ImportError:
-    config.set('weave', 'working', 'False')
+import linalg_cython
 
 
 _scipyversion = np.float64((scipy.__version__).split('.')[:2])
@@ -422,114 +418,33 @@ def DSYR(*args, **kwargs):
 
 def symmetrify(A, upper=False):
     """
-    Take the square matrix A and make it symmetrical by copting elements from the lower half to the upper
+    Take the square matrix A and make it symmetrical by copting elements from
+    the lower half to the upper
 
     works IN PLACE.
 
-    note: tries to use weave, falls back to a slower numpy version
+    note: tries to use cython, falls back to a slower numpy version
     """
-    if config.getboolean('weave', 'working'):
-        try:
-            symmetrify_weave(A, upper)
-        except:
-            print("\n Weave compilation failed. Falling back to (slower) numpy implementation\n")
-            config.set('weave', 'working', 'False')
-            symmetrify_numpy(A, upper)
+    if config.getboolean('cython', 'working'):
+        _symmetrify_cython(A, upper)
     else:
-        symmetrify_numpy(A, upper)
+        _symmetrify_numpy(A, upper)
 
 
-def symmetrify_weave(A, upper=False):
-    """
-    Take the square matrix A and make it symmetrical by copting elements from the lower half to the upper
+def _symmetrify_cython(A, upper=False):
+    return linalg_cython.symmetrify(A, upper)
 
-    works IN PLACE.
-
-
-    """
-    N, M = A.shape
-    assert N == M
-
-    c_contig_code = """
-    int iN;
-    for (int i=1; i<N; i++){
-      iN = i*N;
-      for (int j=0; j<i; j++){
-        A[i+j*N] = A[iN+j];
-      }
-    }
-    """
-    f_contig_code = """
-    int iN;
-    for (int i=1; i<N; i++){
-      iN = i*N;
-      for (int j=0; j<i; j++){
-        A[iN+j] = A[i+j*N];
-      }
-    }
-    """
-
-    N = int(N) # for safe type casting
-    if A.flags['C_CONTIGUOUS'] and upper:
-        weave.inline(f_contig_code, ['A', 'N'], extra_compile_args=['-O3'])
-    elif A.flags['C_CONTIGUOUS'] and not upper:
-        weave.inline(c_contig_code, ['A', 'N'], extra_compile_args=['-O3'])
-    elif A.flags['F_CONTIGUOUS'] and upper:
-        weave.inline(c_contig_code, ['A', 'N'], extra_compile_args=['-O3'])
-    elif A.flags['F_CONTIGUOUS'] and not upper:
-        weave.inline(f_contig_code, ['A', 'N'], extra_compile_args=['-O3'])
-    else:
-        if upper:
-            tmp = np.tril(A.T)
-        else:
-            tmp = np.tril(A)
-        A[:] = 0.0
-        A += tmp
-        A += np.tril(tmp, -1).T
-
-
-def symmetrify_numpy(A, upper=False):
-    """
-    Force a matrix to be symmetric
-    """
+def _symmetrify_numpy(A, upper=False):
     triu = np.triu_indices_from(A,k=1)
     if upper:
         A.T[triu] = A[triu]
     else:
         A[triu] = A.T[triu]
 
-#This function appears to be unused. It's use of weave makes it problematic
-#Commenting out for now
-#def cholupdate(L, x):
-#    """
-#    update the LOWER cholesky factor of a pd matrix IN PLACE
-#
-#    if L is the lower chol. of K, then this function computes L\_
-#    where L\_ is the lower chol of K + x*x^T
-#    """
-#    support_code = """
-#    #include <math.h>
-#    """
-#    code = """
-#    double r,c,s;
-#    int j,i;
-#    for(j=0; j<N; j++){
-#      r = sqrt(L(j,j)*L(j,j) + x(j)*x(j));
-#      c = r / L(j,j);
-#      s = x(j) / L(j,j);
-#      L(j,j) = r;
-#      for (i=j+1; i<N; i++){
-#        L(i,j) = (L(i,j) + s*x(i))/c;
-#        x(i) = c*x(i) - s*L(i,j);
-#      }
-#    }
-#    """
-#    x = x.copy()
-#    N = x.size
-#    weave.inline(code, support_code=support_code, arg_names=['N', 'L', 'x'], type_converters=weave.converters.blitz)
-
 def backsub_both_sides(L, X, transpose='left'):
-    """ Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky"""
+    """
+    Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky
+    """
     if transpose == 'left':
         tmp, _ = dtrtrs(L, X, lower=1, trans=1)
         return dtrtrs(L, tmp.T, lower=1, trans=1)[0].T

GPy/util/linalg_cython.pyx

@@ -0,0 +1,34 @@
+cimport numpy as np
+from cpython cimport bool
+import cython
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.nonecheck(False)
+def symmetrify(np.ndarray[double, ndim=2] A, bool upper):
+    cdef int N = A.shape[0]
+    if not upper:
+        for i in xrange(N):
+            for j in xrange(i):
+                A[j, i] = A[i, j]
+    else:
+        for j in xrange(N):
+            for i in xrange(j):
+                A[j, i] = A[i, j]
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.nonecheck(False)
+def cholupdate(np.ndarray[double, ndim=1] x, np.ndarray[double, ndim=2] L, int N):
+    cdef double r
+    cdef double c
+    cdef double s
+    for j in xrange(N):
+        r = np.sqrt(L[j,j]*L[j,j] + x[j]*x[j])
+        c = r / L[j,j]
+        s = x[j] / L[j,j]
+        L[j,j] = r
+        for i in xrange(j):
+            L[i,j] = (L[i,j] + s*x[i])/c
+            x[i] = c*x[i] - s*L[i,j];
+        r = np.sqrt(L[j,j])

GPy/util/misc.py

@@ -133,10 +133,11 @@ def kmm_init(X, m = 10):
 ### make a parameter to its corresponding array:
 def param_to_array(*param):
     """
-Convert an arbitrary number of parameters to :class:ndarray class objects. This is for
-converting parameter objects to numpy arrays, when using scipy.weave.inline routine.
-In scipy.weave.blitz there is no automatic array detection (even when the array inherits
-from :class:ndarray)"""
+    Convert an arbitrary number of parameters to :class:ndarray class objects.
+    This is for converting parameter objects to numpy arrays, when using
+    scipy.weave.inline routine.  In scipy.weave.blitz there is no automatic
+    array detection (even when the array inherits from :class:ndarray)
+    """
     import warnings
     warnings.warn("Please use param.values, as this function will be deprecated in the next release.", DeprecationWarning)
     assert len(param) > 0, "At least one parameter needed"

setup.py

@@ -23,6 +23,10 @@ ext_mods = [Extension(name='GPy.kern._src.stationary_cython',
                       sources=['GPy/util/choleskies_cython.c'],
                       include_dirs=[np.get_include()],
                       extra_compile_args=compile_flags),
+            Extension(name='GPy.util.linalg_cython',
+                      sources=['GPy/util/linalg_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'],
                       include_dirs=[np.get_include()],