使用嵌套循环时加速 Python

Question

这段代码包含 3 个嵌套循环，在 Matlab 中运行需要 30 秒，而在 Python 中运行则超过 5 分钟。非常感谢任何有关如何使 Python 更快的建议！该代码通过使用 Newmark 数值解来创建响应谱。 我知道矢量化将解决挑战，但我找不到实施它的空间。

from numba import jit
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import timeit

start = timeit.default_timer()
sig=np.loadtxt('LXR1N.txt')
t=sig[:,0]
ag=0.01*sig[:,1]
h=t[1]-t[0]
f=-ag
alfa=1/6; delta=0.5
alfa0=1/(alfa*h**2)
alfa1=delta/(alfa*h)
alfa2=1/(alfa*h)
alfa3=1/(2*alfa)-1
alfa4=delta/alfa-1
alfa5= h/2*(delta/alfa-2)
alfa6= h*(1-delta)
alfa7= delta*h;
D=np.array([0.01,0.04,0.06,0.1])
T=np.zeros(1198,float)
d=np.zeros([len(t)])
len(d)
v=np.zeros([len(t)])
a=np.zeros([len(t)])
fb=np.zeros([len(t)-1])
sd=np.empty((len(D),len(T)))
sa=np.empty((len(D),len(T)))
sv=np.empty((len(D),len(T)))
aex=np.empty((len(D),len(T)))
vex=np.empty((len(D),len(T)))
for i in range(0,len(T)):
        T[i]=0.015+0.005*i
i=None
s=-1
for q in D:
    s+=1
    for i in range(0,len(T)):
        omega=2*np.pi/T[i]
        c=2*q*omega
        k=omega**2
        kb=k+alfa0+alfa1*c
        d[0]=0
        v[0]=0
        a[0]=0
        for j in range(0,len(t)-1):
            fb[j]=f[j]+alfa0*d[j]+alfa2*v[j]+alfa3*a[j]+c*alfa1*d[j]+c*alfa4*v[j]+c*alfa5*a[j] 
            d[j+1]=fb[j]/kb
            a[j+1]=alfa0*d[j+1]-alfa0*d[j]-alfa2*v[j]-alfa3*a[j]
            v[j+1]=v[j]+alfa6*a[j]+alfa7*a[j+1]
        atot=a+ag
        atot1=-omega**2*d-2*q*omega*v
        sd[s,i]=max(abs(d))   
        sa[s,i]=omega**2*sd[s,i]
        sv[s,i]=omega*sd[s,i]
        aex[s,i]=max(abs(atot))
        vex[s,i]=max(abs(v))
del i
del j
plt.figure(1,figsize=(10, 4))
for i in range(sd.shape[0]):
    plt.plot(T,sd[i,:],label=str(D[i]*100)+"%")
    plt.xlabel("Period [s]")
    plt.ylabel("Sd [m]")
plt.xlim([0,6])
plt.ylim([0,0.6])
plt.legend()
plt.savefig('Sd.jpg')
plt.show()
plt.figure(2,figsize=(6, 3))
for j in range(sa.shape[0]):
    plt.plot(T,1/9.81*sa[j,:],label=str(D[j]*100)+"%")
    plt.xlabel("Period [s]")
    plt.ylabel("Sa [g]")
plt.xlim([0,6])
plt.ylim([0,2.5])
plt.legend()
plt.savefig('Sa.png')
plt.show()
stop = timeit.default_timer()
execution_time = stop - start

print("Program Executed in: ",execution_time, 'seconds')

Answer 1

您是否尝试过分析？它应该能够告诉你时间成本在哪里。尝试使用 cProfile 模块，然后使用 snakeViz 进行可视化。分析完代码后，您可以根据其 IO 限制或 CPU 限制来优化某些部分。

1 示例可以在可能的情况下用迭代器替换 for 循环。即

for i in range(0,len(T)):
        T[i]=0.015+0.005*i

到

T = [0.015+0.005*i for i in range(len(T))]

使用 time.time() 为 len(T) = 10000 计时，在我的 PC 上，第二个代码的运行速度快了 0.0009 秒。

您可能还想看看多线程或多处理解决方案。通常这会使代码更加复杂，但允许 python 使用更多可用的处理能力。如果您受 IO 限制，请使用多线程；如果您受 CPU 限制，请使用多处理。

虽然没有分析和/或知道数据集的大小，但很难提出任何可以保证提高性能的更改。

PS。请查看有关 Python 多处理和线程的文档。

https://docs.python.org/3.8/library/multiprocessing.html

https://docs.python.org/3/library/threading.html

Answer 2

编译一下就行了

对不应在纯 Python 代码中使用的代码进行细微更改，例如 np.zeros([len(t)])，仅使用 Numpy 函数的函数部分很容易使用 Numba 进行编译。

hpaulj 已经在 cmets 中解释说，在这种情况下，Matlab 的性能优势与 jit 编译有关。尽管有一些限制（仅支持大多数 Numpy 函数），但在 Python 中也很容易做到这一点。

实施

import numba as nb
import numpy as np

@nb.njit(error_model="numpy")
def func_nb(sig):
    t=sig[:,0]
    ag=0.01*sig[:,1]
    h=t[1]-t[0]
    f=-ag
    alfa=1/6; delta=0.5
    alfa0=1/(alfa*h**2)
    alfa1=delta/(alfa*h)
    alfa2=1/(alfa*h)
    alfa3=1/(2*alfa)-1
    alfa4=delta/alfa-1
    alfa5= h/2*(delta/alfa-2)
    alfa6= h*(1-delta)
    alfa7= delta*h
    D=np.array([0.01,0.04,0.06,0.1])
    T=np.zeros(1198)
    d=np.zeros(len(t))
    v=np.zeros(len(t))
    a=np.zeros(len(t))
    fb=np.zeros((len(t)-1))
    sd=np.empty((len(D),len(T)))
    sa=np.empty((len(D),len(T)))
    sv=np.empty((len(D),len(T)))
    aex=np.empty((len(D),len(T)))
    vex=np.empty((len(D),len(T)))
    for i in range(0,len(T)):
        T[i]=0.015+0.005*i

    for s in range(D.shape[0]):
        q=D[s]
        for i in range(0,len(T)):
            omega=2*np.pi/T[i]
            c=2*q*omega
            k=omega**2
            kb=k+alfa0+alfa1*c
            d[0]=0
            v[0]=0
            a[0]=0
            for j in range(0,len(t)-1):
                fb[j]=f[j]+alfa0*d[j]+alfa2*v[j]+alfa3*a[j]+c*alfa1*d[j]+c*alfa4*v[j]+c*alfa5*a[j] 
                d[j+1]=fb[j]/kb
                a[j+1]=alfa0*d[j+1]-alfa0*d[j]-alfa2*v[j]-alfa3*a[j]
                v[j+1]=v[j]+alfa6*a[j]+alfa7*a[j+1]
            atot=a+ag
            atot1=-omega**2*d-2*q*omega*v
            sd[s,i]=np.max(np.abs(d))   
            sa[s,i]=omega**2*sd[s,i]
            sv[s,i]=omega*sd[s,i]
            aex[s,i]=np.max(np.abs(atot))
            vex[s,i]=np.max(np.abs(v))
    return sd,sa,sv,aex,vex

时间

sig=np.loadtxt("LXR1N.txt")
%timeit sd,sa,sv,aex,vex=func(sig)
#1.33 s ± 4.89 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

#I simply don't want to wait for multiple runs
import time
t1=time.time()
sd,sa,sv,aex,vex=orig()
print(time.time()-t1)
#229.00757384300232