Как вызвать функции numpy/scipy C из Cython напрямую, без накладных вызовов Python?

Я пытаюсь сделать вычисления в Cython, которые в значительной степени полагаются на некоторые математические функции numpy/scipy, такие как numpy.log. Я заметил, что если я повторяю несколько циклов numpy/scipy в цикле в Cython, есть огромные накладные расходы, например:

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

Мои вопросы:

  • Можно ли вызвать функцию C и передать ей массив numpy? Или он может использоваться только для скалярных значений, что потребует от меня написать цикл (например, если я хотел применить его к массиву foo выше).
  • Есть ли аналогичный способ вызова scipy-функций из C напрямую без накладных расходов Python? Каким образом я могу импортировать библиотеку функций Scipy C?

Конкретный пример: скажем, вы хотите называть многие из scipy или numpy полезных статистических функций (например, scipy.stats.*) для скалярных значений внутри цикла for в Cython? Это сумасшедшее, чтобы переопределить все эти функции в 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 в цикле, это будет непомерно медленным.

спасибо.

Ответы

Ответ 1

Вы не можете применять C-функции, такие как log on numpy arrays, а numpy не имеет библиотеки функций C, которую вы можете вызывать из cython.

Функции Numpy уже оптимизированы для вызова массивов numpy. Если у вас нет очень уникального варианта использования, вы не увидите большой выгоды от переопределения функций numpy как функций C. (Возможно, некоторые функции в numpy не реализованы хорошо, в chich случае рассмотрите возможность отправки ваших вложений в виде патчей.) Однако вы действительно хорошо понимаете.

# 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 следует избегать и будет намного медленнее.

Update

Вот код для scipy.stats.norm.pdf, так как вы можете видеть его написанным на питоне с numpy и scipy вызовами. Существует не 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