python - Extending Numpy with C function -

i trying speed numpy code , decided wanted implement 1 particular function code spent of time in c.

i'm rookie in c, managed write function normalizes every row in matrix sum 1. can compile , tested data (in c) , want. @ point proud of myself.

now i'm trying call glorious function python should accept 2d-numpy array.

the various things i've tried are

• swig

• swig + numpy.i

• ctypes

my function has prototype

void normalize_logspace_matrix(size_t nrow, size_t ncol, double mat[nrow][ncol]); 

so takes pointer variable-length array , modifies in place.

i tried following pure swig interface file:

%module c_utils  %{ extern void normalize_logspace_matrix(size_t, size_t, double mat[*][*]); %}  extern void normalize_logspace_matrix(size_t, size_t, double** mat); 

then (on mac os x 64bit):

> swig -python c-utils.i  > gcc -fpic c-utils_wrap.c -o c-utils_wrap.o \      -i/library/frameworks/python.framework/versions/6.2/include/python2.6/ \      -l/library/frameworks/python.framework/versions/6.2/lib/python2.6/ -c c-utils_wrap.c: in function ‘_wrap_normalize_logspace_matrix’: c-utils_wrap.c:2867: warning: passing argument 3 of ‘normalize_logspace_matrix’   incompatible pointer type  > g++ -dynamiclib c-utils.o -o _c_utils.so 

in python following error on importing module:

>>> import c_utils traceback (most recent call last):   file "<stdin>", line 1, in <module> importerror: dynamic module not define init function (initc_utils) 

next tried approach using swig + numpy.i:

%module c_utils  %{ #define swig_file_with_init #include "c-utils.h" %} %include "numpy.i" %init %{ import_array(); %}  %apply ( int dim1, int dim2, data_type* inplace_array2 )         {(size_t nrow, size_t ncol, double* mat)};  %include "c-utils.h" 

however, don't further this:

> swig -python c-utils.i c-utils.i:13: warning 453: can't apply (int dim1,int dim2,data_type *inplace_array2). no typemaps defined. 

swig doesn't seem find typemaps defined in numpy.i, don't understand why, because numpy.i in same directory , swig doesn't complain can't find it.

with ctypes didn't far, got lost in docs pretty since couldn't figure out how pass 2d-array , result back.

so show me magic trick how make function available in python/numpy?

unless have reason not to, should use cython interface c , python. (we starting use cython instead of raw c inside numpy/scipy themselves).

you can see simple example in scikits talkbox (since cython has improved quite bit since then, think write better today).

def cslfilter(c_np.ndarray b, c_np.ndarray a, c_np.ndarray x):     """fast version of slfilter set of frames , filter coefficients.     more precisely, given rank 2 arrays coefficients , input,     computes:      in range(x.shape[0]):         y[i] = lfilter(b[i], a[i], x[i])      useful processing on set of windows variable     filters, e.g. compute lpc residual signal chopped set of     windows.      parameters     ----------         b: array             recursive coefficients         a: array             non-recursive coefficients         x: array             signal filter      note     ----      specialized function, , not handle other types     double, nor initial conditions."""      cdef int na, nb, nfr, i, nx     cdef double *raw_x, *raw_a, *raw_b, *raw_y     cdef c_np.ndarray[double, ndim=2] tb     cdef c_np.ndarray[double, ndim=2] ta     cdef c_np.ndarray[double, ndim=2] tx     cdef c_np.ndarray[double, ndim=2] ty      dt = np.common_type(a, b, x)      if not dt == np.float64:         raise valueerror("only float64 supported now")      if not x.ndim == 2:         raise valueerror("only input of rank 2 support")      if not b.ndim == 2:         raise valueerror("only b of rank 2 support")      if not a.ndim == 2:         raise valueerror("only of rank 2 support")      nfr = a.shape[0]     if not nfr == b.shape[0]:         raise valueerror("number of filters should same")      if not nfr == x.shape[0]:         raise valueerror, \               "number of filters , number of frames should same"      tx = np.ascontiguousarray(x, dtype=dt)     ty = np.ones((x.shape[0], x.shape[1]), dt)      na = a.shape[1]     nb = b.shape[1]     nx = x.shape[1]      ta = np.ascontiguousarray(np.copy(a), dtype=dt)     tb = np.ascontiguousarray(np.copy(b), dtype=dt)      raw_x = <double*>tx.data     raw_b = <double*>tb.data     raw_a = <double*>ta.data     raw_y = <double*>ty.data      in range(nfr):         filter_double(raw_b, nb, raw_a, na, raw_x, nx, raw_y)         raw_b += nb         raw_a += na         raw_x += nx         raw_y += nx      return ty 

as can see, besides usual argument checking in python, same thing (filter_double function can written in pure c in separate library if want to). of course, since compiled code, failing check argument crash interpreter instead of raising exception (there several levels of safety vs speed tradeoffs available recent cython, though).