Wie nummy/scipy C-Funktionen direkt von Cython aufgerufen werden, ohne Python-Aufruf-Overhead?

Question

Ich versuche, Berechnungen in Cython, die sich stark auf einige numpy/scipy mathematische Funktionen wie verlassen. Ich habe bemerkt, dass, wenn ich numpy/scipy Funktionen wiederholt in einer Schleife in Cython nennen, gibt es große Gemeinkosten, zum Beispiel:

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) 

Dies ist sehr teuer, vermutlich weil np.log durch Python geht, anstatt die numpy C-Funktion aufrufen direkt. Wenn ich diese Zeile mit ersetzen:

dann ist es viel schneller. Allerdings, wenn ich versuchen, eine numpy Array übergeben libc.math.log:

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

es diesen Fehler gibt:

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

Meine Fragen sind:

1. Ist es möglich, Rufen Sie die C-Funktion auf und übergeben Sie ihr ein numpiges Array? Oder kann es nur für Skalarwerte verwendet werden, für die ich eine Schleife schreiben müsste (zB wenn ich es auf das obige Array anwenden möchte)
2. Gibt es eine analoge Möglichkeit, scipy-Funktionen von C direkt ohne a aufzurufen Python-Aufwand? Wie kann ich scipys C-Funktionsbibliothek importieren?

Konkretes Beispiel: Angenommen, Sie viele scipy oder der numpy des nützlichen Statistikfunktionen (z scipy.stats.*) auf skalare Werte innerhalb einer for Schleife in Cython nennen wollen? Es ist verrückt, alle diese Funktionen in Cython neu zu implementieren, daher müssen ihre C-Versionen aufgerufen werden. Zum Beispiel alle Funktionen im Zusammenhang mit pdf/cdf und Sampling von verschiedenen statistischen Verteilungen (z. B. siehe http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.pdf.html#scipy.stats.rv_continuous.pdf und http://www.johndcook.com/distributions_scipy.html) Wenn Sie diese Funktionen mit Python Overhead in einer Schleife aufrufen, wird es prohibitiv langsam.

danke.

Answer 1

Sie können keine C-Funktionen wie das Anmelden von numpy-Arrays anwenden, und numpy hat keine C-Funktionsbibliothek, die Sie von Cython aus aufrufen können.

Numpy-Funktionen sind bereits so optimiert, dass sie auf numpy Arrays aufgerufen werden können. Wenn Sie nicht einen sehr speziellen Anwendungsfall haben, werden Sie nicht viel Nutzen daraus ziehen, wenn Sie numpy Funktionen als C-Funktionen neu implementieren. (Es ist möglich, dass einige Funktionen in numpy nicht gut implementiert sind. In diesem Fall sollten Sie Ihre Imports als Patches einreichen.) Sie haben jedoch einen guten Punkt.

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

Im Allgemeinen sollten A und B ähnliche Laufzeiten haben, aber C sollte vermieden werden und wird viel langsamer sein.

aktualisieren

Hier ist der Code für scipy.stats.norm.pdf, wie Sie es geschrieben in Python mit numpy und scipy Anrufe sehen. Es gibt keine C-Version dieses Codes, Sie müssen es "durch Python" nennen. Wenn dies Sie zurückhält, müssen Sie es erneut in C/Cython implantieren, aber zuerst würde ich einige Zeit damit verbringen, den Code sorgfältig zu profilieren, um zu sehen, ob es tiefer hängende Früchte gibt, nach denen zuerst gesucht wird.

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