Расширение Numpy с функцией C

Question

Я пытаюсь ускорить свой код Numpy и решил, что я хотел бы реализовать одну конкретную функцию, где мой код провел большую часть времени в С.

Я на самом деле новичок в C , но мне удалось написать функцию, которая нормализует каждую строку матрицы в сумме до 1. Я могу ее скомпилировать, и я проверил ее с некоторыми данными (в C), и он делает то, что я хочу. В этот момент я очень гордился собой.

Теперь я пытаюсь назвать свою славную функцию из Python, где он должен принять массив 2d-Numpy.

Различные вещи, которые я пробовал это

• SWIG

• SWIG + numpy.i

• ctypes

Моя функция имеет прототип

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

Поэтому он принимает указатель на массив переменной длины и изменяет его на месте.

Я попробовал следующий чистый SWIG интерфейс файла:

%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); 

Тогда я сделал бы (на 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’ from incompatible pointer type 

> g++ -dynamiclib c-utils.o -o _c_utils.so 

В Python я тогда получаю следующую ошибку при импорте мой модуль:

>>> import c_utils 
Traceback (most recent call last): 
    File "<stdin>", line 1, in <module> 
ImportError: dynamic module does not define init function (initc_utils) 

Далее я попробовал этот подход, используя 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" 

Однако, я не получаю дальше, чем это:

> swig -python c-utils.i 
c-utils.i:13: Warning 453: Can't apply (int DIM1,int DIM2,DATA_TYPE *INPLACE_ARRAY2). No typemaps are defined. 

SWIG, кажется, не найти typemaps, определенные в numpy.i, но я не понимаю, почему, потому что numpy.i является в том же каталоге, и SWIG не жалуется, что он не может его найти.

С ctypes я не очень далеко, но потерялся в документах довольно быстро, так как не мог понять, как передать ему 2d-массив, а затем вернуть результат.

Так может кто-нибудь показать мне волшебный трюк, как сделать мою функцию доступной в Python/Numpy?

Answer 1

Если у вас нет по-настоящему уважительной причины, вы должны использовать cython для интерфейса C и python. (Мы начинаем использовать cython вместо raw C внутри numpy/scipy).

Вы можете увидеть простой пример в моих scikits talkbox (так как cython немного улучшился с тех пор, я думаю, что вы могли бы написать его лучше сегодня).

def cslfilter(c_np.ndarray b, c_np.ndarray a, c_np.ndarray x): 
    """Fast version of slfilter for a set of frames and filter coefficients. 
    More precisely, given rank 2 arrays for coefficients and input, this 
    computes: 

    for i in range(x.shape[0]): 
     y[i] = lfilter(b[i], a[i], x[i]) 

    This is mostly useful for processing on a set of windows with variable 
    filters, e.g. to compute LPC residual from a signal chopped into a set of 
    windows. 

    Parameters 
    ---------- 
     b: array 
      recursive coefficients 
     a: array 
      non-recursive coefficients 
     x: array 
      signal to filter 

    Note 
    ---- 

    This is a specialized function, and does not handle other types than 
    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 for 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 a of rank 2 support") 

    nfr = a.shape[0] 
    if not nfr == b.shape[0]: 
     raise ValueError("Number of filters should be the same") 

    if not nfr == x.shape[0]: 
     raise ValueError, \ 
       "Number of filters and number of frames should be the 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 

    for i 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 

Как вы можете видеть, кроме обычного аргумента проверки вы могли бы сделать в Python, это почти то же самое (filter_double это функция, которая может быть записана в чистом C в отдельную библиотеку, если вы хотите) , Конечно, поскольку это скомпилированный код, неспособность проверить ваш аргумент приведет к сбою вашего интерпретатора вместо повышения исключения (хотя есть несколько уровней безопасности и компромисс скорости, доступных с недавним cython).

Answer 2

В SciPy есть учебник по расширению с примером кода для массивов. http://docs.scipy.org/doc/numpy/user/c-info.how-to-extend.html

Answer 3

Во-первых, вы уверены, что писали самый быстрый возможный код?Если по нормализуют значит разделить весь ряд от его суммы, то вы можете написать быстро vectorised код, который выглядит примерно так:

matrix /= matrix.sum(axis=0) 

Если это не то, что вы имели в виду, и вы по-прежнему уверены, что вам нужно быстрое расширение C, я бы настоятельно рекомендовал вам написать его в cython вместо C. Это сэкономит вам все накладные расходы и трудности с переносом кода и позволит вам написать что-то, что похоже на код python, но которое можно заставить запускать как быстро как C в большинстве случаев.

Answer 4

Я согласен с другими, что немного Cython стоит изучать. Но если вы должны написать C или C++, используйте 1d массив, который перекрывает 2d, как это:

// sum1rows.cpp: 2d A as 1d A1 
// Unfortunately 
//  void f(int m, int n, double a[m][n]) { ... } 
// is valid c but not c++ . 
// See also 
// http://stackoverflow.com/questions/3959457/high-performance-c-multi-dimensional-arrays 
// http://stackoverflow.com/questions/tagged/multidimensional-array c++ 

#include <stdio.h> 

void sum1(int n, double x[]) // x /= sum(x) 
{ 
    float sum = 0; 
    for(int j = 0; j < n; j ++ ) 
     sum += x[j]; 
    for(int j = 0; j < n; j ++ ) 
     x[j] /= sum; 
} 

void sum1rows(int nrow, int ncol, double A1[]) // 1d A1 == 2d A[nrow][ncol] 
{ 
    for(int j = 0; j < nrow*ncol; j += ncol ) 
     sum1(ncol, &A1[j]); 
} 

int main(int argc, char** argv) 
{ 
    int nrow = 100, ncol = 10; 
    double A[nrow][ncol]; 
    for(int j = 0; j < nrow; j ++) 
    for(int k = 0; k < ncol; k ++) 
     A[j][k] = (j+1) * k; 

    double* A1 = &A[0][0]; // A as 1d array -- bad practice 
    sum1rows(nrow, ncol, A1); 

    for(int j = 0; j < 2; j ++){ 
     for(int k = 0; k < ncol; k ++){ 
      printf("%.2g ", A[j][k]); 
     } 
     printf("\n"); 
    } 
} 

Добавлено 8 ноября: как вы, наверное, знаете, numpy.reshape может наложить Numpy 2d массив с 1d целью перейти к sum1rows, как это:

import numpy as np 
A = np.arange(10).reshape((2,5)) 
A1 = A.reshape(A.size) # a 1d view of A, not a copy 
# sum1rows(2, 5, A1) 
A[1,1] += 10 
print "A:", A 
print "A1:", A1 

Answer 5

чтобы ответить на реальный вопрос: SWIG не говорит вам, что не может найти никаких typemaps. Он сообщает вам, что он не может применить типовую карту (int DIM1,int DIM2,DATA_TYPE *INPLACE_ARRAY2), потому что нет никакой типовой карты, определенной для DATA_TYPE *. Вы должны сказать ему, вы хотите, чтобы применить его к double*:

%apply (int DIM1, int DIM2, double* INPLACE_ARRAY2) 
     {(size_t nrow, size_t ncol, double* mat)};