python - 在加权网络中查找 "percolation"阈值的算法

Question

我有一些通过转移概率链接的状态(嵌入在转移矩阵中，如马尔可夫链中)。我想通过仅考虑足够高的概率来总结此转换矩阵，以便它们允许从一个状态(〜节点)转到另一个状态(我的转换矩阵中的第一个和最后一个)。一个阈值，这样如果我只考虑更高的概率，我的转换矩阵永远不允许从第一个状态(或节点)移动到最后一个状态(或节点)。

两个问题:

是否有一些知名的库(最好是Python语言)实现了这种方法？我的天真/经验/原型(prototype)方法将是一个循环，它降低阈值，然后检查我是否可以从第一个状态流过过渡矩阵到最后一个状态(距离矩阵中的最佳路径算法？)。但这可能需要非常长的计算时间？

Example 1: 
DistMatrix = [[       'a',   'b',     'c',    'd'], 
             ['a',  0,      0.3,    0.4,    0.1],
             ['b',   0.3,    0,      0.9,    0.2],
             ['c',   0.4,    0.9,    0,      0.7],
             ['d',   0.1,    0.2,    0.7,   0]]  
    states are a,b,c,d. I want to find the value (threshold) that allow to go from a to d (no matter if other states are walked)  
    Naive approach:
    - first loop: threshold 0.9, I get rid of lesser probabilities: I can only connect c and b 
    - second loop: threshold 0.7, I get rid of lesser probabilities: I can only connect c, b, d
    - third loop: threshold 0.4, I get rid of lesser probabilities: I can connect a,c, b, d: here is my threshold: 0.4

--> 当我的转换矩阵有数千个状态时，应该会变得非常复杂？ --> 提出算法？

Example 2:
DistMatrix =
[       'a',   'b',     'c',    'd'],
['a',   0,      0.3,    0.4,    0.7],
['b',   0.3,    0,      0.9,    0.2],
['c',   0.4,    0.9,    0,      0.1],
['d',   0.7,    0.2,    0.1,    0] ] 
states are a,b,c,d. I want to find the value (threshold) that allow to go from a to d (no matter if other states are walked) 
Naive approach:
-first loop: threshold 0.9, I get rid of all others: I can only connect c and b
-second loop: threshold 0.7, I get rid of lesser connexion: I connect b and c, and a and d but because a and d are connected, I have my threshold!

Answer 1

为了扩展 E 先生的建议，这里有两个版本的算法，可以在具有数千个节点的图上正常工作。两个版本都使用Numpy第二个也使用 NetworkX .

您需要删除“a”、“b”、“c”和“d”标识符才能使用 Numpy 数组。通过将节点名称转换为 0 到 len(nodes) 之间的整数，可以轻松完成此操作。您的数组应如下所示

DistMatrix1 = np.array([[0,      0.3,    0.4,    0.1],
                        [0.3,    0,      0.9,    0.2],
                        [0.4,    0.9,    0,      0.7],
                        [0.1,    0.2,    0.7,   0]])

DistMatrix2 = np.array([[0,      0.3,    0.4,    0.7],
                        [0.3,    0,      0.9,    0.2],
                        [0.4,    0.9,    0,      0.1],
                        [0.7,    0.2,    0.1,    0]])

使用numpy.unique获取距离矩阵中所有概率的排序数组。然后，按照 E 先生的建议执行标准二分搜索。在二分搜索的每一步，如果矩阵中的条目低于当前概率，则将其替换为 0。在图上运行广度优先搜索，从第一个节点开始，然后查看是否到达最后一个节点。如果符合，则阈值较高，否则，阈值较低。 bfs代码实际上改编自NetworkX版本。

import numpy as np

def find_threshold_bfs(array):
    first_node = 0
    last_node = len(array) - 1
    probabilities = np.unique(array.ravel())
    low = 0
    high = len(probabilities)

    while high - low > 1:
        i = (high + low) // 2
        prob = probabilities[i]
        copied_array = np.array(array)
        copied_array[copied_array < prob] = 0.0
        if bfs(copied_array, first_node, last_node):
            low = i
        else:
            high = i

    return probabilities[low]


def bfs(graph, source, dest):
    """Perform breadth-first search starting at source. If dest is reached,
    return True, otherwise, return False."""
    # Based on http://www.ics.uci.edu/~eppstein/PADS/BFS.py
    # by D. Eppstein, July 2004.
    visited = set([source])
    nodes = np.arange(0, len(graph))
    stack = [(source, nodes[graph[source] > 0])]
    while stack:
        parent, children = stack[0]
        for child in children:
            if child == dest:
                return True
            if child not in visited:
                visited.add(child)
                stack.append((child, nodes[graph[child] > 0]))
        stack.pop(0)
    return False

速度较慢但较短的版本使用 NetworkX。在二分查找中，不运行 bfs，而是将矩阵转换为 NetworkX 图并检查是否存在从第一个节点到最后一个节点的路径。如果有路径，则阈值较高，如果没有路径，则阈值较低。这很慢，因为 NetworkX 中的所有图形数据结构的效率都比 Numpy 数组低得多。然而，它的优点是可以访问许多其他有用的算法。

import networkx as nx
import numpy as np

def find_threshold_nx(array):
    """Return the threshold value for adjacency matrix in array."""
    first_node = 0
    last_node = len(array) - 1
    probabilities = np.unique(array.ravel())
    low = 0
    high = len(probabilities)

    while high - low > 1:
        i = (high + low) // 2
        prob = probabilities[i]
        copied_array = np.array(array)
        copied_array[copied_array < prob] = 0.0
        graph = nx.from_numpy_matrix(copied_array)
        if nx.has_path(graph, first_node, last_node):
            low = i
        else:
            high = i

    return probabilities[low]

NetworkX 版本在具有一千多个节点的图表上崩溃(在我的笔记本电脑上)。 bfs版本可以轻松找到几千个节点的图的阈值。

代码运行示例如下。

In [5]: from percolation import *

In [6]: print('Threshold is {}'.format(find_threshold_bfs(DistMatrix1)))
Threshold is 0.4

In [7]: print('Threshold is {}'.format(find_threshold_bfs(DistMatrix2)))
Threshold is 0.7

In [10]: big = np.random.random((6000, 6000))

In [11]: print('Threshold is {}'.format(find_threshold_bfs(big)))
Threshold is 0.999766933071

对于时间安排，我得到(在半新的 Macbook Pro 上):

In [5]: smaller = np.random.random((100, 100))

In [6]: larger = np.random.random((800, 800))

In [7]: %timeit find_threshold_bfs(smaller)
100 loops, best of 3: 11.3 ms per loop

In [8]: %timeit find_threshold_nx(smaller)
10 loops, best of 3: 94.9 ms per loop

In [9]: %timeit find_threshold_bfs(larger)
1 loops, best of 3: 207 ms per loop

In [10]: %timeit find_threshold_nx(larger)
1 loops, best of 3: 6 s per loop

希望这有帮助。

更新

我修改了 bfs 代码，以便每当到达目标节点时它就会停止。上述代码和时序已更新。