3

我正在重写分子动力学时间序列的分析代码。由于必须分析大量时间步长（每次模拟运行 150 000 个），因此我的代码尽可能快是非常重要的。

旧代码非常慢（实际上它需要的时间是我的 300 到 500 倍），因为它是为分析几千个 PDB 文件而不是一堆不同的模拟（大约 60 个）而编写的，每个都有150 000 个时间步长。我知道在这种情况下 C 或 Fortran 将是瑞士军刀，但我对 c 的经验是.....

因此，我试图为我的 python 代码使用尽可能多的 numpy/scipy 例程。因为我有使用 mkl 加速分发 anaconda 的许可证，所以这是一个非常显着的加速。

现在我面临一个问题，我希望我能以您理解我的意思的方式解释它。

我有三个数组，每个数组的形状为 (n, 3, 20)。第一行是我的肽的所有残差，通常在 23 到 31 左右。第二行是 xyz 顺序的坐标，第三行是一些特定的时间步长。

现在我正在计算每个时间步的每个残差的扭转。我对于形状为 (n,3,1) 的数组的代码：

def fast_torsion(d1, d2, d3):
    tt = dot(d1, np.cross(d2, d3))
    tb = dot(d1, d1) * dot(d2, d2)
    torsion = np.zeros([len(d1), 1])
    for i in xrange(len(d1)):
        if tb[i] != 0:
            torsion[i] = tt[i]/tb[i]
    return torsion

现在我尝试对具有扩展第三轴的数组使用相同的代码，但与使用 for 循环的原始慢速代码相比，叉积函数产生了错误的值。我用我的大数组尝试了这段代码，它比 for 循环解决方案快大约 10 到 20 倍，比旧代码快大约 200 倍。

我正在尝试的是 np.cross() 仅计算第二个（xyz）轴上的叉积并迭代其他两个轴。在第三轴较短的情况下，它可以正常工作，但对于大阵列，它仅适用于第一个时间步。我也尝试了轴设置，但我没有机会。

如果这是我的问题的唯一解决方案，我也可以使用 Cython 或 numba。

PS对不起我的英语我希望你能理解一切。

4

1 回答 1

3

np.crosshasaxisa和关键字参数axisb，axisc用于选择输入和输出参数中要交叉相乘的向量的位置。我想你想使用：

np.cross(d2, d3, axisa=1, axisb=1, axisc=1)

如果不包含axisc=1，则乘法的结果将位于输出数组的末尾。

此外，您可以torsion通过执行以下操作来避免显式循环数组：

torsion = np.zeros((len(d1), 1)
idx = (tb !=0)
torsion[idx] = tt[idx] / tb[idx]
于 2013-10-26T16:58:22.740 回答