T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i))
Tm и tau - числовые векторы той же длины, которые были предварительно рассчитаны, и желание создать новый вектор T. "i" включен только для указания индексов элементов для желаемого.
Является ли цикл for необходимым для этого случая?
Вы могли бы подумать, что это сработает:
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:])
но это не так: вы не можете фактически сделать рекурсию в numpy таким образом (поскольку numpy вычисляет всю RHS, а затем назначает ее LHS).
Итак, если вы не можете придумать нерекурсивную версию этой формулы, вы застряли в явном цикле:
tt = np.empty(n)
tt[0] = 0.
for i in range(1,n):
tt[i] = Tm[i] + (tt[i-1] - Tm[i])**(-tau[i])
В некоторых случаях такая рекурсия возможна, а именно, если для формулы рекурсии имеется NumPy ufunc, например
c = numpy.arange(10.)
numpy.add(c[:-1], c[1:], c[1:])
Это вычисляет накопленные суммы по c на месте, используя выходной параметр numpy.add. Он не может быть записан как
c[1:] = c[:-1] + c[1:]
потому что теперь результат добавления является временным, который копируется в c[1:] после завершения вычисления.
Самое естественное, что нужно попробовать, это определить свой собственный ufunc:
def f(T, Tm, tau):
return Tm + (T - Tm)**(-tau)
uf = numpy.frompyfunc(f, 3, 1)
Но по причинам, которые выходят за меня, следующее не работает:
uf(T[:-1], Tm[1:], tau[1:], T[1:])
По-видимому, результат не записывается непосредственно в T[1:], а сохраняется во временном и скопированном после завершения операции. Даже если бы это сработало, я бы не ожидал ускорения от этого по сравнению с обычным циклом, так как ему нужно вызвать функцию Python для каждой записи.
Если вы действительно хотите избежать цикла Python, вам, вероятно, нужно обратиться за Cython или ctypes.
Ясно, что где-то должен быть петля. Думаю, ваш вопрос заключается в том, как сделать цикл внутри numpy, а не внутри самого Python. Если это реальный вопрос, то как насчет некоторого творческого использования numpy.fromiter?
Я наткнулся на этот старый вопрос и публикую свои выводы и векторизованное решение вопроса!
Принятый ответ я бы выполнил следующим образом:
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:]
И векторизованное решение:
def calc_vect(Tm_, tau_):
return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:])
Для проверки:
print(calc_py(Tm, tau)-calc_vect(Tm, tau))
Дает выход (предположительно плавающие ошибки ошибки?):
[ -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]
Однако можно использовать Numba для "компиляции" версии python, получающей ту же производительность, что и векторизация:
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 Результаты трех решений:
Python: 118.34983052355095
Vectorized: 0.9753564222721991
Numba: 0.6740745080564352
Я решил создать несколько тестов исходной проблемы, потому что мне также нужно было вычислить что-то рекурсивно, что невозможно было векторизовать. Я использовал abs() для базы экспонентов, потому что результат undefined, когда база отрицательна. Пример вычисления путем нормального цикла массива 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
Результаты состоят в том, что компиляция функции в виде модуля Cython очень быстро выполняется как в исполнении (в 14 раз быстрее, чем любая реализация Python), так и при написании кода.
Другие интересные факты: Различные реализации Numpy немного отличаются от самой быстрой реализации с использованием нотации a.item() и a.itemset(). Странно, но добавление списков работает наравне с добавлением предварительно распределенных массивов Numpy. Воспроизведение с memviews в Cython не дало значительного повышения производительности кода. Код Fortran был очень коротким, но почти невозможно отладить, и в итоге Cython появился быстрее Fortran с f2py, хотя и незначительно. C потребовалось много времени, потому что у него слишком много шаблонов, и в конце он работает с той же скоростью, что и Cython. Чистый C с чистыми C-массивами был в два раза быстрее, чем что-либо, вызванное Python.
Я не являюсь профессионалом C/С++, поэтому, возможно, можно написать более быструю программу.
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:
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, скомпилированная через 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
Чистый 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);
}
Чтобы построить ответ NPE, я согласен, что где-то должен быть цикл. Возможно, ваша цель - избежать накладных расходов, связанных с циклом Python for? В этом случае numpy.fromiter выбивает цикл for, но только немного:
Используя очень простое рекуррентное соотношение,
x[i+1] = x[i] + 0.1
Я получаю
#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
Интересно, что предварительное выделение массива numpy приводит к существенной потере производительности. Это для меня загадка, хотя я бы предположил, что должно быть больше накладных расходов, связанных с доступом к элементу массива, чем с добавлением в список.
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