Estendentesi Numpy con funzione C

Question

Sto cercando di accelerare il mio codice Numpy e ha deciso che volevo implementare una funzione particolare in cui il mio codice trascorso la maggior parte del tempo in C.

In realtà sono un principiante in C, ma sono riuscito a scrivere la funzione che normalizza ogni riga in una matrice di somma a 1. posso compilare e ho provato con alcuni dati (in C) e si fa quello Voglio. A quel punto ero molto orgoglioso di me stesso.

Ora sto provando a chiamare il mio glorioso funzione da Python dove dovrebbe accettare una matrice 2D-Numpy.

Le varie cose che ho provato sono

• SWIG

• SWIG + numpy.i

• ctypes

La mia funzione è il prototipo

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

Quindi ci vuole un puntatore a una matrice di lunghezza variabile e lo modifica in posizione.

Ho provato il seguente file di interfaccia SWIG puro:

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

Poi mi avrebbe fatto (su Mac OS X a 64 bit):

> 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

In Python ho quindi ottenere il seguente errore di importare il mio modulo:

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

Poi ho provato questo approccio utilizzando 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"

Comunque, non ho ricevuto oltre questo:

> 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 non sembra trovare il typemaps definiti numpy.i, ma non capisco il motivo per cui, a causa numpy.i si trova nella stessa directory e SWIG non si lamenta che non può trovarlo.

Con ctypes non ho avuto molto lontano, ma si sono persi nella documentazione abbastanza rapidamente dal momento che non riuscivo a capire come passare un 2d-array e quindi ottenere la schiena risultato.

Quindi qualcuno potrebbe mostrarmi il trucco magico come rendere la mia funzione disponibile in Python / Numpy?

Answer 1

Se non avete una buona ragione per non farlo, si dovrebbe usare Cython per interfacciarsi C e Python. (Si stanno iniziando a utilizzare Cython invece di greggio C all'interno NumPy / SciPy stessi).

È possibile vedere un semplice esempio nei miei scikits talkbox (dal Cython è migliorata un bel po 'da allora , penso che si potrebbe scrivere meglio di oggi).

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

Come si può vedere, oltre al solito argomento controllando fareste in python, è quasi la stessa cosa (filter_double è una funzione che può essere scritto in puro C in una libreria separata, se si desidera). Naturalmente, dal momento che è stato compilato il codice, non riuscendo a controllare il vostro argomento andrà in crash il vostro interprete, invece di sollevare un'eccezione (ci sono diversi livelli di sicurezza vs compromessi velocità disponibili con la recente Cython, però).

Answer 2

Per rispondere alla vera domanda: SWIG non vi dice che non può trovare nessuna typemaps. Ti dice che non può applicare la (int DIM1,int DIM2,DATA_TYPE *INPLACE_ARRAY2) typemap, che è perché non c'è typemap definito per DATA_TYPE *. È necessario dire che si desidera applicare a una double*:

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

Answer 3

In primo luogo, sei sicuro che si stesse scrivendo il codice NumPy più veloce possibile? Se per normalizzare si divide media tutta la fila per la sua somma, allora si può scrivere codice veloce Vectorised che sembra qualcosa di simile:

matrix /= matrix.sum(axis=0)

Se questo non è quello che si aveva in mente e vi sono ancora sicuro che hai bisogno di un'estensione C veloce, che vi consiglio vivamente che si scrive in Cython al posto di C. Questo vi farà risparmiare tutto l'overhead e le difficoltà nel codice avvolgimento e si permettono di scrivere qualcosa che assomiglia codice Python, ma che possono essere fatte a correre veloce come C nella maggior parte dei casi.

Answer 4

Sono d'accordo con gli altri che un po 'di Cython è ben apprendimento vale la pena. Ma se è necessario scrivere C o C ++, utilizzare un array 1D che si sovrappone il 2d, in questo modo:

// 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" );
    }
}

Aggiunto 8 Nov: come probabilmente sapete, numpy.reshape può sovrapporre una matrice 2D NumPy in vista 1d per passare a sum1rows, in questo modo:

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

SciPy ha un tutorial prolunga con codice di esempio per le matrici. http://docs.scipy.org/doc /numpy/user/c-info.how-to-extend.html