根据索引和值填充 scipy / numpy 矩阵

Question

我有一张节点图，每个节点代表大脑中大约 100 个体素。我将图划分为社区，但现在我需要制作一个相关矩阵，其中节点中的每个体素都连接到同一社区中节点中的每个体素。换句话说，如果节点 1 和 2 在同一个社区中，我需要在节点 1 中的每个体素和节点 2 中的每个体素之间的矩阵中有一个 1。下面的代码需要很长时间。有谁知道如何加快速度？

for edge in combinations(graph.nodes(),2):
    if partition.get_node_community(edge[0]) == partition.get_node_community(edge[1]): # if nodes are in same community
        voxels1 = np.argwhere(flat_parcel==edge[0]+1) # this is where I find the voxels in each node, and I get the indices for the matrix where I want them.
        voxels2 = np.argwhere(flat_parcel==edge[1]+1)
        for voxel1 in voxels1:
            voxel_matrix[voxel1,voxels2] = 1

感谢您的回复，我认为最简单和最快的解决方案是将最后一个循环替换为 voxel_matrix[np.ix_(体素 1, 体素 2)] = 1

Answer 1

这是我希望对您有用的一种方法。这对我的机器来说是一种延伸——即使存储体素邻接矩阵的两个副本（使用 dtype=bool）也会将我的（有点旧的）桌面推到其内存容量的边缘。但我假设您有一台能够处理至少两个 (300 * 100) ** 2 = 900 MB 数组的机器——否则，在此阶段之前您可能会遇到运行问题。我的桌面大约需要 30 分钟才能处理 30000 个体素。

这假设 voxel_communities 是一个简单数组，其中包含索引 i 处每个体素的社区标签。听起来您可以很快生成它。它还假定体素仅存在于一个节点中。

def voxel_adjacency(voxel_communities):
    n_voxels = voxel_communities.size
    comm_labels = sorted(set(voxel_communities))
    comm_counts = [(voxel_communities == l).sum() for l in comm_labels]

    blocks = numpy.zeros((n_voxels, n_voxels), dtype=bool)
    start = 0
    for c in comm_counts:
        blocks[start:start + c, start:start + c] = 1
        start += c

    ix = numpy.empty_like(voxel_communities)
    ix[voxel_communities.argsort()] = numpy.arange(n_voxels)
    blocks[:] = blocks[ix,:]
    blocks[:] = blocks[:,ix]
    return blocks

这里有一个简单的解释。这使用反向索引技巧将对角块数组的列和行重新排序为所需的矩阵。

    n_voxels = voxel_communities.size
    comm_labels = sorted(set(voxel_communities))
    comm_counts = [(voxel_communities == l).sum() for l in comm_labels]

    blocks = numpy.zeros((n_voxels, n_voxels), dtype=bool)
    start = 0
    for c in comm_counts:
        blocks[start:start + c, start:start + c] = 1
        start += c

这些行用于构造初始块矩阵。例如，假设您有六个体素和三个社区，每个社区包含两个体素。那么初始块矩阵将如下所示：

array([[ True,  True, False, False, False, False],
       [ True,  True, False, False, False, False],
       [False, False,  True,  True, False, False],
       [False, False,  True,  True, False, False],
       [False, False, False, False,  True,  True],
       [False, False, False, False,  True,  True]], dtype=bool)

这与所需的邻接矩阵基本相同在体素已按社区成员资格排序之后。所以我们需要反转排序。我们通过构造一个逆 argsort 数组来做到这一点。

    ix = numpy.empty_like(voxel_communities)
    ix[voxel_communities.argsort()] = numpy.arange(n_voxels)

现在ix用作索引时将反转排序过程。并且由于这是一个对称矩阵，我们可以分别对列和行进行反向排序操作：

    blocks[:] = blocks[ix,:]
    blocks[:] = blocks[:,ix]
    return blocks

这是它为小输入生成的结果示例：

>>> voxel_adjacency(numpy.array([0, 3, 1, 1, 0, 2]))
array([[ True, False, False, False,  True, False],
       [False,  True, False, False, False, False],
       [False, False,  True,  True, False, False],
       [False, False,  True,  True, False, False],
       [ True, False, False, False,  True, False],
       [False, False, False, False, False,  True]], dtype=bool)

在我看来，这与所建议的 voxel_matrix[np.ix_(voxels1, voxels2)] = 1 非常相似，只是它一次完成，而不是跟踪每个可能的节点组合。

可能有更好的解决方案，但这至少应该是一种改进。

此外，请注意，如果您可以简单地接受体素的新排序作为规范，那么此解决方案将变得与创建块数组一样简单！这需要大约 300 毫秒的时间。

根据索引和值填充 scipy / numpy 矩阵

Fill scipy / numpy matrix based on indices and values

python

numpy

matrix

scipy