Ist eine "für" Schleife notwendig, wenn Elemente des a-Numpy-Vektors vom vorherigen Element abhängig sind?

T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i)) 

Tm und Tau sind numpy Vektoren der gleichen Länge, die zuvor berechnet wurden und der Wunsch ist, einen neuen Vektor T zu erstellen. Das "i" ist nur enthalten, um die Element-Index für das gewünschte zu zeigen.

Ist eine für die Schleife notwendige Schleife?

6 Solutions collect form web for “Ist eine "für" Schleife notwendig, wenn Elemente des a-Numpy-Vektors vom vorherigen Element abhängig sind?”

Sie könnten denken, das würde funktionieren:

 import numpy as np n = len(Tm) t = np.empty(n) t[0] = 0 # or whatever the initial condition is t[1:] = Tm[1:] + (t[0:n-1] - Tm[1:])**(-tau[1:]) 

Aber das tut es nicht: man kann auf diese Weise nicht rekursionieren (da numpy die ganze RHS berechnet und dann dem LHS zuweist).

Also, wenn du nicht mit einer nicht rekursiven Version dieser Formel kommen kannst, steckst du mit einer expliziten Schleife:

 tt = np.empty(n) tt[0] = 0. for i in range(1,n): tt[i] = Tm[i] + (tt[i-1] - Tm[i])**(-tau[i]) 

In manchen Fällen ist es möglich, diese Art von Rekursion zu haben – nämlich wenn es einen NumPy ufunc für die Rekursionsformel gibt, zB

 c = numpy.arange(10.) numpy.add(c[:-1], c[1:], c[1:]) 

Dies berechnet die akkumulativen Summen über c an Ort und Stelle mit dem Ausgangsparameter von numpy.add . Es kann nicht geschrieben werden als

 c[1:] = c[:-1] + c[1:] 

Denn jetzt ist das Ergebnis der Addition eine vorübergehende, die nach der Berechnung in c[1:] kopiert wird.

Die natürlichste Sache, um jetzt zu versuchen ist, Ihre eigenen ufunc zu definieren:

 def f(T, Tm, tau): return Tm + (T - Tm)**(-tau) uf = numpy.frompyfunc(f, 3, 1) 

Aber aus Gründen, die über mich hinausgehen, funktioniert das nicht:

 uf(T[:-1], Tm[1:], tau[1:], T[1:]) 

Anscheinend wird das Ergebnis nicht direkt an T[1:] , sondern nach Beendigung der Operation temporär gespeichert und kopiert. Auch wenn es funktionierte, würde ich hier keine Beschleunigung erwarten, verglichen mit einer gewöhnlichen Schleife, da es für jeden Eintrag eine Python-Funktion anrufen muss.

Wenn du die Pythonschleife wirklich vermeiden willst, musst du wohl für Cython oder Ctypes gehen.

Es ist klar, dass es irgendwo eine Schleife geben muss. Ich denke, deine Frage ist, wie man die Schleife passieren in numpy statt in Python richtig. Wenn das die wahre Frage ist, dann wie über irgendeinen kreativen Gebrauch von numpy.fromiter ?

Ich habe beschlossen, einige Benchmarks des ursprünglichen Problems zu schaffen, weil ich auch etwas rekursiv berechnen musste, das nicht vektorisiert werden konnte. Ich habe abs () für die Exponentenbasis verwendet, weil das Ergebnis undefiniert ist, wenn die Basis negativ ist. Beispiel für die Berechnung durch normales Looping von Numpy Array:

 alen=100000 Tm = np.random.normal(size=alen).astype('float64') tau = np.random.normal(size=alen).astype('float64') def rec_numpy_loop(Tm,tau,alen): T=np.empty(alen) T[0]=0.0 for i in range(1,alen): T[i] = Tm[i] + (abs(T[i-1] - Tm[i]))**(-tau[i]) return T 

Die Ergebnisse sind, dass die Kompilierung einer Funktion als Cython-Modul ist sehr schnell sowohl in der Ausführung (14 mal schneller als jede Python-Implementierung) und das Schreiben des Codes.

Weitere interessante Fakten: Verschiedene Numpy-Implementierungen variieren leicht mit der schnellsten Implementierung mit a.item () und a.itemset () Notation. Merkwürdig, aber anhängend, listet Arbeiten auf Par mit Anhängen an vorzugeordnete Numpy Arrays. Das Spielen mit Memviews in Cython gab keinen signifikanten Schub für Code-Performance. Fortran Code war sehr kurz, aber fast unmöglich zu debuggen und am Ende Cython erschien schneller als Fortran mit f2py, wenn auch geringfügig. C-Erweiterung hat viel Zeit zu schreiben, weil es zu viel Kesselplatte hat und am Ende funktioniert es mit der gleichen Geschwindigkeit wie Cython. Pure C mit reinen C Arrays war zweimal schneller als alles, was von Python genannt wurde.

Ich bin kein Pro von C / C ++, vielleicht kann man ein schnelleres Programm schreiben.

 Looping over Numpy array: In [57]: %timeit -o rec_numpy_loop(Tm,tau,alen) 10 loops, best of 3: 87.9 ms per loop Using a.item() and a.itemset() with Numpy: In [58]: %timeit -o rec_numpy_loop_item(Tm,tau,alen) 10 loops, best of 3: 74.3 ms per loop Using np.fromiter() from Numpy: In [59]: %timeit -o rec_numpy_iter(Tm,tau,alen) 10 loops, best of 3: 78.1 ms per loop Using Numpy arrays with data but appending to a List: In [60]: %timeit -o rec_py_loop(Tm,tau,alen) 10 loops, best of 3: 91 ms per loop Calling a function compiled to Cython: In [61]: %timeit -o rec_cy_loop(Tm,tau,alen) 100 loops, best of 3: 6.46 ms per loop Using Memviews in Cython: In [62]: %timeit -o rec_cy_loop_memviews(Tm,tau,alen) 100 loops, best of 3: 6.26 ms per loop Using Memviews both for looping and as input variables in Cython: In [63]: %timeit -o rec_cy_loop_memviews_w_input(Tm,tau,alen) 100 loops, best of 3: 6.53 ms per loop Calling a fortran function compiled using f2py: In [64]: %timeit -o rec_fortran(Tm,tau,alen) 100 loops, best of 3: 6.78 ms per loop Calling a function compiled as C extension of Python using Numpy arrays: In [65]: %timeit -o rec_c(Tm,tau,alen) 100 loops, best of 3: 6.22 ms per loop Doing everything in C using dynamic C arrays: 1000 loops,best of 3: 2.751 ms per loop 

Python-Code:

 import timeit import numpy as np from rec_cy_loop import rec_cy_loop #python setup_rec_cy_loop.py build_ext --inplace from rec_cy_loop_memviews import rec_cy_loop_memviews from rec_cy_loop_memviews_w_input import rec_cy_loop_memviews_w_input from rec_c import rec_c #python setup.py build_ext --inplace from rec_fortran import rec_fortran #f2py -c rec.f95 -m rec_fortran alen=100000 Tm = np.random.normal(size=alen).astype('float64') tau = np.random.normal(size=alen).astype('float64') def rec_numpy_loop(Tm,tau,alen): T=np.empty(alen) T[0]=0.0 for i in range(1,alen): T[i] = Tm[i] + (abs(T[i-1] - Tm[i]))**(-tau[i]) return T def rec_numpy_loop_item(Tm,tau,alen): T=np.empty(alen) Ti=T.item Tis=T.itemset Tmi=Tm.item taui=tau.item Tis(0,0.0) for i in range(1,alen): Tis(i,Tmi(i) + (abs(Ti(i-1) - Tmi(i)))**(-taui(i))) return T def it(Tm,tau): T=0.0 i=0 while True: yield T i+=1 T=Tm[i] + (abs(T - Tm[i]))**(-tau[i]) def rec_numpy_iter(Tm,tau,alen): return np.fromiter(it(Tm,tau), np.float64, alen) def rec_py_loop(Tm,tau,alen): T = [0.0] for i in range(1,alen): T.append(Tm[i] + (abs(T[i-1] - Tm[i]))**(-tau[i])) return np.array(T) T=rec_numpy_loop(Tm,tau,alen) Titer=rec_numpy_loop_item(Tm,tau,alen) np.sum(np.abs(Titer-T)) Titer=rec_numpy_iter(Tm,tau,alen) np.sum(np.abs(Titer-T)) Titer=rec_py_loop(Tm,tau,alen) np.sum(np.abs(Titer-T)) Titer=rec_cy_loop(Tm,tau,alen) np.sum(np.abs(Titer-T)) Titer=rec_cy_loop_memviews(Tm,tau,alen) np.sum(np.abs(Titer-T)) Titer=rec_cy_loop_memviews_w_input(Tm,tau,alen) np.sum(np.abs(Titer-T)) Titer=rec_fortran(Tm,tau,alen) np.sum(np.abs(Titer-T)) Titer=rec_c(Tm,tau,alen) np.sum(np.abs(Titer-T)) %timeit -o rec_numpy_loop(Tm,tau,alen) %timeit -o rec_numpy_loop_item(Tm,tau,alen) %timeit -o rec_numpy_iter(Tm,tau,alen) %timeit -o rec_py_loop(Tm,tau,alen) %timeit -o rec_cy_loop(Tm,tau,alen) %timeit -o rec_cy_loop_memviews(Tm,tau,alen) %timeit -o rec_cy_loop_memviews_w_input(Tm,tau,alen) %timeit -o rec_fortran(Tm,tau,alen) %timeit -o rec_c(Tm,tau,alen) 

Cython rec_cy_loop:

 cdef extern from "math.h": double exp(double m) import cython import numpy as np cimport numpy as np from numpy cimport ndarray @cython.boundscheck(False) @cython.wraparound(False) @cython.infer_types(True) def rec_cy_loop(ndarray[np.float64_t, ndim=1] Tm,ndarray[np.float64_t, ndim=1] tau,int alen): cdef ndarray[np.float64_t, ndim=1] T=np.empty(alen) cdef int i T[0]=0.0 for i in range(1,alen): T[i] = Tm[i] + (abs(T[i-1] - Tm[i]))**(-tau[i]) return T 

Cython rec_cy_loop_memviews:

 cdef extern from "math.h": double exp(double m) import cython import numpy as np cimport numpy as np from numpy cimport ndarray @cython.boundscheck(False) @cython.wraparound(False) @cython.infer_types(True) def rec_cy_loop_memviews(ndarray[np.float64_t, ndim=1] Tm,ndarray[np.float64_t, ndim=1] tau,int alen): cdef ndarray[np.float64_t, ndim=1] T=np.empty(alen) cdef int i cdef double[::1] T2 = T cdef double[::1] Tm2 = Tm cdef double[::1] tau2 = tau T2[0]=0.0 for i in range(1,alen): T2[i] = Tm2[i] + (abs(T2[i-1] - Tm2[i]))**(-tau2[i]) return T2 

Cython rec_cy_loop_memviews_w_input:

 cdef extern from "math.h": double exp(double m) import cython import numpy as np cimport numpy as np from numpy cimport ndarray @cython.boundscheck(False) @cython.wraparound(False) @cython.infer_types(True) def rec_cy_loop_memviews_w_input(double[::1] Tm,double[::1] tau,int alen): cdef double[::1] T=np.empty(alen) cdef int i T[0]=0.0 for i in range(1,alen): T[i] = Tm[i] + (abs(T[i-1] - Tm[i]))**(-tau[i]) return T 

Fortran-Funktion kompiliert über f2py:

 subroutine rec_fortran(tm,tau,alen,result) integer*8, intent(in) :: alen real*8, dimension(alen), intent(in) :: tm real*8, dimension(alen), intent(in) :: tau real*8, dimension(alen) :: res real*8, dimension(alen), intent(out) :: result res(1)=0 do i=2,alen res(i) = tm(i) + (abs(res(i-1) - tm(i)))**(-tau(i)) end do result=res end subroutine rec_fortran 

Reine C:

 #include <stdio.h> #include <math.h> #include <stdlib.h> #include <windows.h> #include <sys\timeb.h> double randn() { double u = rand(); if (u > 0.5) { return sqrt(-1.57079632679*log(1.0 - pow(2.0 * u - 1, 2))); } else { return -sqrt(-1.57079632679*log(1.0 - pow(1 - 2.0 * u,2))); } } void rec_pure_c(double *Tm, double *tau, int alen, double *T) { for (int i = 1; i < alen; i++) { T[i] = Tm[i] + pow(fabs(T[i - 1] - Tm[i]), (-tau[i])); } } int main() { int N = 100000; double *Tm= calloc(N, sizeof *Tm); double *tau = calloc(N, sizeof *tau); double *T = calloc(N, sizeof *T); double time = 0; double sumtime = 0; for (int i = 0; i < N; i++) { Tm[i] = randn(); tau[i] = randn(); } LARGE_INTEGER StartingTime, EndingTime, ElapsedMicroseconds; LARGE_INTEGER Frequency; for (int j = 0; j < 1000; j++) { for (int i = 0; i < 3; i++) { QueryPerformanceFrequency(&Frequency); QueryPerformanceCounter(&StartingTime); rec_pure_c(Tm, tau, N, T); QueryPerformanceCounter(&EndingTime); ElapsedMicroseconds.QuadPart = EndingTime.QuadPart - StartingTime.QuadPart; ElapsedMicroseconds.QuadPart *= 1000000; ElapsedMicroseconds.QuadPart /= Frequency.QuadPart; if (i == 0) time = (double)ElapsedMicroseconds.QuadPart / 1000; else { if (time > (double)ElapsedMicroseconds.QuadPart / 1000) time = (double)ElapsedMicroseconds.QuadPart / 1000; } } sumtime += time; } printf("1000 loops,best of 3: %.3f ms per loop\n",sumtime/1000); free(Tm); free(tau); free(T); } 

Um auf die Antwort von NPE aufzubauen, stimme ich zu, dass es irgendwo eine Schleife geben muss. Vielleicht ist dein Ziel, den Overhead zu vermeiden, der mit einem Python für Loop verbunden ist? In diesem Fall schlägt numpy.fromiter eine for-Schleife aus, aber nur ein wenig:

Mit der sehr einfachen Rekursionsrelation,

 x[i+1] = x[i] + 0.1 

Ich bekomme

 #FOR LOOP def loopit(n): x = [0.0] for i in range(n-1): x.append(x[-1] + 0.1) return np.array(x) #FROMITER #define an iterator (a better way probably exists -- I'm a novice) def it(): x = 0.0 while True: yield x x += 0.1 #use the iterator with np.fromiter def fi_it(n): return np.fromiter(it(), np.float, n) %timeit -n 100 loopit(100000) #100 loops, best of 3: 31.7 ms per loop %timeit -n 100 fi_it(100000) #100 loops, best of 3: 18.6 ms per loop 

Interessanterweise führt die Vorverteilung eines knöchernen Arrays zu einem erheblichen Leistungsverlust. Dies ist ein Rätsel für mich, obwohl ich vermute, dass es mehr Overhead mit dem Zugriff auf ein Array-Element als mit Anhängen zu einer Liste verbunden sein muss.

 def loopit(n): x = np.zeros(n) for i in range(n-1): x[i+1] = x[i] + 0.1 return x %timeit -n 100 loopit(100000) #100 loops, best of 3: 50.1 ms per loop 

Ich stolperte über diese alte Frage und sende meine Erkenntnisse und eine vektorisierte Lösung auf die Frage!

Die akzeptierte Antwort würde ich so umsetzen:

 import numpy as np np.random.seed(0) n = 100000 Tm = np.random.uniform(1, 10, size=n).astype('f') tau = np.random.uniform(-1, 0, size=n).astype('f') def calc_py(Tm_, tau_): tt = np.empty(len(Tm_)) tt[0] = Tm_[0] for i in range(1, len(Tm_)): tt[i] = (Tm_[i] - (tt[i-1] + Tm_[i])**(-tau_[i])) return tt[1:] 

Und eine vektorisierte Lösung:

 def calc_vect(Tm_, tau_): return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:]) 

Zur Überprüfung:

 print(calc_py(Tm, tau)-calc_vect(Tm, tau)) 

Gibt die Ausgabe (vermutlich schwimmende Präzisionsfehler?):

 [ -2.39640069e-06 0.00000000e+00 -1.22639676e-11 -3.82019749e-09 -3.43394937e-06 0.00000000e+00 -1.64249059e-10 -1.27897692e-13 -6.96935984e-07 -1.13686838e-13 -7.81597009e-14 -1.56319402e-13 0.00000000e+00 -1.06581410e-14 7.70565954e-07 -3.68179363e-07 -1.42108547e-14 -2.67318768e-06 0.00000000e+00 0.00000000e+00 -2.74236818e-06 4.36147587e-07 -2.05780665e-07 -5.14934904e-08] 

Jedoch kann man Numba benutzen, um die Python-Version zu kompilieren, die die gleiche Leistung wie die Vektorisierung erhält:

 from numba import jit, float32 @jit(float32[:](float32[:], float32[:]), nopython=False, nogil=True) def calc(Tm_, tau_): tt = np.empty(len(Tm_)) tt[0] = Tm_[0] for i in range(1,len(Tm_)): tt[i] = (Tm_[i] - (tt[i-1] + Tm_[i]) ** (-tau_[i])) return tt[1:] 

Timeit Ergebnisse aus den drei Lösungen:

 Python: 118.34983052355095 Vectorized: 0.9753564222721991 Numba: 0.6740745080564352 
  • Erstellen Sie Numpy Array ohne Aufzählung Array
  • Wie kann ich einen Vektor auf eine bestimmte Länge mit numpy auflösen und / oder abschneiden?
  • Erzeugen zufälliger korrelierter x- und y-Punkte mit Numpy
  • `Numpy.diff` und` scipy.fftpack.diff` geben unterschiedliche Ergebnisse bei der Unterscheidung
  • Python Pip nicht in der Lage, vcvarsall.bat zu finden
  • Systemfehler beim Ausführen von Unterprozessen mit Multiprocessing
  • Kann PySpark mit numpy Arrays arbeiten?
  • Plotten einer 2D-Wärmekarte mit Matplotlib
  • Matplotlib unerwartete Ergebnisse Polarplot
  • Migration zu numpy api 1.7
  • Wie zu optimieren ändern den Wert von 3d numpy.array, wenn eine Bedingung erfüllen
  • Python ist die beste Programmiersprache der Welt.