как использовать scipy.integrate.quadpack (или другой c / fortran в scipy) напрямую как c из cython - PullRequest
Новая
       4

как использовать scipy.integrate.quadpack (или другой c / fortran в scipy) напрямую как c из cython

3 голосов
/ 03 апреля 2012

У Numpy есть базовый pxd, который объявляет свой интерфейс c к cython.Существует ли такой pxd для компонентов scipy (особенно scipy.integrate.quadpack)?

В качестве альтернативы, может ли кто-нибудь привести пример прямой связи с функциями c / fortran, включенными в scipy из cython?До сих пор я всегда использовал pyximport ... будет ли это работать здесь, или мне придется связывать, используя distutils (или make ..)?

Спасибо!

---- ОБНОВЛЕНИЕ ----

Приведенный ниже код Cython компилируется;тем не менее, я получаю

ImportError: Building module failed: ['ImportError: dlopen(/Users/shauncutts/.pyxbld/lib.macosx-10.7-intel-2.7/factfiber/stat/pmodel/c/meer.so, 2): Symbol not found: _DQAGSE\n Referenced from: /Users/shauncutts/.pyxbld/lib.macosx-10.7-intel-2.7/factfiber/stat/pmodel/c/meer.so\n Expected in: flat namespace\n in /Users/shauncutts/.pyxbld/lib.macosx-10.7-intel-2.7/factfiber/stat/pmodel/c/meer.so\n']

Так что я думаю, что сосредоточил свою проблему на том, чтобы на самом деле «указать» мою ссылку на fortran QDAGSE на _quadpack.so в scipy.integrate.(NB пробовал и строчную версию.) ... без необходимости взламывать scipy, поместив туда pxd.Есть ли способ, которым я могу это сделать?Возможно динамически?(Могу ли я загрузить .so динамически и привязать соответствующий указатель на глобальную переменную?) 

cdef extern from "stdlib.h":
    void free(void* ptr)
    void* malloc(size_t size)
    void* realloc(void* ptr, size_t size)


ctypedef double (*qagfunc)( double* )

cdef extern:
    # in scipy.integrate._quadpack
    cdef void dqagse(
        qagfunc f, double *a, double *b, double *epsabs, double *epsrel,
        int *limit, 
        double *result, double *abserr, int *neval, int *ier,
        double *alist, double *blist, double *rlist, double *elist,
        int *iord, int *last
        )

class QAGError( ValueError ):
    code = None

cdef double qags( 
    qagfunc quad_function, double a, double b,
    double epsabs=1.49e-8, double epsrel=1.49e-8, 
    int limit = 50 
    ):

    '''
    wrapper for QUADPACK quags/quagse
    '''
    cdef double result, abserr
    cdef int neval = 0, ier = 0, last = 0
    cdef int *iord = <int *>malloc( limit * sizeof( double ) )
    cdef double *alist = <double *>malloc( limit * sizeof( double ) )
    cdef double *blist = <double *>malloc( limit * sizeof( double ) )
    cdef double *rlist = <double *>malloc( limit * sizeof( double ) )
    cdef double *elist = <double *>malloc( limit * sizeof( double ) )

    try:
        DQAGSE(
            quad_function, &a, &b, &epsabs, &epsrel, &limit, 
            &result, &abserr, &neval, &ier, 
            alist, blist, rlist, elist, iord, &last);

    finally:
        free( iord )
        free( alist )
        free( blist )
        free( rlist )
        free( elist )

    if ier > 0:
        raise QAGError( QUAGS_IER_MAP[ ier ] )
    return result

1 Ответ

3 голосов
/ 03 апреля 2012

Благодаря: http://codespeak.net/pipermail/cython-dev/2009-October/007239.html

Адаптация, кажется, работает следующее: 

from os import path
import ctypes

from scipy.integrate import quadpack

ctypedef double (*qagfunc)( double* )

ctypedef void (*quadpack_qagse_fp)( 
    qagfunc f, double *a, double *b, double *epsabs, double *epsrel,
    int *limit, 
    double *result, double *abserr, int *neval, int *ier,
    double *alist, double *blist, double *rlist, double *elist,
    int *iord, int *last
    )

quadpack_path = path.dirname( quadpack.__file__ )
quadpack_ext_path = path.join( quadpack_path, '_quadpack.so' ) 
quadpack_lib = ctypes.CDLL( quadpack_ext_path )

cdef quadpack_qagse_fp quadpack_qagse
quadpack_qagse = ( 
    <quadpack_qagse_fp*><size_t> ctypes.addressof( quadpack_lib.dqagse_ ) )[ 0 ]

С этим добавленным кодом, приведенный выше код, кажется, работает ... хорошо, за исключением того, что досадноДостаточно, ответ немного отличается от Python Quad, а также Quadpack dblquad (я оцениваю двойной интеграл).Я предполагаю, что мне либо нужно больше разрешения, либо нужно указать параметр со значением не по умолчанию.(Возможно, он также не будет портировать на Windows? Но это просто означало бы изменение ".so" на ".dll")

Это дает мне чуть менее 5-кратную скорость по сравнению с предыдущим кодом, который использовал обратный вызов Python, нов противном случае был оптимизирован настолько, насколько я мог.

Для дополнительного кредита ... как мне переправить дополнительные параметры в обратный вызов c?Я полагаю, это невозможно без глобальных переменных?или я мог бы использовать плетение?

В любом случае - спасибо за ваше внимание - это неясный пример, но он имеет более общую применимость к вызову необъявленных функций c в модулях расширения из cython.Если у кого-то есть более изящный способ сделать это, я буду рад услышать.

...