Est-il possible de vectoriser calcul récursif d'un tableau NumPy où chaque élément dépend de la précédente?
https://stackoverflow.com/questions/4407984
-
08-10-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianQuestion
T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i))
Tm et tau sont des vecteurs numpy de la même longueur qui ont été précédemment calculé, et le désir est de créer un nouveau vecteur T. Le i est inclus uniquement pour indiquer l'indice d'élément pour ce qui est souhaité.
est une boucle nécessaire pour ce cas?
La solution
Vous pourriez penser cela fonctionnerait:
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:])
mais il ne fait pas:. Vous ne pouvez pas faire réellement une récursion en numpy cette façon (depuis numpy l'ensemble calcule RHS et attribue ensuite à la LHS)
Donc, sauf si vous pouvez trouver une version non récurrente de cette formule, vous êtes coincé avec une boucle explicite:
tt = np.empty(n)
tt[0] = 0.
for i in range(1,n):
tt[i] = Tm[i] + (tt[i-1] - Tm[i])**(-tau[i])
Autres conseils
Dans certains cas, il est possible d'avoir ce genre de récursivité -. À savoir s'il y a un ufunc NumPy pour la formule de récursion, par exemple
c = numpy.arange(10.)
numpy.add(c[:-1], c[1:], c[1:])
calcule les sommes cumulées sur c en place, en utilisant le paramètre de sortie de numpy.add. Il ne peut pas être écrit comme
c[1:] = c[:-1] + c[1:]
parce que maintenant le résultat de l'addition est un temporaire qui est copié dans c[1:] après le calcul est terminé.
La chose la plus naturelle est d'essayer maintenant de définir votre propre ufunc:
def f(T, Tm, tau):
return Tm + (T - Tm)**(-tau)
uf = numpy.frompyfunc(f, 3, 1)
Mais pour des raisons qui sont au-delà de moi, ce qui suit ne fonctionne pas:
uf(T[:-1], Tm[1:], tau[1:], T[1:])
Il semble que le résultat ne soit pas directement écrit à T[1:], mais stocké dans un temporaire et copié après la fin de l'opération. Même si cela a fonctionné, je ne s'attendre à aucune accélération de ce rapport à une boucle ordinaire, car il a besoin d'appeler une fonction Python pour chaque entrée.
Si vous voulez vraiment éviter la boucle Python, vous avez probablement besoin d'aller Cython ou ctypes.
J'effectué quelques repères et en 2018 en utilisant Numba est le premier peuple d'options devraient essayer d'accélérer les fonctions récursives Numpy (proposition Aronstef). Numba est déjà pré-installé dans le package Anaconda et a l'un des temps les plus rapides (environ 20 fois plus rapide que tout Python). En 2018 Python prend en charge les annotations @numba sans étapes supplémentaires (au moins les versions 3.6 et 3.7). Voici deux points de référence: une réalisée sur 20/10/2018 et l'autre sur 18/05/2016
.Et, comme mentionné par Jaffe, en 2018, il est encore impossible de vectoriser fonctions récursives. J'ai vérifié la vectorisation par Aronstef et il ne fonctionne pas.
Repères triés par temps d'exécution:
-----------------------------------
|Variant |2018-10 |2016-05 |
-----------------------------------
|Pure C | na | 2.75 ms|
|C extension | na | 6.22 ms|
|Cython float32 | 1.01 ms| na |
|Cython float64 | 1.05 ms| 6.26 ms|
|Fortran f2py | na | 6.78 ms|
|Numba float32 | 2.81 ms| na |
|(Aronstef) | | |
|Numba float64 | 5.28 ms| na |
|Append to list |48.2 ms|91.0 ms|
|Using a.item() |58.3 ms|74.4 ms|
|np.fromiter() |60.0 ms|78.1 ms|
|Loop over Numpy|71.9 ms|87.9 ms|
|(Jaffe) | | |
|Loop over Numpy|74.4 ms| na |
|(Aronstef) | | |
-----------------------------------
code correspondant est fourni à la fin de la réponse.
Je ne vérifiait pas pur C en 2018, mais je suppose qu'il est toujours le plus rapide en fonction de l'indice de référence précédent.
J'ai aussi n'a pas vérifié l'extension C en 2018 et je pense qu'il a presque en même temps que Cython sur la base de référence précédent.
Fortran a été très difficile à déboguer et de compiler, donc je n'ai pas vérifié la version f2py en 2018. Et ce fut pire que Cython de toute façon.
J'ai la configuration suivante en 2018:
Processor: Intel i7-7500U 2.7GHz
Versions:
Python: 3.7.0
Numba: 0.39.0
Cython: 0.28.5
Numpy: 1.15.1
recommandé Numba code à l'aide float32 (de Aronstef):
@numba.jit("float32[:](float32[:], float32[:])", nopython=False, nogil=True)
def calc_py_jit32(Tm_, tau_):
tt = np.empty(len(Tm_),dtype="float32")
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:]
Tous les autres codes:
Création de données (comme Aronstef + Mike T commentaire):
np.random.seed(0)
n = 100000
Tm = np.cumsum(np.random.uniform(0.1, 1, size=n).astype('float64'))
tau = np.random.uniform(-1, 0, size=n).astype('float64')
ar = np.column_stack([Tm,tau])
Tm32 = Tm.astype('float32')
tau32 = tau.astype('float32')
Tm_l = list(Tm)
tau_l = list(tau)
Le code en 2016 était légèrement différente comme je fonction abs () pour éviter nans et non la variante de Mike T. En 2018, la fonction est exactement la même que OP (Affiche originale) écrit.
Cython float32 en utilisant Jupyter %% magique. La fonction peut être utilisée directement dans Python. Cython a besoin d'un compilateur C ++ dans lequel Python a été compilé. Installation de la bonne version de compilateur Visual C ++ (pour Windows) pourrait être problématique:
%%cython
import cython
import numpy as np
cimport numpy as np
from numpy cimport ndarray
cdef extern from "math.h":
np.float32_t exp(np.float32_t m)
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.infer_types(True)
@cython.initializedcheck(False)
def cy_loop32(np.float32_t[:] Tm,np.float32_t[:] tau,int alen):
cdef np.float32_t[:] T=np.empty(alen, dtype=np.float32)
cdef int i
T[0]=0.0
for i in range(1,alen):
T[i] = Tm[i] + (T[i-1] - Tm[i])**(-tau[i])
return T
Cython float64 en utilisant Jupyter %% magique. La fonction peut être utilisée directement dans Python:
%%cython
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)
@cython.initializedcheck(False)
def cy_loop(double[:] Tm,double[:] tau,int alen):
cdef double[:] T=np.empty(alen)
cdef int i
T[0]=0.0
for i in range(1,alen):
T[i] = Tm[i] + (T[i-1] - Tm[i])**(-tau[i])
return T
Numba float64:
@numba.jit("float64[:](float64[:], float64[:])", nopython=False, nogil=True)
def calc_py_jit(Tm_, tau_):
tt = np.empty(len(Tm_),dtype="float64")
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:]
Append à la liste . Solution la plus rapide non compilé:
def rec_py_loop(Tm,tau,alen):
T = [Tm[0]]
for i in range(1,alen):
T.append(Tm[i] - (T[i-1] + Tm[i])**(-tau[i]))
return np.array(T)
Utilisation a.item ():
def rec_numpy_loop_item(Tm_,tau_):
n_ = len(Tm_)
tt=np.empty(n_)
Ti=tt.item
Tis=tt.itemset
Tmi=Tm_.item
taui=tau_.item
Tis(0,Tm_[0])
for i in range(1,n_):
Tis(i,Tmi(i) - (Ti(i-1) + Tmi(i))**(-taui(i)))
return tt[1:]
np.fromiter ():
def it(Tm,tau):
T=Tm[0]
i=0
while True:
yield T
i+=1
T=Tm[i] - (T + Tm[i])**(-tau[i])
def rec_numpy_iter(Tm,tau,alen):
return np.fromiter(it(Tm,tau), np.float64, alen)[1:]
Boucle sur Numpy (basé sur l'idée de l'Jaffe):
def rec_numpy_loop(Tm,tau,alen):
tt=np.empty(alen)
tt[0]=Tm[0]
for i in range(1,alen):
tt[i] = Tm[i] - (tt[i-1] + Tm[i])**(-tau[i])
return tt[1:]
Boucle sur Numpy (le code de Aronstef). Sur mon ordinateur float64 est le type par défaut pour np.empty.
def calc_py(Tm_, tau_):
tt = np.empty(len(Tm_),dtype="float64")
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:]
pur C sans utiliser Python du tout. Version de l'année 2016 (avec fonction de fabs ()):
#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);
}
Fortran f2py. La fonction peut être utilisée à partir Python. Version de l'année 2016 (avec fonction de abs ()):
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
Mise à jour: 21-10-2018 J'ai corrigé ma réponse en fonction des commentaires.
Il est possible de vectoriser des opérations sur des vecteurs aussi longtemps que le calcul n'est pas récursive. En raison d'une opération récursive dépend de la valeur calculée précédente, il est impossible de processus parallèle de l'opération. Cela ne permet donc pas de travail:
def calc_vect(Tm_, tau_):
return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:])
Depuis (traitement série / boucle) est nécessaire, la meilleure performance est obtenue en déplaçant le plus près possible de code machine optimisé, donc Numba et Cython sont les meilleures réponses ici.
Une approche Numba peut être permette d'obtenir comme suit:
init_string = """
from math import pow
import numpy as np
from numba import jit, float32
np.random.seed(0)
n = 100000
Tm = np.cumsum(np.random.uniform(0.1, 1, size=n).astype('float32'))
tau = np.random.uniform(-1, 0, size=n).astype('float32')
def calc_python(Tm_, tau_):
tt = np.empty(len(Tm_))
tt[0] = Tm_[0]
for i in range(1, len(Tm_)):
tt[i] = Tm_[i] - pow(tt[i-1] + Tm_[i], -tau_[i])
return tt
@jit(float32[:](float32[:], float32[:]), nopython=False, nogil=True)
def calc_numba(Tm_, tau_):
tt = np.empty(len(Tm_))
tt[0] = Tm_[0]
for i in range(1, len(Tm_)):
tt[i] = Tm_[i] - pow(tt[i-1] + Tm_[i], -tau_[i])
return tt
"""
import timeit
py_time = timeit.timeit('calc_python(Tm, tau)', init_string, number=100)
numba_time = timeit.timeit('calc_numba(Tm, tau)', init_string, number=100)
print("Python Solution: {}".format(py_time))
print("Numba Soltution: {}".format(numba_time))
comparaison TimeIt des fonctions Python et Numba:
Python Solution: 54.58057559299999
Numba Soltution: 1.1389029540000024
Pour construire sur la réponse de NPE, je suis d'accord qu'il doit y avoir une part de la boucle. Peut-être que votre objectif est d'éviter la surcharge associée à un Python pour la boucle? Dans ce cas, numpy.fromiter ne bat une boucle, mais seulement par un peu:
En utilisant la relation de récursion très simple,
x[i+1] = x[i] + 0.1
Je reçois
#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
Il est intéressant de pré-attribution d'un tableau de résultats numpy une perte substantielle de la performance. Ceci est un mystère pour moi, même si je suppose qu'il doit y avoir plus de frais généraux associés à l'accès à un élément de tableau que de annexant à une liste.
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
