Алгоритм нахождения порога перколяции во взвешенной сети

Чтобы расширить то, что предложил г-н Э., вот две версии алгоритма, которые прилично работают на графах с несколькими тысячами узлов. Обе версии используют 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, чтобы получить отсортированный массив всех вероятностей в матрице расстояний. Затем выполните стандартный бинарный поиск, как предложил г-н Э. На каждом шаге бинарного поиска замените элементы в матрице на 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, чтобы он останавливался всякий раз, когда достигается узел назначения. Код и тайминги выше были обновлены.

Однако не уверен, что правильно интерпретирую ваш вопрос:

Предположим, у вас есть пороговое значение кандидата, и вы хотите определить, существует ли путь между a и d. Вы можете проверить, какие узлы доступны из a, выполнив простой поиск в глубину на пороговом графе и проверив, был ли посещен желаемый конечный узел d.

Чтобы действительно найти порог, который вы знаете, он должен быть между нулем и максимальной вероятностью перехода на вашем графике (здесь 0,9). Таким образом, вы можете выполнить бинарный поиск порога, на каждом этапе используя поиск в глубину, чтобы проверить, есть ли у вас путь между a и d.