python - Np.cross 产生错误的结果，寻找一个可行的替代方案

Question

我正在重写分子动力学时间序列的分析代码。由于必须分析大量的时间步长(每次模拟运行 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。

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

Answer 1

np.cross有axisa、axisb 和axisc 关键字参数来选择输入和输出的位置arguments 是要交叉相乘的向量。我想你想使用:

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]