我试图用scipy.integrate.odeint
模块来解决一个相对较大的ODE系统。我已经实现了代码,我可以正确地求解方程。但这个过程非常缓慢。我分析了代码,我意识到几乎大部分的计算时间都花在计算或生成ODE系统本身上,sigmoid函数也很昂贵,但我想我必须接受它。下面是我正在使用的一段代码:
def __sigmoid(self, u):
# return .5 * ( u / np.sqrt(u**2 + 1) + 1)
return 0.5 + np.arctan(u) / np.pi
def __connectionistModel(self, g, t):
"""
Returning the ODE system
"""
g_ia_s = np.zeros(self.nGenes * self.nCells)
for i in xrange(0, self.nCells):
for a in xrange(0, self.nGenes):
g_ia = self.Params["R"][a] *\
self.__sigmoid( sum([ self.Params["W"][b + a*self.nGenes]*g[self.nGenes*i + b] for b in xrange(0, self.nGenes) ]) +\
self.Params["Wm"][a]*self.mData[i] +\
self.Params["h"][a] ) -\
self.Params["l"][a] * g[self.nGenes*i + a]
# Uncomment this part for including the diffusion
if i == 0:
g_ia += self.Params["D"][a] * ( - 2*g[self.nGenes*i + a] + g[self.nGenes*(i+1) + a] )
elif i == self.nCells-1:
g_ia += self.Params["D"][a] * ( g[self.nGenes*(i-1) + a] - 2*g[self.nGenes*i + a] )
else:
g_ia += self.Params["D"][a] * ( g[self.nGenes*(i-1) + a] - 2*g[self.nGenes*i + a] + g[self.nGenes*(i+1) + a] )
g_ia_s[self.nGenes*i + a] = g_ia
return g_ia_s
def solve(self, inp):
g0 = np.zeros(self.nGenes * self.nCells)
t = np.arange(self.t0, self.t1, self.dt)
self.integratedExpression = odeint(self.__connectionistModel, g0, t, full_output=0)
return self.integratedExpression
正如您在每次迭代中看到的,我应该生成nCells*nGenes(100*3=300)方程并将其传递给odeint
。虽然我不确定,但我想生成方程比求解它们要贵得多。在我的实验中,求解整个系统需要7秒,其中包括1秒的{
我想知道有没有办法我可以改进一下?我试图使用SymPy来定义一个符号ODE系统,并将符号方程传递给odeint
,但是它不能正常工作,因为你不能真正定义一个符号数组,以后你可以像数组一样访问这些符号。在
在最坏的情况下,我必须处理它或使用一个Cython来加快解决过程,但我想确保我做的是正确的,没有办法改善它。在
提前谢谢你的帮助。在
[Update]:分析结果
^{pr2}$[更新2]:我公开了代码:pyStGRN
矢量化,矢量化,然后再矢量化一些。并使用有助于矢量化的数据结构。在
函数
__connectionistModel
使用了大量的访问模式A[i*m+j]
,这相当于在一个总共有m
列的2D数组中访问i
和列j
。这表明二维数组是存储数据的正确方法。通过使用NumPy切片表示法和向量化,我们可以从函数中消除i
上的循环,如下所示:据我所知,它返回与原始
__connectionistModel
相同的值。另外,该函数现在更加紧凑。我只对这个函数进行了优化,使其具有与原始函数相同的输入和输出,但是为了获得额外的性能,您可能需要将数据组织在NumPy数组中,而不是列表中,以避免每次调用时从列表到数组的转换。我相信还有其他一些小的性能调整。在不管怎样,原始代码给了我这些分析结果(在这里插入强制性的“我的电脑比你的电脑快”吹嘘):
^{pr2}$使用
__connectionistModel_vec
,我得到:相关问题 更多 >
编程相关推荐