如何直接从Cython调用numpy/scipy C函数，而不需要Python调用开销

Question

我正试图在Cython中进行计算，这些计算严重依赖于一些numpy/scipy数学函数，如numpy.log。我注意到，如果我在Cython的循环中反复调用numpy/scipy函数，会有巨大的开销，例如:

import numpy as np
cimport numpy as np
np.import_array()
cimport cython
def myloop(int num_elts):
   cdef double value = 0
   for n in xrange(num_elts):
     # call numpy function
     value = np.log(2)

这是非常昂贵的，可能是因为np.log通过Python而不是直接调用numpy C函数。如果我用

替换这一行

from libc.math cimport log
...
# calling libc function 'log'
value = log(2)

就快多了。但是，当我尝试将numpy数组传递给libc.math.log:

时

cdef np.ndarray[long, ndim=1] foo = np.array([1, 2, 3])
log(foo)

给出如下错误:

TypeError: only length-1 arrays can be converted to Python scalars

我的问题是:

1. 是否有可能调用C函数并传递它一个numpy数组?或者它只能用于标量值，这需要我写一个循环(例如，如果我想把它应用到上面的foo数组)。
有没有类似的方法可以直接从C调用scipy函数而不需要Python开销?我如何导入scipy的C函数库?

具体的例子:假设你想在Cython的for循环中调用许多scipy或numpy的有用统计函数(例如scipy.stats.*)的标量值?在Cython中重新实现所有这些函数是疯狂的，所以必须调用它们的C版本。例如，所有与pdf/cdf相关的函数和来自各种统计分布的采样(例如，参见http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.pdf.html#scipy.stats.rv_continuous.pdf和http://www.johndcook.com/distributions_scipy.html)如果你在循环中调用Python开销的这些函数，它将会非常慢。

谢谢。

Answer 1

您不能在numpy数组上应用log等C函数，并且numpy没有可以从cython调用的C函数库。

Numpy函数已经优化为可以在Numpy数组上调用。除非您有非常独特的用例，否则您不会看到将numpy函数重新实现为C函数的好处。(numpy中的某些函数可能没有很好地实现，在这种情况下，请考虑将导入作为补丁提交。)但是你提出了一个很好的观点。

# A
from libc.math cimport log
for i in range(N):
    r[i] = log(foo[i])
# B
r = np.log(foo)
# C
for i in range(n):
    r[i] = np.log(foo[i])

一般来说，A和B应该有相似的运行时间，但应该避免使用C，因为C会慢得多。

更新

下面是scipy.stats.norm.pdf的代码，正如您所看到的，它是用numpy和scipy调用用python编写的。这段代码没有C版本，你必须调用它"通过python"。如果这是阻碍你的原因，你需要在C/Cython中重新植入它，但首先我会花一些时间非常仔细地分析代码，看看是否有更容易实现的目标。

def pdf(self,x,*args,**kwds):
    loc,scale=map(kwds.get,['loc','scale'])
    args, loc, scale = self._fix_loc_scale(args, loc, scale)
    x,loc,scale = map(asarray,(x,loc,scale))
    args = tuple(map(asarray,args))
    x = asarray((x-loc)*1.0/scale)
    cond0 = self._argcheck(*args) & (scale > 0)
    cond1 = (scale > 0) & (x >= self.a) & (x <= self.b)
    cond = cond0 & cond1
    output = zeros(shape(cond),'d')
    putmask(output,(1-cond0)+np.isnan(x),self.badvalue)
    if any(cond):
        goodargs = argsreduce(cond, *((x,)+args+(scale,)))
        scale, goodargs = goodargs[-1], goodargs[:-1]
        place(output,cond,self._pdf(*goodargs) / scale)
    if output.ndim == 0:
        return output[()]
    return output