我有一个常微分方程要解决，它是用于心脏细胞建模的菲茨休南云方程。我制作了一个使用欧拉方法求解两个常微分方程的代码。所以我有这个：
import numpy as np
from numba import jitclass
from numba import int32, float64
import matplotlib.pyplot as plt
import time
spec = [('V_init' ,float64),
('a' ,float64),
('b' ,float64),
('g',float64),
('dt' ,float64),
('NbODEs',int32),
('dydx' ,float64[:]),
('y' ,float64[:]) ]
@jitclass(spec, )
class FHNfunc:
def __init__(self,):
self.V_init = .04
self.a= 0.25
self.b=0.001
self.g = 0.003
self.dt = .01
self.NbODEs = 2
self.dydx =np.zeros(self.NbODEs, )
self.y =np.zeros(self.NbODEs, )
def Eul(self):
self.deriv()
self.y += (self.dydx * self.dt)
def deriv(self , ):
self.dydx[0]= self.V_init - self.y[0] *(self.a-(self.y[0]))*(1-(self.y[0]))-self.y[1]
self.dydx[1]= self.b * self.y[0] - self.g * self.y[1]
return
FH = FHNfunc()
dt = .001
tp = np.linspace(0, 1000, num = int((1000)/dt))
V = np.zeros(len(tp), )
W = np.zeros(len(tp), )
t0 = time.time()
for idx, t in enumerate(tp):
FH.Eul()
V[idx] = FH.y[0]
W[idx] = FH.y[1]
print(time.time()- t0)
plt.subplots()
plt.plot(tp,V)
plt.plot(tp,W)
plt.show()
我尝试的是使用该numba jitclass来提高 FHN 常微分方程求解的时间性能，但我希望它没有那么有帮助。 对于这个例子，代码给了我 11.44 秒而不使用 jitclass (当我评论时@jitclass(spec, ))和 6.14 秒使用 jitclass 使用。我不是在抱怨获得两倍的计算时间，但我期望更多。我知道我可以在类中集成 for 循环，但我需要在外面。 因此，我正在寻找解决方案来进一步改善此示例的计算时间。
编辑：这次我尝试在类外使用 jit 实现 ODE 函数：
__author__ = 'Maxime'
import numpy as np
from numba import jitclass, jit
from numba import int32, float64
import matplotlib.pyplot as plt
import time
spec = [('V_init' ,float64),
('a' ,float64),
('b' ,float64),
('g',float64),
('dt' ,float64),
('NbODEs',int32),
('dydx' ,float64[:]),
('time' ,float64[:]),
('V' ,float64[:]),
('W' ,float64[:]),
('y' ,float64[:]) ]
# @jitclass(spec, )
class FHNfunc:
def __init__(self,):
self.V_init = .04
self.a= 0.25
self.b=0.001
self.g = 0.003
self.dt = .001
self.NbODEs = 2
self.dydx =np.zeros(self.NbODEs )
self.y =np.zeros(self.NbODEs )
def Eul(self):
self.deriv()
self.y += (self.dydx * self.dt)
def deriv(self):
# self.dydx[0]= self.V_init - self.y[0] *(self.a-(self.y[0]))*(1-(self.y[0]))-self.y[1]
# self.dydx[1]= self.b * self.y[0] - self.g * self.y[1]
self.dydx[0]= fV(self.V_init,self.y[0],self.y[1],self.a)
self.dydx[1]= fW(self.y[0],self.y[1],self.b,self.g)
return
@jit(float64(float64, float64, float64, float64))
def fV(V_init,y0,y1,a):
return V_init - y0 *(a-(y0))*(1-(y0))-y1
@jit(float64(float64, float64, float64, float64))
def fW(y0,y1,b,g):
return b * y0 - g * y1
FH = FHNfunc()
dt = .001
tp = np.linspace(0, 1000, num = int((1000)/dt))
V = np.zeros(len(tp), )
W = np.zeros(len(tp), )
t0 = time.time()
for idx, t in enumerate(tp):
FH.Eul()
V[idx] = FH.y[0]
W[idx] = FH.y[1]
print(time.time()- t0)
plt.subplots()
plt.plot(tp,V)
plt.plot(tp,W)
plt.show()
但在这种情况下，我根本没有时间改进：11.4 秒。
为什么我不能在类中拥有积分循环
当我有多个模型并且想要它们之间的耦合时，我需要在 FHN 实例之间传递变量。例如：
__author__ = 'Maxime'
import numpy as np
from numba import jitclass, jit, njit
from numba import int32, float64
import matplotlib.pyplot as plt
import time
spec = [('V_init' ,float64),
('a' ,float64),
('b' ,float64),
('g',float64),
('dt' ,float64),
('NbODEs',int32),
('dydx' ,float64[:]),
('time' ,float64[:]),
('V' ,float64[:]),
('W' ,float64[:]),
('y' ,float64[:]) ]
@jitclass(spec, )
class FHNfunc:
def __init__(self,):
self.V_init = .04
self.a= 0.25
self.b=0.001
self.g = 0.003
self.dt = .001
self.NbODEs = 2
self.dydx =np.zeros(self.NbODEs )
self.y =np.zeros(self.NbODEs )
def Eul(self):
self.deriv()
self.y += (self.dydx * self.dt)
def deriv(self):
self.dydx[0]= self.V_init - self.y[0] *(self.a-(self.y[0]))*(1-(self.y[0]))-self.y[1]
self.dydx[1]= self.b * self.y[0] - self.g * self.y[1]
return
FH1 = FHNfunc()
FH2 = FHNfunc()
FH2.V_init=0.
dt = .001
tp = np.linspace(0, 1000, num = int((1000)/dt))
V1 = np.zeros(len(tp), )
V2 = np.zeros(len(tp), )
W1 = np.zeros(len(tp), )
W2 = np.zeros(len(tp), )
t0 = time.time()
for idx, t in enumerate(tp):
FH1.Eul()
FH2.V_init=FH1.V_init
FH2.Eul()
V1[idx] = FH1.y[0]
W1[idx] = FH1.y[1]
V2[idx] = FH2.y[0]
W2[idx] = FH2.y[1]
print(time.time()- t0)
plt.figure
plt.subplot(211)
plt.plot(tp,V1)
plt.plot(tp,W1)
plt.subplot(212)
plt.plot(tp,V2)
plt.plot(tp,W2)
plt.show()
在这种情况下，我不知道如何将 numpy 与在实例之间传递的变量一起使用。此外，对于此示例，所有 nstance 都属于同一类，但在我的完整模型中，我有 8 个不同的类来表示属于系统的不同类型的模型。
@max9111的答案
所以我用njit测试了它，两个神经元连接在一起，效果很好：
__author__ = 'Maxime'
import numpy as np
from numba import jitclass, jit, njit
from numba import int32, float64
import matplotlib.pyplot as plt
import time
spec = [('V_init' ,float64),
('a' ,float64),
('b' ,float64),
('g',float64),
('dt' ,float64),
('NbODEs',int32),
('dydx' ,float64[:]),
('time' ,float64[:]),
('V' ,float64[:]),
('W' ,float64[:]),
('y' ,float64[:]) ]
@jitclass(spec, )
class FHNfunc:
def __init__(self,):
self.V_init = .04
self.a= 0.25
self.b=0.001
self.g = 0.003
self.dt = .001
self.NbODEs = 2
self.dydx =np.zeros(self.NbODEs )
self.y =np.zeros(self.NbODEs )
def Eul(self,):
self.deriv()
self.y += (self.dydx * self.dt)
def deriv(self,):
self.dydx[0]= self.V_init - self.y[0] *(self.a-(self.y[0]))*(1-(self.y[0]))-self.y[1]
self.dydx[1]= self.b * self.y[0] - self.g * self.y[1]
return
@njit(fastmath=True)
def solve2(FH1,FH2,tp):
V1 = np.zeros(len(tp), )
V2 = np.zeros(len(tp), )
W1 = np.zeros(len(tp), )
W2 = np.zeros(len(tp), )
for idx, t in enumerate(tp):
FH1.Eul()
FH2.V_init=FH1.V_init
FH2.Eul()
V1[idx] = FH1.y[0]
W1[idx] = FH1.y[1]
V2[idx] = FH2.y[0]
W2[idx] = FH2.y[1]
return V1,W1,V2,W2
if __name__ == "__main__":
#with njit and jiclass
FH1 = FHNfunc()
FH2 = FHNfunc()
FH2.V_init=0.
dt = .001
tp = np.linspace(0, 1000, num = int((1000)/dt))
t0 = time.time()
[V1,W1,V2,W2] = solve2(FH1,FH2,tp)
print(time.time()- t0)
plt.figure()
plt.subplot(211)
plt.plot(tp,V1)
plt.plot(tp,W1)
plt.subplot(212)
plt.plot(tp,V2)
plt.plot(tp,W2)
#with jitclass only
FH1 = FHNfunc()
FH2 = FHNfunc()
FH2.V_init=0.
dt = .001
tp = np.linspace(0, 1000, num = int((1000)/dt))
t0 = time.time()
V1 = np.zeros(len(tp), )
V2 = np.zeros(len(tp), )
W1 = np.zeros(len(tp), )
W2 = np.zeros(len(tp), )
for idx, t in enumerate(tp):
FH1.Eul()
FH2.V_init=FH1.V_init
FH2.Eul()
V1[idx] = FH1.y[0]
W1[idx] = FH1.y[1]
V2[idx] = FH2.y[0]
W2[idx] = FH2.y[1]
print(time.time()- t0)
plt.figure()
plt.subplot(211)
plt.plot(tp,V1)
plt.plot(tp,W1)
plt.subplot(212)
plt.plot(tp,V2)
plt.plot(tp,W2)
plt.show()
有了这个，我有 1.8 秒的所有优化(NJIT 和 jitclass)和模型的两个实例。我只有 12.4 秒的 jitclass 和 21.7 秒，根本没有 numba。所以是 12 倍，一点也不差。 感谢@max9111提供的解决方案。
这一切都与函数内联和 LLVM 优化有关
所有函数都是非常原始的(关于计算时间)。因此，numba 在这里唯一能做的就是内联此函数并缓存已编译的函数，以避免下一次调用时的编译开销。
您的 Jitclass 基准测试有一个主要问题。您正在从未编译的代码调用基元函数 1000000 次。(意味着至少 1000000 次函数调用)。这应该看起来像：
使用抖动类Example_1
import numpy as np
from numba import jitclass,njit
from numba import int32, float64
import matplotlib.pyplot as plt
import time
spec = [('V_init' ,float64),
('a' ,float64),
('b' ,float64),
('g',float64),
('dt' ,float64),
('NbODEs',int32),
('dydx' ,float64[:]),
('y' ,float64[:]) ]
@jitclass(spec)
class FHNfunc:
def __init__(self,):
self.V_init = .04
self.a= 0.25
self.b=0.001
self.g = 0.003
self.dt = .001
self.NbODEs = 2
self.dydx =np.zeros(self.NbODEs, )
self.y =np.zeros(self.NbODEs, )
def Eul(self):
self.deriv()
self.y += (self.dydx * self.dt)
def deriv(self , ):
self.dydx[0]= self.V_init - self.y[0] *(self.a-(self.y[0]))*(1-(self.y[0]))-self.y[1]
self.dydx[1]= self.b * self.y[0] - self.g * self.y[1]
return
@njit(fastmath=True)
def solve(FH,dt,tp):
V = np.zeros(len(tp), )
W = np.zeros(len(tp), )
for idx, t in enumerate(tp):
FH.Eul()
V[idx] = FH.y[0]
W[idx] = FH.y[1]
return V,W
if __name__ == "__main__":
FH = FHNfunc()
dt = .001
tp = np.linspace(0, 1000, num = int((1000)/dt))
t1=time.time()
[V,W]=solve(FH,dt,tp)
print(time.time()-t1)
plt.subplots()
plt.plot(tp,V)
plt.plot(tp,W)
plt.show()
这给出了大约 0.4 秒的运行时间。
Example_2和 3
import numpy as np
import numba as nb
import matplotlib.pyplot as plt
import time
@nb.njit(fastmath=True,cache=True)
def Eul(V_init,y,a,g,dt,dydx):
deriv(V_init,y,a,b,g,dydx)
y += (dydx * dt)
@nb.njit(fastmath=True,cache=True)
def deriv(V_init,y,a,b,g,dydx):
dydx[0]= fV(V_init,y[0],y[1],a)
dydx[1]= fW(y[0],y[1],b,g)
@nb.njit(fastmath=True,cache=True)
def fV(V_init,y0,y1,a):
return V_init - y0 *(a-(y0))*(1-(y0))-y1
@nb.njit(fastmath=True,cache=True)
def fW(y0,y1,b,g):
return b * y0 - g * y1
@nb.njit(fastmath=True,cache=True)
def solving_1(V_init,y,a,g,dt,tp):
V = np.empty(tp.shape[0],dtype=y.dtype)
W = np.empty(tp.shape[0],dtype=y.dtype)
dydx=np.empty(2,dtype=np.float64)
for idx, t in enumerate(tp):
Eul(V_init,y,a,g,dt,dydx)
V[idx] = y[0]
W[idx] = y[1]
return V,W
@nb.njit(fastmath=True,cache=True)
def solving_2(V_init,y,a,g,dt,tp):
V = np.empty(tp.shape[0],dtype=y.dtype)
W = np.empty(tp.shape[0],dtype=y.dtype)
dydx=np.empty(2,dtype=y.dtype)
for idx, t in enumerate(tp):
dydx[0]=V_init - y[0] *(a-(y[0]))*(1-(y[0]))-y[1]
dydx[1]=b * y[0] - g * y[1]
y[0] += (dydx[0] * dt)
y[1] += (dydx[1] * dt)
V[idx] = y[0]
W[idx] = y[1]
return V,W
if __name__ == "__main__":
V_init = .04
a= 0.25
b=0.001
g = 0.003
dt = .001
dt = .001
tp = np.linspace(0, 1000, num = int((1000)/dt))
y=np.zeros(2,dtype=np.float64)
t1=time.time()
[V,W]=solving_2(V_init,y,a,g,dt,tp)
print(time.time()-t1)
plt.subplots()
plt.plot(tp,V)
plt.plot(tp,W)
plt.show()
我在这里测试了两个变体。所有工作都在一个功能中，并拆分为服务器功能。这为solving_1提供了 0.17 秒，为 solving_2 提供了 0.06 秒。
我并不感到惊讶，jitclass approuch 有点慢(不支持缓存，相当新的功能)，但我没想到在方法solving_1和solving_2中看到性能的第二个因素，如果有人使用一些内存复制，它们甚至更大，这也没有优化掉。