scipy - 如何使用 scipy.integrate.quadpack(或 scipy 中的其他 c/fortran)直接作为来自 cython 的 c

Question

Numpy 有一个基本的 pxd，声明它的 c 接口(interface)到 cython。是否有用于 scipy 组件(尤其是 scipy.integrate.quadpack)的 pxd？

或者，有人可以举例说明如何直接链接到 cython 的 scipy 中包含的 c/fortran 函数吗？到目前为止，我一直使用 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) :未找到符号:_DQAGSE\n 引用自:/Users/shauncutts/.pyxbld/lib.macosx-10.7-intel-2.7/factfiber/stat/pmodel/c/meer.so\n 预期位于:平面命名空间\n/Users/shauncutts/.pyxbld/lib.macosx-10.7-intel-2.7/factfiber/stat/pmodel/c/meer.so\n']

所以我认为我的问题集中在实际上将我对 fortran QDAGSE 的引用“指向”scipy.integrate 中的 _quadpack.so。 (N.B. 也尝试了小写版本。)...而不必通过在其中放置 pxd 来破解 scipy。有什么办法可以做到这一点？也许是动态的？ (我可以动态加载 .so 并将适当的函数指针绑定(bind)到全局变量吗？)

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

Answer 1

感谢: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”)

与之前使用 python 回调的代码相比，它使我的速度提高了将近 5 倍，但在其他方面已经尽可能地进行了优化。

对于额外的功劳...我如何将额外的参数偷偷带入 c 回调？我认为没有全局变量是不可能的？或者我可以使用编织？

无论如何——感谢您的关注——这是一个晦涩难懂的示例，但对于从 cython 调用扩展模块中未声明的 c 函数具有更普遍的适用性。如果有人有更优雅的方法，我很想听听。