时间平均msd的优化答案

【问题标题】：Optimization of time average msd时间平均msd的优化
【发布时间】：2021-10-27 08:28:25
【问题描述】：

我尝试模拟 nT = 100 条轨道，每条 N = 10*10**6 步 dt = 0.02，然后计算以下定义的时间平均 MSD：

def calc_msd_1D(x, nLags):
  
    
    N = len(x)
    inv_sq_np = 1./np.sqrt(N)
    msd = np.zeros(nLags)

    
    for delta in range(0, nLags):
        r = 0;
    
        #msd_array = np.zeros(N)
        for i in range(N-(delta)):
            r += (x[i+delta] - x[i])**2
        msd[delta] = 1/(N-delta) * r
    
    
   # msd[0] -= 2*np.random.normal(0,1)**2   
    #msd[1:] += 2*np.random.normal(0,1)**2   
        


    return msd

因此使用类类型的结构为每个轨迹计算 MSD：

class Trajectory_Analysis_MSD: 
    
    def __init__(self,X, Y, nP, dT):
       
        
        # save parameters
        self.dT = dT
        self.X = X
        self.Y = Y
        self.nP = nP
        
  

    
    def getMSD(self,nLags): 
    
        # initialize memory
        self.MSD_x = np.zeros(nLags)
        self.MSD_y = np.zeros(nLags)
        
            
        # calculate the correlations for components
        self.MSD_x = calc_msd_1D(self.X, nLags)
        self.MSD_y= calc_msd_1D(self.Y, nLags)
            
        # calculate the msd
        self.msd= (self.MSD_x + self.MSD_y)

不幸的是，计算时间非常昂贵，目前我必须对多达 50000 个点的轨迹进行采样以存储 msd（平均时间约为 26 分钟）。有没有办法可以计算每个轨迹的整个轨道的时间平均 msd？可能没有保存每个数据点？

【问题讨论】：

你能举例说明这些类的用法吗？此外，getMSD() 中的“初始化内存”代码毫无用处

标签： python statistics implementation

【解决方案1】：

这个for循环：

r = 0
for i in range(N-delta):
    r += (x[i+delta] - x[i])**2

对于大 N 将非常慢，因为它在纯 Python 中迭代。我猜这（或您在其他地方拥有的类似代码）是您的瓶颈。

尝试向量化你的代码，这样所有的内部循环都在 numpy 中运行，而不是 Python：

r = np.sum((x[delta:N] - x[0:N-delta])**2)

您甚至不需要N 变量：

msd[delta] = np.mean((x[delta:] - x[:-delta])**2)

在这种情况下，也许你可以使用像np.correlate 这样的现成函数。

【讨论】：