基于指标值填充枕木/块状矩阵

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_matrixnp.ix_(voxels1，voxels2) =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)

在我看来，这与pv.所建议的pv.非常类似，只是它一次性完成，而不是跟踪每个可能的节点组合。

也许有一个更好的解决办法，但这至少应该是一个改进。

另外，请注意，如果您可以简单地接受新的体素排序为规范，那么这个解决方案就像创建块数组一样简单！这需要300毫秒左右。