Расширение Numpy с помощью функции C

Question

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

Я на самом деле новичок в 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, похоже, не находит таблицы типов, определенные в numpy.i, но я не понимаю, почему, поскольку numpy.i находится в том же каталоге, и SWIG не жалуется, что не может его найти.

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

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

Answer 1

Если у вас нет веской причины не делать этого, вам следует использовать cython для взаимодействия с C и python.(Мы начинаем использовать cython вместо необработанного 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

Чтобы ответить на реальный вопрос: SWIG не говорит вам, что не может найти никаких типографских карт.Он говорит вам, что не может применить карту типов (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)};

Answer 3

Я согласен с другими, что маленький 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 4

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

matrix /= matrix.sum(axis=0)

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

Answer 5

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