如何直接从Cython调用numpy/scipy的C函数,避免Python调用开销?

17

我正在尝试在Cython中进行计算,这些计算严重依赖于一些像numpy.log这样的numpy/scipy数学函数。我注意到,如果我在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 函数。如果我将那一行替换为:

This is very expensive, presumably because np.log goes through Python rather than calling the numpy C function directly. If I replace that line with:

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数组)。
  2. 是否有一种类似的方法可以直接从C中调用scipy函数而不需要Python开销?如何导入scipy的C函数库?

具体例子:假设你想在Cython中的for循环内部对标量值调用许多scipy或numpy有用的统计函数(例如scipy.stats.*),重新实现所有这些函数是疯狂的,因此必须调用它们的C版本。例如,与pdf/cdf相关的所有函数以及从各种统计分布进行抽样的函数(例如,请参见http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.pdf.html#scipy.stats.rv_continuous.pdfhttp://www.johndcook.com/distributions_scipy.html)。如果在循环中使用Python开销调用这些函数,它将变得非常慢。

谢谢。


scipy.stats 的概率密度函数等函数主要是用 Python 实现的。你可以通过一次处理多个数字来避免开销。 - pv.
我刚发现这篇文章,因为我有一个非常相似的问题。我已经将我的循环使用cython进行了优化,但在循环内部还有numpy和scipy的调用,所以我没有看到比默认Python更快的速度。考虑到重新编写numpy和scipy函数是不切实际的,这似乎使得Cython的价值主张比我最初希望的要少得多。 - indigoblue
1个回答

3

您无法在numpy数组上应用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,因为它会慢得多。
更新
以下是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

2
但是如果我需要在无法矢量化的代码中编写for循环,并且我想在该循环中使用某些numpy函数(在标量值而非向量上),该怎么办?我需要在Cython中重新实现这些numpy函数吗?这似乎非常愚蠢。我只能选择选项C来完成,但正如您所指出的,这不是一个好选择。 - user248237
在这种情况下,您最好的选择可能是在Cython中重新实现这些函数并使用选项A。如果您不介意我问一下,在循环内部您需要调用哪些NumPy函数? - Bi Rico
2
与统计学相关的各种函数。我编辑了我的答案以提供更多细节,例如所有与概率密度函数/累积分布函数有关的函数以及从各种统计分布中进行抽样(例如,请参见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)。似乎不应该使用libc.math的基元重新实现所有这些pdf,只是为了将其转换为Cython ... 必须有更好的方法吧? - user248237
你可能在numpy中有很多逐元素ufunc的访问权限。 - Mad Physicist

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接