Come chiamare direttamente le funzioni C di numpy/scipy da Cython, senza sovraccarico di chiamata Python?

Question

Sto tentando di eseguire calcoli in Cython che si basano in gran parte su alcune funzioni matematiche numpy/scipy come numpy.log. Ho notato che se chiamo NumPy funzioni/SciPy ripetutamente in un ciclo in Cython, ci sono enormi costi fissi, ad esempio:

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) 

Questo è molto costoso, presumibilmente perché np.log attraversa Python piuttosto che chiamare la funzione numpy C direttamente. Se sostituisco la riga con:

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

quindi è molto più veloce. Tuttavia, quando si tenta di passare un array di NumPy a libc.math.log:

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

dà questo errore:

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

Le mie domande sono:

1. E 'possibile chiama la funzione C e passa una matrice numpy? Oppure può essere usato solo su valori scalari, il che richiederebbe di scrivere un ciclo (ad esempio se volessi applicarlo all'array foo sopra.)
2. C'è un modo analogo di chiamare le funzioni scipy da C direttamente senza Python in testa? Come posso importare la libreria delle funzioni C di Scipy?

Esempio concreto: dire che si desidera chiamare molti dei SciPy di ​​o statistiche utili funzioni di NumPy (ad esempio scipy.stats.*) su valori scalari all'interno di un ciclo for in Cython? È assurdo dover reimplementare tutte quelle funzioni in Cython, quindi è necessario chiamare le loro versioni C. Ad esempio, tutte le funzioni relative a pdf/cdf e campionamento da varie distribuzioni statistiche (ad esempio, vedere http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.pdf.html#scipy.stats.rv_continuous.pdf e http://www.johndcook.com/distributions_scipy.html) Se si richiamano queste funzioni con l'overhead di Python in un ciclo, sarà proibitivamente lento.

grazie.

Answer 1

Non è possibile applicare funzioni C come accedere a serie di numeri e numpy non ha una libreria di funzioni C che è possibile chiamare da cython.

Le funzioni di Numpy sono già ottimizzate per essere chiamate su array numpy. A meno che non si abbia un caso d'uso davvero unico, non si vedranno molti benefici dal reimplementare le funzioni di Numpy come funzioni C. (È possibile che alcune funzioni in numpy non siano implementate correttamente, nel caso in cui si consideri la possibilità di inviare le proprie importazioni come patch.) Tuttavia si fa notare un buon punto.

# 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]) 

In generale, A e B dovrebbe avere tempi di esecuzione simili, ma C dovrebbero essere evitati e sarà molto più lento.

Aggiornamento

Ecco il codice per scipy.stats.norm.pdf, come si può vedere è scritto in python con NumPy e chiamate SciPy. Non esiste una versione C di questo codice, devi chiamarlo "attraverso python". Se questo è ciò che ti trattiene, dovrai re-innestarlo in C/Cython, ma prima passerei del tempo a profilare il codice molto attentamente per vedere se ci sono dei frutti in sospensione più bassi.

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