0

我正在尝试使用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 秒odeint和 6秒__ConnectionistModel。

我想知道是否有办法可以改善这一点？我尝试使用 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

4

1 回答 1

4

向量化，向量化，然后再向量化一些。并使用有助于矢量化的数据结构。

该函数__connectionistModel使用了大量的访问模式A[i*m+j]，相当于访问一共有列的二维数组中的行i和列。这表明二维数组是存储数据的正确方法。我们可以通过使用 NumPy 切片符号和向量化来消除函数中的循环，如下所示：jmi

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}
于 2014-06-24T21:26:43.580 回答