如何提高这段代码的速度？（用scipy.integrate.odeint求解ODE）
我正在尝试用 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 秒，其中 odeint 花了 1 秒，而 __ConnectionistModel 花了 6 秒。
我在想有没有办法可以改善这个情况？我尝试使用 SymPy 来定义一个符号 ODE 系统，并将符号方程传递给 odeint，但这并没有正常工作，因为你不能真正定义一个符号数组，然后像数组一样访问它。
在最糟糕的情况下，我只能接受这个问题，或者使用 Cython 来加速解决过程，但我想确认我是否做对了，并且没有其他改进的办法。
谢谢你们的帮助。
[更新]：性能分析结果，
ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 7.915 7.915 grnTest.py:1(<module>)
1 0.000 0.000 7.554 7.554 grn.py:83(solve)
1 0.000 0.000 7.554 7.554 odepack.py:18(odeint)
1 0.027 0.027 7.554 7.554 {scipy.integrate._odepack.odeint}
1597 5.506 0.003 7.527 0.005 grn.py:55(__connectionistModel)
479100 1.434 0.000 1.434 0.000 grn.py:48(__sigmoid)
479102 0.585 0.000 0.585 0.000 {sum}
1 0.001 0.001 0.358 0.358 grn.py:4(<module>)
2 0.001 0.001 0.207 0.104 __init__.py:10(<module>)
27 0.014 0.001 0.185 0.007 __init__.py:1(<module>)
7 0.006 0.001 0.106 0.015 __init__.py:2(<module>)
[更新 2]：我将代码公开了：pyStGRN
1 个回答
要多多使用向量化，向量化，再向量化。而且要使用那些能帮助你进行向量化的数据结构。
这个函数 __connectionistModel 经常用到一种访问方式 A[i*m+j]，这其实就是在访问一个二维数组中的第 i 行和第 j 列，假设这个数组总共有 m 列。这说明用二维数组来存储数据是个不错的选择。我们可以通过使用 NumPy 的切片语法来消除函数中对 i 的循环，进行向量化，代码如下：
def __connectionistModel_vec(self, g, t):
"""
Returning the ODE system
"""
g_ia_s = np.zeros(self.nGenes * self.nCells)
g_2d = g.reshape((self.nCells, self.nGenes))
W = np.array(self.Params["W"])
mData = np.array(self.mData)
g_ia_s = np.zeros((self.nCells, self.nGenes))
for a in xrange(0, self.nGenes):
g_ia = self.Params["R"][a] *\
self.__sigmoid( (W[a*self.nGenes:(a+1)*self.nGenes]*g_2d).sum(axis=1) +\
self.Params["Wm"][a]*mData +\
self.Params["h"][a] ) -\
self.Params["l"][a] * g_2d[:,a]
g_ia[0] += self.Params["D"][a] * ( - 2*g_2d[0,a] + g_2d[1,a] )
g_ia[-1] += self.Params["D"][a] * ( g_2d[-2,a] - 2*g_2d[-1,a] )
g_ia[1:-1] += self.Params["D"][a] * ( g_2d[:-2,a] - 2*g_2d[1:-1,a] + g_2d[2:,a] )
g_ia_s[:,a] = g_ia
return g_ia_s.ravel()
据我所见，这样的写法返回的值和原来的 __connectionistModel 是一样的。而且，这个函数现在变得更简洁了。我只是优化了这个函数，让它的输入和输出和原来的保持一致，但如果想要更好的性能，你可能需要把数据组织成 NumPy 数组，而不是列表，这样就能避免每次调用时从列表转换成数组的过程。我相信还有其他一些小的性能优化可以做。
总之，原始代码给我的性能分析结果是这样的（这里可以插入一句“我的电脑比你的快”来炫耀一下）：
ncalls tottime percall cumtime percall filename:lineno(function)
1597 3.648 0.002 5.250 0.003 grn.py:52(__connectionistModel)
479100 0.965 0.000 0.965 0.000 grn.py:48(__sigmoid)
479100 0.635 0.000 0.635 0.000 {sum}
1 0.017 0.017 5.267 5.267 {scipy.integrate._odepack.odeint}
1598 0.002 0.000 0.002 0.000 {numpy.core.multiarray.zeros}
使用 __connectionistModel_vec 后，我得到了：
ncalls tottime percall cumtime percall filename:lineno(function)
1597 0.175 0.000 0.247 0.000 grn.py:79(__connectionistModel_vec)
4791 0.031 0.000 0.031 0.000 grn.py:48(__sigmoid)
4800 0.021 0.000 0.021 0.000 {method 'reduce' of 'numpy.ufunc' objects}
1 0.018 0.018 0.265 0.265 {scipy.integrate._odepack.odeint}
3197 0.013 0.000 0.013 0.000 {numpy.core.multiarray.array}