2

Numpy コードを高速化しようとしていますが、コードがほとんどの時間を C で費やしている特定の関数を実装したいと決めました。

私は実際にはCの新人ですが、行列のすべての行を正規化して合計が1になる関数を書くことができました.コンパイルでき、(Cで)いくつかのデータでテストしましたが、それは私が望むことを行います. その時点で、私は自分自身をとても誇りに思っていました。

今、私は 2d-Numpy 配列を受け入れる必要がある Python から素晴らしい関数を呼び出そうとしています。

私が試したさまざまなことは、

  • 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 64ビットで）次のようにします。

> 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 を見つけられないようですが、同じディレクトリにあり、SWIG が見つからないと文句を言わないnumpy.iため、その理由がわかりません。numpy.i

ctypes の場合、私はそれほど遠くまでは行きませんでしたが、2 次元配列を渡して結果を返す方法がわからなかったので、すぐにドキュメントで迷子になりました。

それで、誰かが私の関数を Python/Numpy で利用できるようにする魔法のトリックを教えてくれませんか?

4

5 に答える 5

8

本当に正当な理由がない限り、C と python のインターフェイスとして cython を使用する必要があります。(numpy/scipy 内で raw C の代わりに cython を使用し始めています)。

私の scikits トークボックスで簡単な例を見ることができます(それ以来 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 では、いくつかのレベルの安全性と速度のトレードオフが利用可能です)。

于 2010-12-01T02:31:53.907 に答える
3

本当の質問に答えるには: SWIG はタイプマップが見つからないとは言いません。typemap を適用できないことを示しています(int DIM1,int DIM2,DATA_TYPE *INPLACE_ARRAY2)。これは、 typemap が定義されていないためDATA_TYPE *です。に適用することを伝える必要がありますdouble*。

%apply ( int DIM1, int DIM2, double* INPLACE_ARRAY2 ) 
       {(size_t nrow, size_t ncol, double* mat)};
于 2011-05-30T18:58:16.833 に答える
2

まず、可能な限り最速の numpy コードを書いていたと確信していますか? 正規化とは、行全体をその合計で割ることを意味する場合、次のような高速なベクトル化コードを記述できます。

matrix /= matrix.sum(axis=0)

これが意図していたものではなく、それでも高速な C 拡張が必要であると確信している場合は、Cではなくcythonで記述することを強くお勧めします。これにより、コードをラッピングする際のすべてのオーバーヘッドと困難が節約され、 Python コードのように見えますが、ほとんどの状況で C と同じくらい速く実行できるものを作成する必要があります。

于 2010-12-01T02:31:05.803 に答える
2

少しの Cython は学ぶ価値があるという他の人たちの意見に同意します。ただし、C または C++ を記述する必要がある場合は、次のように、2 次元をオーバーレイする 1 次元配列を使用します。

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

11 月 8 日に追加: ご存じのとおり、numpy.reshapeは numpy の 2 次元配列を に渡す 1 次元ビューでオーバーレイできます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
于 2010-12-06T11:53:47.203 に答える
1

SciPy には、配列のサンプル コードを含む拡張チュートリアルがあります。 http://docs.scipy.org/doc/numpy/user/c-info.how-to-extend.html

于 2010-12-01T00:38:10.743 に答える