T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i))
Tm y tau son vectores numpy de la misma longitud que se han calculado previamente y se desea crear un nuevo vector T. La "i" se incluye solo para indicar los índices del elemento para lo que se desea.¿Es necesario un ciclo "for" si los elementos del vector numpy dependen del elemento anterior?
¿Es necesario un lazo para este caso?
Se podría pensar que esto funcionaría:
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:])
pero no es así: no se puede hacer realidad la recursividad en numpy esta manera (ya numpy calcula todo el lado derecho y luego lo asigna a la LHS) .
Así que a menos que uno se puede topar con una versión no recursiva de esta fórmula, que está pegado con un bucle explícito:
tt = np.empty(n)
tt[0] = 0.
for i in range(1,n):
tt[i] = Tm[i] + (tt[i-1] - Tm[i])**(-tau[i])
Está claro que tiene que haber un bucle en algún lugar. Supongo que tu pregunta es cómo hacer que el ciclo se produzca dentro de numpy en lugar de dentro de Python. Si esa es la verdadera pregunta, ¿qué tal un uso creativo de numpy.fromiter?
En algunos casos es posible tener este tipo de recursividad - es decir, si hay un NumPy ufunc para la fórmula de recursión, por ejemplo
c = numpy.arange(10.)
numpy.add(c[:-1], c[1:], c[1:])
Este calcula las sumas acumulativas más de c en su lugar, utilizando el parámetro de salida de numpy.add. No se puede escribir como
c[1:] = c[:-1] + c[1:]
porque ahora el resultado de la suma es un temporal que se copia en c[1:] después del cálculo está terminado.
Lo más natural intentar ahora es definir su propio ufunc:
def f(T, Tm, tau):
return Tm + (T - Tm)**(-tau)
uf = numpy.frompyfunc(f, 3, 1)
Pero por razones que están más allá de mí, lo siguiente no funciona:
uf(T[:-1], Tm[1:], tau[1:], T[1:])
Al parecer, el resultado es no escrito directamente en T[1:], sino que se almacena temporalmente y se copia después de la finalización de la operación. Incluso si funcionó, no esperaría ninguna aceleración de esto en comparación con un bucle ordinario, ya que necesita llamar a una función de Python para cada entrada.
Si realmente desea evitar el bucle de Python, probablemente necesite ir a Cython o ctypes.
Para construir sobre la respuesta de NPE, estoy de acuerdo en que tiene que haber un bucle en alguna parte. Tal vez su objetivo es evitar la sobrecarga asociada con Python for loop? En ese caso, numpy.fromiter no superó a un bucle, pero sólo por un poco:
Usando la sencilla relación de recurrencia,
x[i+1] = x[i] + 0.1
consigo
#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
Curiosamente, antes de la asignación de un numpy resultados de la matriz en una pérdida sustancial en el rendimiento. Esto es un misterio para mí, aunque supongo que debe haber más sobrecarga asociada con el acceso a un elemento de matriz que con la adición a una lista.
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
he decidido crear algunos puntos de referencia del problema original porque yo también necesitaba para calcular de forma recursiva algo que no podía ser vectorizado. Utilicé abs() para la base del exponente porque el resultado no está definido cuando la base es negativa. Ejemplo de cálculo a través de un bucle normal de la gama Numpy:
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
Los resultados son que la compilación de una función como un módulo Cython es muy rápido tanto en ejecución (14 veces más rápido que cualquier implementación de Python) y escribir el código.
Otros datos de interés: Varias implementaciones de Numpy varían ligeramente con la implementación más rápida utilizando la notación a.item() y a.itemset(). Extraño, pero agregar listas funciona a la par con las matrices numeradas preasignadas. Jugar con memviews en Cython no dio ningún impulso significativo al rendimiento del código. El código de Fortran era muy corto pero casi imposible de depurar y, al final, Cython apareció más rápido que Fortran con f2py, aunque marginalmente. La extensión C tomó mucho tiempo para escribir porque tiene demasiadas repeticiones y al final funciona con la misma velocidad que Cython. Pure C con matrices C puras fue dos veces más rápido que cualquier otro llamado desde Python.
No soy un profesional de C/C++ así que tal vez uno puede escribir un programa más rápido. código
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:
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:
rec_cy_loop_memviewscdef 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:
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
función Fortran compilado a través de 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
Pure 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);
}
Hermoso. Extremadamente útil al ver la opción de Cython también. – Chuck
me encontré con este vieja pregunta y estoy publicar mis hallazgos y una solución vectorizado a la pregunta!
La respuesta aceptada me gustaría poner en práctica la siguiente manera:
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:]
y una solución vectorizado:
def calc_vect(Tm_, tau_):
return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:])
Para la verificación:
print(calc_py(Tm, tau)-calc_vect(Tm, tau))
proporciona la salida (Presumiblemente flotante errores precicion ?):
[ -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]
Sin embargo se puede utilizar Numba a "compilar" la versión pitón obtener el mismo rendimiento que la vectorización:
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 Los resultados de las tres soluciones:
Python: 118.34983052355095
Vectorized: 0.9753564222721991
Numba: 0.6740745080564352
suena como un gran candidato para una lista de comprensión, pero no puedo intentar escribir una ahora. Me interesaría ver lo que otros piensan. – duffymo
Si 'tau' es un vector, ¿debería ser indexado por' i' también? – mtrw
¿Cuál es la condición de frontera? I.E., ¿qué sucede cuando 'i = 0'? –