tnt*_*tnt 30 python numpy vectorization
T(i) = Tm(i) + (T(i-1)-Tm(i))**(-tau(i))
Tm并且tau是先前已经计算过的具有相同长度的NumPy向量,并且期望创建新向量T.在i它只是为了说明什么是想要的元素索引.
这种情况需要for循环吗？
And*_*ffe 17
您可能认为这会起作用:
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])
Sve*_*ach 10
在某些情况下,可能会有这种递归 - 即如果递归公式有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.
我执行了一些基准测试,并且在2018年使用Numba是人们应该尝试加速Numpy中的递归函数的第一个选项(Aronstef的提议).Numba已经预装在Anaconda软件包中,并且速度最快(比任何Python快20倍).在2018年,Python支持@numba注释而无需额外的步骤(至少版本3.6和3.7).以下是两个基准:一个在2018-10-20执行,另一个在2016-05-18执行.
而且,正如Jaffe所说,在2018年,仍然无法对递归函数进行矢量化.我检查了Aronstef的矢量化,它不起作用.
按执行时间排序的基准:
-----------------------------------
|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) | | |
-----------------------------------
答案的最后提供了相应的代码.
我没有在2018年检查Pure C,但我认为它仍然是基于之前基准测试的最快速度.
我也没有在2018年检查C扩展,我认为它与基于之前基准的Cython几乎相同.
Fortran非常难以调试和编译,因此我没有在2018年检查f2py版本.而且它比Cython更糟糕.
我在2018年进行了以下设置:
Processor: Intel i7-7500U 2.7GHz
Versions:
Python: 3.7.0
Numba: 0.39.0
Cython: 0.28.5
Numpy: 1.15.1
使用float32(来自Aronstef)推荐的Numba代码:
@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:]
所有其他代码:
数据创建(如Aronstef + Mike T评论):
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)
2016年的代码略有不同,因为我使用abs()函数来防止nans而不是Mike T.的变体.在2018年,该函数与OP(原始海报)写的完全相同.
使用Jupyter %% magic的Cython float32.该功能可直接使用Python.Cython需要一个编译Python的C++编译器.安装正确版本的Visual C++编译器(适用于Windows)可能会出现问题:
%%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
使用Jupyter %% magic的Cython float64.该功能可直接用于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:]
附加到列表.最快的非编译解决方案:
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)
使用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:]
绕过Numpy(根据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:]
绕过Numpy(Aronstef的代码).在我的电脑上float64是默认类型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:]
纯C完全不使用Python.2016年版(带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.功能可以使用Python.2016年版(带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
更新:21-10-2018 我根据评论更正了我的答案.
只要计算不是递归的,就可以对向量进行矢量化操作.因为递归操作取决于先前计算的值,所以不可能并行处理操作.因此,这不起作用:
def calc_vect(Tm_, tau_):
return Tm_[1:] - (Tm_[:-1] + Tm_[1:]) ** (-tau_[1:])
由于(串行处理/循环)是必要的,因此通过尽可能接近优化的机器代码来获得最佳性能,因此Numba和Cython是这里的最佳答案.
Numba方法可以实现如下:
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))
Python和Numba函数的Timeit比较:
Python Solution: 54.58057559299999
Numba Soltution: 1.1389029540000024
| 归档时间： |
|
| 查看次数： |
5532 次 |
| 最近记录： |