Ускорение 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

Просто скомпилируйте его

С небольшими изменениями кода, которые также не должны использоваться в чистом коде Python, например np.zeros([len(t)]) ту часть вашей функции, которая использует только функции Numpy, легко компилировать с помощью Numba.

hpaulj уже объяснил в комментариях, что выигрыш в производительности от 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

Answer 2

Вы пробовали профилирование? Он должен быть в состоянии сказать вам, где стоимость времени. Попробуйте использовать модуль cProfile, а затем визуализируйте с помощью snakeViz. После того, как вы профилировали свой код, вы можете оптимизировать некоторые разделы в зависимости от того, связан ли он с IO или с процессором. Пример

1 может заменить циклы на итераторы, где это возможно. Т.е.

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, на моем ПК второй код выполняется на 0,0009 секунд быстрее.

Возможно, вы захотите взглянуть на многопоточные или многопроцессорные решения. Как правило, это сделает код намного более сложным, но позволит Python использовать больше вашей вычислительной мощности. Используйте многопоточность, если вы связаны с IO, и многопроцессорность, если вы связаны с процессором.

Хотя без профилирования и / или знания размера набора данных довольно сложно предложить какие-либо изменения, которые гарантировали бы повышение производительности.

PS. Изучите документы по многопроцессорности и многопоточности Python.

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

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