【问题标题】：Structuring a list of objects in python for vectorization: Can a list of structures(objects) be vectorized, or are explicit arrays required在python中构建对象列表以进行矢量化：可以对结构（对象）列表进行矢量化，还是需要显式数组
【发布时间】：2019-01-11 18:38:08
【问题描述】：

分子模拟中的能量计算本质上充满了“for”循环。传统上，每个原子/分子的坐标存储在数组中。数组很容易矢量化，但结构很适合编码。将分子视为单独的对象，每个对象都有自己的坐标和其他属性，就簿记而言非常方便和清晰。

我使用的是 Python 3.6 版

我的问题是，当我使用对象数组时，我无法弄清楚如何对计算进行矢量化......似乎无法避免 for 循环。我是否有必要使用数组来利用 numpy 并矢量化我的代码？

这是一个使用数组的 python 示例（代码的第 121 行），并显示了快速（numpy）和慢速（'normal'）python 能量计算。

https://github.com/Allen-Tildesley/examples/blob/master/python_examples/mc_lj_module.py

使用 numpy 加速方法计算速度要快得多，因为它是矢量化的。

如果我不使用数组，而是使用一组对象，每个对象都有自己的坐标，我将如何向量化能量计算？这似乎需要使用较慢的 for 循环。

这是一个简单的示例代码，其中包含一个运行缓慢的 for 循环版本，以及一个无效的向量化尝试：

import numpy as np
import time

class Mol:  
    num = 0    
    def __init__(self, r):
        Mol.num += 1
        self.r       = np.empty((3),dtype=np.float_)
        self.r[0]     = r[0]
        self.r[1]     = r[1] 
        self.r[2]     = r[2]
    """ Alot more useful things go in here in practice"""

################################################
#                                              #
#               Main Program                   #
#                                              #
################################################
L = 5.0            # Length of simulation box (arbitrary)
r_cut_box_sq = L/2 # arbitrary cutoff - required
mol_list=[]
nmol = 1000    # number of molecules
part = 1    # arbitrary molecule to interact with rest of molecules

""" make 1000 molecules (1 atom per molecule), give random coordinates """
for i in range(nmol):
    r = np.random.rand(3) * L
    mol_list.append( Mol( r ) )

energy = 0.0

start = time.time()
################################################
#                                              #
#   Slow but functioning loop                  #
#                                              #
################################################
for i in range(nmol):
    if i == part:
        continue

    rij = mol_list[part].r - mol_list[i].r
    rij = rij - np.rint(rij/L)*L                # apply periodic boundary conditions
    rij_sq = np.sum(rij**2)  # Squared separations

    in_range = rij_sq < r_cut_box_sq                
    sr2      = np.where ( in_range, 1.0 / rij_sq, 0.0 )
    sr6  = sr2 ** 3
    sr12 = sr6 ** 2
    energy  += sr12 - sr6                    

end = time.time()
print('slow: ', end-start)
print('energy: ', energy)

start = time.time()
################################################
#                                              #
#   Failed vectorization attempt               #
#                                              #
################################################


    """ The next line is my problem, how do I vectorize this so I can avoid the for loop all together?
Leads to error AttributeError: 'list' object has no attribute 'r' """

""" I also must add in that part cannot interact with itself in mol_list"""
rij = mol_list[part].r - mol_list[:].r
rij = rij - np.rint(rij/L)*L                # apply periodic boundary conditions
rij_sq = np.sum(rij**2) 

in_range = rij_sq < r_cut_box_sq
sr2      = np.where ( in_range, 1.0 / rij_sq, 0.0 )
sr6  = sr2 ** 3
sr12 = sr6 ** 2
energy  = sr12 - sr6                    

energy = sum(energy)
end = time.time()
print('faster??: ', end-start)
print('energy: ', energy)

最后

如果在能量计算中，任何可能的解决方案都会受到影响，有必要遍历每个分子中的每个原子，现在每个分子的原子数超过 1 个，并且并非所有分子都具有相同数量的原子，因此有用于分子-分子相互作用的双 for 循环，而不是当前使用的简单的对-对相互作用。

【问题讨论】：

您希望矢量化哪些具体计算
最后一个for循环中的能量计算。能量是所有对相互作用的总和......我想同时计算所有对相互作用，然后在最后总结它们，而不是循环遍历所有对相互作用
您能否让您的示例仅包含您寻求矢量化的部分...这对未来的用户根本没有用。您的问题很可能是伪装成新问题的重复问题。
第一个 for 循环工作并为任何工作的矢量化提供了一个基线进行比较。没有一段代码可以测试是没有意义的。一切都标明了。
numpy 中最快的计算是那些使用编译代码的计算，使用数字 dtypes - 元素乘法和矩阵乘法，各种 ufunc，以及避免复制的索引。 numba 和 cython 可用于编译更多的迭代计算。但作为一般规则，这些工具不能很好地处理面向对象的代码。对象数组需要与对象列表相同的迭代（并且速度稍慢）。而cython 无法将计算简化为纯 C。

标签： python arrays numpy vectorization

【解决方案1】：

使用 itertools 库可能是前进的方向。假设你将一对分子的能量计算封装在一个函数中：

def calc_pairwise_energy((mol_a,mol_b)):
    # function takes a 2 item tuple of molecules
    # energy calculating code here
    return pairwise_energy

然后您可以使用 itertools.combinations 来获取所有分子对和 python 的内置列表推导（下面最后一行 [ ] 内的代码）：

from itertools import combinations
pairs = combinations(mol_list,2)
energy = sum( [calc_pairwise_energy(pair) for pair in pairs] )

当我意识到我没有正确回答您的问题时，我回到了这个答案。根据我已经发布的内容，成对能量计算函数看起来像这样（我对您的代码进行了一些优化）：

def calc_pairwise_energy(molecules):
    rij = molecules[0].r - molecules[1].r
    rij = rij - np.rint(rij/L)*L
    rij_sq = np.sum(rij**2)  # Squared separations
    if rij_sq < r_cut_box_sq:
        return (rij_sq ** -6) - (rij_sq ** - 3)
    else:
        return 0.0

而在单个调用中执行所有成对计算的矢量化实现可能如下所示：

def calc_all_energies(molecules):
    energy = 0
    for i in range(len(molecules)-1):
        mol_a = molecules[i]
        other_mols = molecules[i+1:]
        coords = np.array([mol.r for mol in other_mols])
        rijs = coords - mol_a.r
        # np.apply_along_axis replaced as per @hpaulj's  comment (see below)
        #rijs = np.apply_along_axis(lambda x: x - np.rint(x/L)*L,0,rijs)
        rijs = rijs - np.rint(rijs/L)*L
        rijs_sq = np.sum(rijs**2,axis=1)
        rijs_in_range= rijs_sq[rijs_sq < r_cut_box_sq]
        energy += sum(rijs_in_range ** -6 - rijs_in_range ** -3)
    return energy

这要快得多，但这里还有很多需要优化的地方。

【讨论】：

@T Burgis 我一直认为第二种方法，矢量化方法 (calc_all_energies) 快 10 倍多一点！如果您玩得更久，发现更多优化，请告诉我。
apply_along_axis 只是迭代输入数组的其他维度。这是一个便利功能，但不是真正的“矢量化”。
@hpaulj - 我想你可以分两步矢量化它，它可能会更快。先生成np.rint(x/L)*L的向量，然后从rijs向量中减去？
@hpaulj - 我看了看，apply_along_axis 可以替换为 rijs = rijs - np.rint(rijs/L)*L

【解决方案2】：

如果您想以坐标作为输入来计算能量，我假设您正在寻找成对距离。为此，您应该查看 SciPy 库。具体来说，我会看scipy.spatial.distance.pdist。可以在here找到文档。

【讨论】：