Очень торопящиеся попасть из пункта А в пункт Б

Люди пытались найти более быстрые способы передвижения на протяжении всей своей истории. Появление качественной дорожной системы в римской империи в своё время привело к её расцвету, но со временем выяснилось, что и в продуманных дорожных системах бывают забавные изъяны, как например в небезызвестной задачи о кёнигсбергских мостах, считающейся отправной точкой возникновения теории графов. Неудивительно и то, что с развитием вычислительной техники логистические задачи стали одними из первых, над которыми трудились первопроходцы компьютерных наук. Задача о кратчайшем пути -- одна из них, звучит достаточно просто: есть несколько городов и дорог, соединяющих пару городов между собой, мы хотим попасть из города А в город Б пройдя при этом минимальное расстояние. Первый системный подход к этой задаче был описан в работе Эгервари в 1931г., спустя 25 лет Эдсгер Дейкстра придумал алгоритм, который сейчас является частью любого уважающего себя базового курса алгоритмов на графах. На нём же, будем честны, заканчиваются знания о кратчайших путях у большинства профессиональных разработчиков, ибо сценариев, где реализации с википедии/stackoverflow будет не хватать, крайне мало.

Жено Эгервари, один из первых теоретиков задачи о кратчайшем пути. Венгерский метод на самом деле метод Эгервари-Кёнига-Куна

Может показаться, что на самом деле просто не было существенного прогресса с 60х годов, так как Дейкстра предоставил почти асимптотически оптимальный алгоритм решения задачи. На самом деле нет, прогресс был и придумали много чего интересного, хоть и действительно с того времени фокус сместился на другие задачи. Приглашаю под кат если интересно узнать что такого напридумывали, что используется в современных логистических системах, почему меня огорчает отсутствие учёта флага единства в HOMM3 при расчёте пути, ну и наконец, что за мужики на картинке выше рядом с Дейкстрой?

Преамбула

Для начала определим задачу формально. Дан ориентированный граф , где -- множество вершин размера , -- множество рёбер размера , каждое ребро имеет длину , -- длина ребра из в если такое присутствует в графе.

Путь из в -- это последовательность вершин соседние вершины соединены ребром. Длина пути -- сумма длин ребёр этого пути. Задача о кратчайшем пути: для каждой пары источник-назначение найти путь наименьшей длины. Далее будем обозначать за длину кратчайшего пути из в .

Есть три основные разновидности задачи о кратчайших путях:

  • SSSP (single source shortest path): найти кратчайшие пути от одной вершины до всех остальных

  • APSP (all pair shortest path): найти кратчайшие пути от всех вершин до всех

  • P2P (point to point): найти кратчайший путь от одной вершины до другой

Если для задачи SSSP для каждой вершины зафиксировать по одному кратчайшему пути от истока до этой вершины и обозначить последнее ребро в выбранном кратчайшем пути из в , то подмножество рёбер образует дерево кратчайших путей -- структуру, позволяющую легко извлечь каждый отдельный путь при этом не храня их все непосредственно.

В дальнейшем в статье я буду использовать математические обозначения для псевдокода где обозначает присвоение в переменную , обращение к элементу ассоциативного массива/словаря обозначается круглыми скобками вместо более привычных квадратных , но также будет присутствовать и более привычный код на Python, для избежания лишнего ООП граф представляется в виде 

Dict[int, Dict[int, float]]

В частности graph[u][v] соответствует длине.

Ну и наконец будут картинки с визуализацией на OSM (надеюсь после этого все будут знать как в Санкт-Петербурге дойти от Владимирской до Юбилейного)

Базовый пример в OSM
import osmnx

G = osmnx.graph_from_point((59.93893094417527, 30.32268115454809), dist=2500)
s = 365356365
t = 8518165879
shortest_path = osmnx.routing.shortest_path(G, s, t, weight='length')
osmnx.plot.plot_graph_route(G,
                            shortest_path,
                            route_color='g',
                            node_size=2,
                            node_color='red',
                            edge_color='b',
                            edge_linewidth=0.2,
                            bgcolor='white')

Классические алгоритмы: Scanning method, Форд-Беллман, Уоршелл-Флойд, Дейкстра

Лестер Форд младший, его вы можете знать по алгоритмам Форда-Беллмана и Форда-Фалкерсона

Один из первых значимых алгоритмов для SSSP был предложен Фордом в его работах по сетевым потокам. Из-за его простоты называли его обычно либо scanning method, либо labelling method (примечание: к сожалению в рускоязычной литературе этот метод практически не встречается, поэтому оставил англоязычное название). Метод хранит два непересекающихся множества вершин: помеченные,и обработанные, для каждой вершины хранится текущая оценка минимального расстояния и последнее ребро в минимальном пути . Изначально есть только одна помеченная вершина и ни одной обработанных, расстояния инициализируются как . Далее происходит следующее

  • : если то обновить и поместить вершину в помеченные

  • : применить для и поместить вершинув обработанные

  • Основной цикл: пока есть помеченная вершина просканировать её

Код scanning method
from collections import Counter

def relax(graph: Dict[int, Dict[int, float]],
          u: int,
          v: int,
          distances: Dict[int, float],
          prev: Dict[int, int],
          labelled: Set,
          stats: Counter = None):
    if stats is not None:
        stats["relax_calls"] += 1
    if v not in distances or distances[v] > distances[u] + graph[u][v]:
        if stats is not None:
            stats["distance_updates"] += 1
        distances[v] = distances[u] + graph[u][v]
        prev[v] = u
        labelled.add(v)

def scan(graph: Dict[int, Dict[int, float]],
         node: int,
         distances: Dict[int, float],
         prev: Dict[int, int],
         labelled: Set,
         stats: Counter = None):
    if stats is not None:
        stats["scan_calls"] += 1
        
    for other, length in graph[node].items():
        relax(graph, node, other, distances, prev, labelled, stats)
    # labelled.remove(node)

def scanning_method(graph: Dict[int, Dict[int, float]],
                    source: int,
                    target: int = None,
                    stats: Counter = None):
    distances = {source: 0}
    prev = dict()
    labelled = {source}
    while len(labelled) > 0:
        u = labelled.pop()
        if u == target:
            break
        scan(graph, u, distances, prev, labelled, stats)

    return distances, prev

Этот метод гарантированно сходится при условии существования решения, однако при неудачном выборе порядка вершин для сканирования имеет экспоненциальную сложность. Одной из уникальных особенностей этого метода является то, что подавляющее число появившихся после методов являются его частными случаями и различаются только способом выбора вершины для сканирования.

Доказательство корректности scanning method

  1. Расстояния не увеличиваются по ходу работы алгоритма

  2. Если в какой-то момент , то существует путь от до с расстоянием . При отсутствии отрицательных циклов этот путь простой (вершины в нем не повторяются)

  3. Остановка: каждый раз когда мы помечаем вершину текущий кратчайший путь до неё изменяется, так как путей конечное число, то и количество раз когда вершина оказывается помеченной тоже ограничено

  4. Если вершина находится в обработанных, то для всех соседних вершин так как после последнего сканирования неравенство выполнялось, а могло только уменьшиться

  5. Для выполняется так как при обновлении расстояния меняется и образуя равенство, а может только уменьшиться

  6. Корректность: допустим после остановки алгоритма для нескольких вершин расстояние посчиталось неправильно. Тогда существует пара вершин, что для первой расстояние неправильное, для второй -- правильное (вторая должна найтись ибо для расстояние посчитано правильно) и при этом минимальный путь от до проходит через . Тогда что противоречит пункту 4 критерию остановки.

Ричард Беллман, отец динамического программирования

Корректность и гарантия остановки -- это хорошо, но хочется при этом гарантии, что алгоритм не будет работать слишком долго, да и вообще хочется понять насколько быстро будет работать алгоритм. Во многом первичный прогресс как по анализу сложности, так и по разработке эффективных методов во многом основан на динамическом программировании -- общей техники для ряда оптимизационных задач, предложенный Ричардом Беллманом. Алгоритм Форда-Беллмана -- первый из алгоритмов, построенном на этом принципе, он делится на несколько стадий. На стадии вычисляются расстояния на путях, состоящих из не более, чем рёбер. Если обозначить за минимальное расстояние от до используя не более, чем шагов, то при и, самое главное, . Сам алгоритм по сути представляет собой последовательное (рекуррентное) вычисление начиная с малых , каждая стадия занимает , всего стадий , что даёт оценку Последний штрих алгоритма -- мало того, что посчитав полностью значения для можно забыть о значениях для так еще и можно это делать в одном массиве не храня отдельно расстояния для двух последовательных стадий алгоритма.

Алгоритма Форда-Беллмана
def bellman_ford_moore(graph: Dict[int, Dict[int, float]],
                       source: int) -> Tuple[Dict[int, float], Dict[int, int]]:
    distances = {source: 0}
    prev = dict()
    changed = True
    for k in graph:
        if not changed:
            break
        changed = False
        for u, adjacent in graph.items():
            if u not in distances:
                continue
            for v, l in adjacent.items():
                if v not in distances or distances[v] > distances[u] + l:
                    distances[v] = distances[u] + l
                    prev[v] = u
                    changed = True

    return distances, prev

Позднее Мур заметил, что можно избавиться от лишних действий переписав его в виде scanning method с FIFO порядком обработки вершин (иногда этот метод называют shortest path faster algorithm SPFA , в англоязычной литературе алгоритм обычно обозначается как BFM от Bellman-Ford-Moore).

Эдсгер Дейкстра, автор самого распространённого алгоритма для SSSP

Довольно быстро идеи динамического программирования подхватили и остальные, больше всего выделяются два: алгоритм Дейкстры для SSSP и алгоритм Флойда-Уоршелла для APSP. Дейкстра заметил следующее: если отсортировать все вершины графа в порядке отдаления от, т.е. последовательность такова, что кратчайшие расстояния удовлетворяют, то в случае отсутствия рёбер отрицательного веса кратчайший путь от до может проходить только по вершинам . Если обозначить -- длину кратчайшего пути, имеющего в качестве промежуточных вершин только , то и при этом . Таким образом чтобы получить из достаточно просканировать вершину расстояние до которой корректно посчитано на шаге и при этом она определяется как вершина с наименьшим расстоянием еще не отсканированных на предыдущих шагах. В терминах сканирующего метода алгоритм Дейкстры -- это выбор вершины с наименьшим текущим расстоянием для очередного сканирования. При наличии способа хранения текущих расстояний с возможностью уменьшения и извлечения минимума со сложностью алгоритм Дейкстры являлся бы оптимальным для SSSP так как работал бы за .

Роберт Тарьян, один из самых главных выдумщиков структур данных, в том числе кучи Фибоначчи

К сожалению такой структуры так и не придумали, из наиболее применимых на практике -- двоичная куча, с которой алгоритм Дейкстры работает за и куча Фибоначчи с декрементом за амортизированную константу и удалением минимума за логарифм, что дает .

Алгоритм Дейкстры
def dijkstra(graph: Dict[int, Dict[int, float]],
             source: int,
             target: int = None,
             stats: Counter = None) -> Tuple[Dict[int, float], Dict[int, int]]:
    distances = {source: 0}
    prev = dict()
    labelled = {source}
    while len(labelled) > 0:
        u = min(labelled, key=lambda x: distances[x])
        labelled.remove(u)
        if u == target:
            break
        scan(graph, u, distances, prev, labelled, stats)

    return distances, prev

Роберт Флойд, придумал самый короткий алгоритм для APSP

Еще одним интересным результатом является алгоритм Флойда-Уоршелла также основывающийся на динамическом программировании. Если предположить, что все вершины пронумерованы и обозначить -- минимальное расстояние от до , имеющее в качестве промежуточных вершин только вершины с номерами не более , то получается, что и при этом .

Немного про линейное программирование

Слева Джордж Данциг, будучи студентом случайно доказал две нерешенные проблемы в статистике, думая, что это ДЗ, придумал симплекс-метод и ... алгоритм Дейкстры. Справа Александр Мадри, не так давно представил серию статей с использованием IPM для APSP.

Задача о кратчайших путях является частным случаем задачи потока минимальной стоимости, которая в свою очередь является частным случаем задачи линейного программирования. Долгое время основным методом решения общих задач линейного программирования был симплекс метод, который не был эффективен для таких задач как поиск кратчайшего пути был слишком громоздким и все перечисленные выше методы были лучше него. Однако сначала в 70х-80х годах Нестеровым и Немировским был разработан метод внутренней точки (interior point method, IPM), который был быстр на практике, и при имел полиномиальную оценку сложности. Позднее этот метод стал основным для решения не только задач линейного программирования, но и более общих. Совсем недавно на его основе были получены ряд улучшений для потоковых задач и, косвенно, для APSP.

Эти работы касаются продвинутой математики, мы её опустим в этой статье, но пару вещей стоит разобрать. Во-первых, если -- матрица инцидентности графа, -- вектор такой, что , то P2P задача о кратчайшем пути из в формулируется в виде следующей задачи линейного программирования

Интерпретация этой задачи следующая: мы пытаемся пустить единицу потока из в , пропустить единицу потока через стоит , оптимальное решение -- пустить поток по кратчайшему пути. Двойственная задача в терминах линейного программирования имеет следующий вид

Интерпретация этой задачи следующая: представьте, что ваш граф состоит из небольших шариков (вершины), связанных нитками (рёбра) соответствующей длины. Если потянуть влево, а вправо, на какое максимальное расстояние растянется конструкция? Точно также на расстояние минимального пути от до . Величины , удовлетворяющие ограничениям в двойственной задаче обычно называются допустимыми потенциалами и играют важную роль в последующих алгоритмах.

Ключевой трюк с потенциалами: рассмотрим модифицированные длины

Длина пути в терминах нового расстояния

Таким образом пути, соединяющие одну и ту же пару вершин, отличаются на константу, из чего следует, что кратчайший путь относительно является также и кратчайшим путём относительно . Такие длины называются приведёнными. Интересен следующий факт: если -- допустимый потенциал по отношению к , то приведённые длины неотрицательны. Более того, кратчайшие расстояния от одной вершины до других является допустимыми потенциалами (так как является решением двойственной задачи указанной выше). На этих свойствах построен алгоритм Джонсона для APSP: применить алгоритм Форда-Беллмана для одной из вершин, для всех остальных применить алгоритм Дейкстры использую приведённые длины относительно потенциалов, полученных из расстояний до первой вершины. Этот алгоритм выигрывает у Уоршелла-Флойда в случае, когда .

Усеченный и двунаправленный алгоритмы Дейкстры для P2P

Несмотря на то, что оценка алгоритма Дейкстры с идеальной структурой данных для извлечения минимума имеет неулучшаемую оценку для SSSP, в случае P2P есть одно простое улучшение, позволяющее не просматривать граф целиком: если нам нужно найти кратчайший путь от до , то можно остановиться если на очередной итерации мы сканируем .

Пример на OSM
def graph_from_osmnx(G):
    graph = dict()
    for n, adj in G.adjacency():
        if n not in graph:
            graph[n] = dict()
        for e, eattr in adj.items():
            for _, iattr in eattr.items():
                if e not in graph[n] or graph[n][e] > iattr["length"]:
                    graph[n][e] = iattr["length"]
    return graph

g = graph_from_osmnx(G)
distances, prev = dijkstra(g, s, t)

def print_route_and_visited(graph, source, target, distances, prev):
    route = [target]
    cur = target
    while cur != source:
        cur = prev[cur]
        route.append(cur)
    return osmnx.plot.plot_graph_route(G.subgraph(list(distances.keys())), list(reversed(route)), route_color='g', node_size=2, node_color='r', edge_color='b', edge_linewidth=0.2, bgcolor='white', close=True);

print_route_and_visited(G, s, t, distances, prev)
Вершина, находящиеся дальше от точки назначения, не рассматривались

Двунаправленный алгоритм Дейкстры -- применение принципа meet-in-the-middle: пытаемся параллельно искать пути с двух концов, от в исходном графе и от в графе, в котором развернуты рёбра, останавливаемся когда на очереди сканирования вершина, которая уже была отсканирована в обратном порядке. Общий смысл в том, что обычно количество вершин на расстоянии от больше, чем количество вершин на расстоянии от или . У двунаправленного Дейкстры есть важный нюанс: если в какой-то момент мы нашли вершину , которую просканировали в прямом графе от и в обратном графе от , то все вершины на кратчайшем пути мы уже нашли, но не обязана лежать на кратчайшем пути. В примере ниже таковой вершиной будет являться , но кратчайший путь -- .

Двунаправленный алгоритм Дейкстры не так прост как кажется

Обозначим за минимальные расстояния до вершин в прямом и обратном проходе соответственно. Если же мы на каждом шаге отслеживаем величину

То указанный выше критерий остановки становится корректным, т.е. указанная величина будет соответствовать кратчайшему пути.

Доказательство

Предположим, что -- пройденные в прямом проходе, -- пройденные в обратном проходе. Из критерия остановки следует, что существует вершина, который присутствует в обоих, обозначим её за . Кратчайший путь должен начинаться с нескольких и заканчиваться несколькими . Предположим, что в любом кратчайшем пути есть хотя бы одна вершина не из этого списка, но тогда для кратчайших расстояний выполняется и получаем противоречие. Таким образом есть хотя бы один кратчайший путь, состоящий из просканированных вершин, а указанная выше величина -- минимальная длина таких путей.

Код двунаправленного алгоритма Дейкстры
def reversed_graph(graph: Dict[int, Dict[int, float]]) -> Dict[int, Dict[int, float]]:
    rev_graph = {v: dict() for v in graph}
    for v, adj in graph.items():
        for u, l in adj.items():
            rev_graph[u][v] = l
    return rev_graph

def bidir_dijkstra(graph: Dict[int, Dict[int, float]],
                   source: int,
                   target: int = None,
                   stats: Counter=None):
    rev_graph = reversed_graph(graph)
    distances_f = {source: 0}
    prev_f = dict()
    distances_b = {target: 0}
    prev_b = dict()
    scanned_f = set()
    scanned_b = set()
    
    labelled_f = {source}
    labelled_b = {target}

    best = float("inf")
    best_f = source
    best_b = target

    u = source
    v = target
    
    while len(labelled_f) > 0 or len(labelled_b):
        if distances_f[u] < distances_b[v]:
            if u in scanned_b:
                break
            labelled_f.remove(u)
            scan(graph, u, distances_f, prev_f, labelled_f, stats)
            for p, l in graph[u].items():
                if p in distances_b and distances_f[u] + distances_b[p] + l < best:
                    best_f = u
                    best_b = p
                    best =  distances_f[u] + distances_b[p] + l
            scanned_f.add(u)
            u = min(labelled_f, key=lambda x: distances_f[x])
        else:
            if v in scanned_f:
                break
            labelled_b.remove(v)
            scan(rev_graph, v, distances_b, prev_b, labelled_b, stats)
            for p, l in rev_graph[v].items():
                if p in distances_f and distances_b[v] + distances_f[p] + l < best:
                    best_f = p
                    best_b = v
                    best =  distances_b[v] + distances_f[p] + l
            scanned_b.add(v)
            v = min(labelled_b, key=lambda x: distances_b[x])

    best_path = [best_f]
    cur = best_f
    while cur != source:
        cur = prev_f[cur]
        best_path.append(cur)
    best_path = list(reversed(best_path))

    cur = best_b
    if cur != best_path[-1]:
        best_path.append(cur)
    while cur != target:
        cur = prev_b[cur]
        best_path.append(cur)
    
    return distances_f, distances_b, prev_f, prev_b, best_path, best

Пример на OSM
distances_f, distances_b, prev_f, prev_b, best_path, best = bidir_dijkstra(g, s, t)
osmnx.plot.plot_graph_route(G.subgraph(list(distances_f.keys()) + list(distances_b.keys())),
                            best_path,
                            route_color='g',
                            node_size=2,
                            node_color='r',
                            edge_color='b',
                            edge_linewidth=0.2,
                            bgcolor='white')

A*

Питер Харт и Нильс Нильсон. Придумали как улучшить алгоритм Дейкстры и назвали его A*

Как мы обсудили раньше, задачу о кратчайших путях можно эквивалентным образом преобразовать с помощью потенциалов. Джонсон использовал это для того, чтобы была возможность применить алгоритм Дейкстры к APSP c отрицательными рёбрами. Для APSP/SSSP этот трюк к сожалению кроме неотрицательности рёбер ничего дать не может. А вот с P2P ситуация иная: при использовании разных потенциалов усеченный алгоритм Дейкстры может просматривать разное количество вершин.

Теорема 1. Если допустимые потенциалы таковы, что и, то множество пройденных алгоритмом Дейкстры вершин с потенциалами является подмножеством вершин, пройденных алгоритмом Дейкстры с потенциалами .

Доказательство

По сути нужно доказать, что если длина кратчайшего пути меньше длины кратчайшего пути относительно то она меньше и относительно Если вершина была пройдена с потенциалами , то обозначив минимальный путь, а за минимальный путь получаем

Это в частности доказывает, что 1) если применить нетривиальные потенциалы к SSSP с неотрицательными весами, то усеченный Дейкстра с потенциалами будет не хуже, чем обычный; 2) Чем лучше потенциалы, тем лучше алгоритм.

Отметим, что при если -- допустимые потенциалы и то является оценкой снизу минимального расстояние от до (но не наоборот). Алгоритм Дейкстры с приведенными длинами эквивалентен альтернативному правилу выбора очередной вершины для сканирования: вместо используется -- оригинальная формулировка алгоритма A*.

Пример с OSM, потенциалы взяты как евклидово (если быть точнее геодезическое) расстояние
import math
import geopy.distance

def geo_dist(graph, u, v):
    coords_1 = (graph.nodes[u]['y'], graph.nodes[u]['x'])
    coords_2 = (graph.nodes[v]['y'], graph.nodes[v]['x'])

    return geopy.distance.distance(coords_1, coords_2).m

def a_star(graph: Dict[int, Dict[int, float]],
           source: int,
           target: int,
           graphnx,
           stats: Counter=None) -> Tuple[Dict[int, float], Dict[int, int]]:
    distances = {source: 0}
    prev = dict()
    potentials = dict() 
    labelled = {source}
    while len(labelled) > 0:
        u = -1
        best = -1
        for v in labelled:
            if v not in potentials:
                potentials[v] = geo_dist(graphnx, v, target)
            if u == -1 or distances[v] + potentials[v] < best:
                best = potentials[v] + distances[v]
                u = v
        labelled.remove(u)
        if u == target:
            break
        scan(graph, u, distances, prev, labelled, stats)

    return distances, prev

distances, prev = a_star(g, s, t, G)
print_route_and_visited(G, s, t, distances, prev)

Насколько же хорош может быть A*? Как оказывается, с идеальными потенциалами мы можем получить оптимальный алгоритм.

Теорема 2. Если -- длина кратчайшего пути, то A* c потенциалами просканирует только вершины, лежащие на каком-либо кратчайшем пути.

Доказательство

Если лежит на кратчайшем пути, то его приведенная стоимость равна нулю, если же не лежит, то строго больше нуля. А* в первую очередь обойдет эти рёбра и дойдет по ним до

Несмотря на то, что теорема может показаться бессмысленной, так как если у нас уже есть кратчайшие расстояние до то и задачи поиска кратчайшего пути нет. На самом деле не совсем, во-первых есть ситуации когда сами расстояния есть, но пути нет, но с помощью А* его можно быстро найти. Во-вторых, учитывая теорему 1 мы получаем, что чем точнее потенциалы, тем ближе мы к теоретически оптимальному алгоритму, который обходит только вершины на кратчайшем пути.

Двунаправленный A*

Хироши Имай, единственный соавтор двунаправленного А*, фото которого удалось найти

Как и с двунаправленным Дейкстрой есть небольшая дополнительная подковырка в случае с А*. Дело в том, что если взять два произвольных потенциала для прямого и обратного прохода, то критерий остановки двунаправленного Дейкстры становится некорректным. Критерий будет корректным, если потребовать от потенциалов для прямого и обратного прохода . Если у нас есть потенциалы то удовлетворяют указанным условиям и являются допустимыми потенциалами

Код двунаправленного А* и пример с OSM
def bidir_a_star(graph: Dict[int, Dict[int, float]],
                 source: int,
                 target: int,
                 graphnx,
                 stats: Counter=None):
    rev_graph = reversed_graph(graph)
    distances_f = {source: 0}
    prev_f = dict()
    distances_b = {target: 0}
    prev_b = dict()
    scanned_f = set()
    scanned_b = set()
    
    labelled_f = {source}
    labelled_b = {target}

    best = float("inf")
    best_f = source
    best_b = target

    u = source
    v = target
    potentials = dict()
    
    while len(labelled_f) > 0 or len(labelled_b):
        if distances_f[u] < distances_b[v]:
            if u in scanned_b:
                break
            labelled_f.remove(u)
            scan(graph, u, distances_f, prev_f, labelled_f, stats)
            for p, l in graph[u].items():
                if p in distances_b and distances_f[u] + distances_b[p] + l < best:
                    best_f = u
                    best_b = p
                    best =  distances_f[u] + distances_b[p] + l
            scanned_f.add(u)
            for x in labelled_f:
                if x not in potentials:
                    potentials[x] = 0.5 * (geo_dist(graphnx, x, target) - geo_dist(graphnx, x, source))
            u = min(labelled_f, key=lambda x: distances_f[x] + potentials[x])
        else:
            if v in scanned_f:
                break
            labelled_b.remove(v)
            scan(rev_graph, v, distances_b, prev_b, labelled_b, stats)
            for p, l in rev_graph[v].items():
                if p in distances_f and distances_b[v] + distances_f[p] + l < best:
                    best_f = p
                    best_b = v
                    best =  distances_b[v] + distances_f[p] + l
            scanned_b.add(v)
            for x in labelled_b:
                if x not in potentials:
                    potentials[x] = 0.5 * (geo_dist(graphnx, x, target) - geo_dist(graphnx, x, source))
            v = min(labelled_b, key=lambda x: distances_b[x] - potentials[x])

    best_path = [best_f]
    cur = best_f
    while cur != source:
        cur = prev_f[cur]
        best_path.append(cur)
    best_path = list(reversed(best_path))

    cur = best_b
    if cur != best_path[-1]:
        best_path.append(cur)
    while cur != target:
        cur = prev_b[cur]
        best_path.append(cur)
    
    return distances_f, distances_b, prev_f, prev_b, best_path, best
fig, ax = osmnx.plot.plot_graph_route(G.subgraph(list(distances_f.keys()) + list(distances_b.keys())),
                                      best_path,
                                      route_color='g',
                                      node_size=2,
                                      node_color='r',
                                      edge_color='b',
                                      edge_linewidth=0.2,
                                      bgcolor='white')

Бонус: несколько кратчайших путей и граф почти оптимальных путей

Есть две интересные вспомогательные вещи, которые можно делать с помощью потенциалов и двунаправленного прохода.

Первая -- найтилучших путей. Как уже обсуждали раньше, если запустить Дейкстру с потенциалами в точности соответствующими кратчайшему расстоянию до то он просканирует только вершины на кратчайшем пути. Если же чуть модифицировать процедуру позволив вершины повторно обрабатываться с бОльшим весом, то получается алгоритм, находящий лучших путей, а не один. Вот например старая статейка от разработчиков майкрософт, там этот метод используется для быстрого поиска самых частотных продолжений запросов

Код
import heapq

def k_shortest_path_looply(graph: Dict[int, Dict[int, float]],
                           source: int,
                           target: int,
                           k: int,
                           min_loop_length: int = 1) -> List[List[int]]:
    distances_b, _ = dijkstra(reversed_graph(graph), target)
    result = []
    sorted_adjacency = dict() # Dict[List[Tuple[int, float]]]
    entries = [(source, -1, 0)] # vertex_id, parent entry, sorted adjacency position
    path_queue = [(0, 0)] # reduced length, -entry_id 

    while len(result) < k:
        length, entry_id_rev = heapq.heappop(path_queue)
        entry_id = -entry_id_rev
        v, parent_entry, pos = entries[entry_id]
        if v == target:
            cur = entry_id
            path = [target]
            while cur != 0:
                cur = entries[cur][1]
                path.append(entries[cur][0])
            result.append(list(reversed(path)))

    
        if v not in sorted_adjacency:
            sorted_adjacency[v] = list(sorted(graph[v].items(), key=lambda x: x[1] + distances_b[x[0]]))
        if parent_entry != -1 and pos < len(sorted_adjacency[entries[parent_entry][0]]) - 1:
            u, length = sorted_adjacency[entries[parent_entry][0]][pos + 1]
            cur = parent_entry            
            if not has_cycle:    
                entries.append((u, parent_entry, pos + 1))
                heapq.heappush(path_queue, (length + distances_b[u] - distances_b[entries[parent_entry][0]], 1 - len(entries)))

        if len(sorted_adjacency) > 0:
            u, length = sorted_adjacency[v][0]
            has_cycle = False
            cur = entry_id
            for i in range(min_loop_length - 1):
                if cur == 0:
                    break
                if entries[cur][0] == u:
                    has_cycle = True
                    break
                cur = entries[cur][1]
            if not has_cycle:
                entries.append((u, entry_id, 0))
                heapq.heappush(path_queue, (length + distances_b[u] - distances_b[v], 1 - len(entries)))

    return result

К сожалению эта версию допускает циклы и часто на реальных графах, побочные пути -- это кратчайший путь с небольшим циклом.

Вторая -- найти подграф, который содержит все пути, которые хуже, чем кратчайший, не больше, чем на Если -- прямые и обратные расстояния, то критерий того, что ребро находится в подграфе

Код и пример OSM
def best_path_subgraph(graph: Dict[int, Dict[int, float]],
                       source: int,
                       target: int,
                       beam: float) -> List[Tuple[int, int]]:
    distances_f, _ = dijkstra(graph, source)
    distances_b, _ = dijkstra(reversed_graph(graph), target)
    best = distances_f[target]

    edges = []
    for u, adj in graph.items():
        if u not in distances_f:
            continue
        for v, length in adj.items():
            if v not in distances_b:
                continue
            if distances_f[u] + length + distances_b[v] < best + beam:
                edges.append((u, v))
    return edges

  edges = best_path_subgraph(g, s, t, 500)
  fig, ax = osmnx.plot.plot_graph(G.edge_subgraph([(u, v, 0) for u, v in edges]),
                                node_size=2,
                                node_color='r',
                                edge_color='b',
                                edge_linewidth=0.2,
                                bgcolor='white')

ALT

Эндрю Гольдберг, живой классик и по совместительству главный ответственный за логистику в амазон.

Последний в этой статье (но не последний известный науке) алгоритм -- продолжение концепции потенциалов и A*. В примерах я использовал естественным образом возникающее евклидово расстояние в качестве потенциалов, но остаётся вопрос, можно ли придумать потенциалы лучше? На практике PTP задача применяется в формате, что у нас есть очень много запросов на большом графе, предпосчитать их возможности нет, но есть возможность сделать какой-то более простой предпосчет. В целом можно считать, что в P2P мы на самом деле хотим структуру данных, которая тратит времени при инициализации и находить кратчайший путь между двумя вершинами за A* в целом является таковым, но при этом не делает никакого предпосчета. Возникает ощущение, что все-таки за счет какого-то предподсчета можно получить выигрышь на производительности запросов. И вот мы приходит к алгоритму ALT (A*, landmarks, triangle inequality): идея в том, чтобы предпосчитать кратчайшие пути до какой-то конкретной вершины и использовать потенциал из неравенства треугольника

Такие вершины называют ориентирами (landmarks, не уверен есть ли официальный рускоязычный термин). Хорошие ориентиры получаются из вершин, находящихся в отдалении. В примерах, которые я рассматривал, хорошим ориентиром будет вершина на границе области. Последний штрих -- это возможность использовать несколько ориентиров

Теорема 3. Если -- допустимые потенциалы, то -- тоже допустимый потенциал.

Доказательство

применяем максимум к левой и правой части, получаем то, что нужно.

Код
def alt(graph: Dict[int, Dict[int, float]],
        source: int,
        target: int,
        landmark_distances,
        stats: Counter=None):
    rev_graph = reversed_graph(graph)
    distances_f = {source: 0}
    prev_f = dict()
    distances_b = {target: 0}
    prev_b = dict()
    scanned_f = set()
    scanned_b = set()
    
    labelled_f = {source}
    labelled_b = {target}

    best = float("inf")
    best_f = source
    best_b = target

    u = source
    v = target
    potentials = dict()
    
    while len(labelled_f) > 0 or len(labelled_b):
        if distances_f[u] < distances_b[v]:
            if u in scanned_b:
                break
            labelled_f.remove(u)
            scan(graph, u, distances_f, prev_f, labelled_f, stats)
            for p, l in graph[u].items():
                if p in distances_b and distances_f[u] + distances_b[p] + l < best:
                    best_f = u
                    best_b = p
                    best =  distances_f[u] + distances_b[p] + l
            scanned_f.add(u)
            for x in labelled_f:
                if x not in potentials:
                    potential_f = 0
                    potential_b = 0
                    for p in landmark_distances:
                        potential_f = max(potential_f, landmark_distances[p][target] - landmark_distances[p][x])
                        potential_b = max(potential_b, landmark_distances[p][x] - landmark_distances[p][source])                  
                    potentials[x] = 0.5 * (potential_f - potential_b)
                        
            u = min(labelled_f, key=lambda x: distances_f[x] + potentials[x])
        else:
            if v in scanned_f:
                break
            labelled_b.remove(v)
            scan(rev_graph, v, distances_b, prev_b, labelled_b, stats)
            for p, l in rev_graph[v].items():
                if p in distances_f and distances_b[v] + distances_f[p] + l < best:
                    best_f = p
                    best_b = v
                    best =  distances_b[v] + distances_f[p] + l
            scanned_b.add(v)
            for x in labelled_b:
                if x not in potentials:
                    potential_f = 0
                    potential_b = 0
                    for p in landmark_distances:
                        potential_f = max(potential_f, landmark_distances[p][target] - landmark_distances[p][x])
                        potential_b = max(potential_b, landmark_distances[p][x] - landmark_distances[p][source])                  
                    potentials[x] = 0.5 * (potential_f - potential_b)
            v = min(labelled_b, key=lambda x: distances_b[x] - potentials[x])

    best_path = [best_f]
    cur = best_f
    while cur != source:
        cur = prev_f[cur]
        best_path.append(cur)
    best_path = list(reversed(best_path))

    cur = best_b
    if cur != best_path[-1]:
        best_path.append(cur)
    while cur != target:
        cur = prev_b[cur]
        best_path.append(cur)
    
    return distances_f, distances_b, best_path, best

Пример на OSM
Точки по периметру -- использованные ориентиры

Бонус: pathfinding в играх

Battle for Wesnoth

В опенсорс игре Battle for Wesnoth юниты передвигаются по гексагональному миру, каждый юнит может передвинуться на несколько хексов за ход, некоторые хексы требуют большое очков передвижения. В общем, ничего особенного, для поиска пути используется А*.

Heroes of might and magic 3

Heroes of might and magic 3 карта делится на квадратные тайлы, передвигаться можно либо по горизонтали (100 очков передвижения), либо по диагонали (141 очков передвижения). Также есть штрафы за местность. В общем то тоже ничего особенного, но в игре есть объекты, которые увеличивают количество очков передвижения (например флаг единства) и pathfinding не учитывает эти объекты, а казалось бы всего-то надо научиться решать задачу о кратчайших путях с отрицательными рёбрами.

Warcraft 2000

Стратегии реального времени как оказывается тоже для поиска пути делят карту на квадратные тайлы и используют специальные версии A*. В более старых стратегиях это приводило к тому, что юниты передвигались горизонтально-вертикально-диагонально, но не по прямой.

Dota 2

В общем то где-то между вторым и третьим варкрафтом таки научились ходить по прямой, зачастую просто с помощью сглаживания пути после его нахождения с помощью А*.

Итог

В статью не вошли свежие методы такие как contraction hierarchies и hub labeling потому что я заколебался писать. Для показанных методов есть статистика производительности по число вызванных scan, relax и количество обновлений кратчайших путей. Эта статистика не может считаться полноценным бенчмарком, но в целом даёт хороший ориентир

+----------------+-------------+------------------+------------+
|     Method     | Relax calls | Distance updates | Scan calls |
+----------------+-------------+------------------+------------+
|     basic      |   1060170   |      437690      |   385509   |
|    dijkstra    |    81090    |      34065       |   29730    |
| bidir_dijkstra |    52701    |      22416       |   19307    |
|     a_star     |    13100    |       5904       |    4764    |
|  bidir_a_star  |     8858    |       4160       |    3203    |
|      alt       |     1231    |       671        |    463     |
+----------------+-------------+------------------+------------+
Список литературы

Bellman, R. (1958). On a routing problem. Quarterly of Applied Mathematics, 16(1), 87–90.

Cohen, M. B., Mądry, A., Sankowski, P., & Vladu, A. (2017). Negative-weight shortest paths and unit capacity minimum cost flow in õ (m 10/7 log w) time*. Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms, 752–771.

Dantzig, G. B. (1962). Linear programming and extensions.

Dial, R. B. (1969). Algorithm 360: shortest-path forest with topological ordering [H]. Commun. ACM, 12(11), 632–633. https://doi.org/10.1145/363269.363610

Dijkstra, E. W. (1959). A note on two problems in connexion with graphs. Numerische Mathematik, 50, 269–271.

Doran, J. (1967). An approach to automatic problem-solving. Machine Intelligence, 1, 105–127.

Edmonds, J., & Karp, R. M. (1972). Theoretical improvements in algorithmic efficiency for network flow problems. Journal of the ACM (JACM), 19(2), 248–264.

Egerváry, E. (1931). On combinatorial properties of matrices.

Floyd, R. W. (1962). Algorithm 97: Shortest path. Commun. ACM, 5(6), 345. https://doi.org/10.1145/367766.368168

Ford, L. R. (1956). Network flow theory.

Ford, L. R., & Fulkerson, D. (1962). Flows in Networks (Vol. 56). Princeton University Press.

Goldberg, A. (2005). Basic shortest path algorithms. DIKU Summer School on Shortest Paths. Microsoft Research. Silicon Valey.

Goldberg, A. V., & Harrelson, C. (2005). Computing the shortest path: A search meets graph theory. SODA, 5, 156–165.

Hart, P. E., Nilsson, N. J., & Raphael, B. (1968). A Formal Basis for the Heuristic Determination of Minimum Cost Paths. IEEE Transactions on Systems Science and Cybernetics, 4(2), 100–107. https://doi.org/10.1109/TSSC.1968.300136

Hsu, B.-J., & Ottaviano, G. (2013). Space-efficient data structures for top-k completion. Proceedings of the 22nd International Conference on World Wide Web, 583–594.

Ikeda, T., & Imai, H. (1994). Fast a algorithms for multiple sequence alignment. Genome Informatics, 5, 90–99.

Johnson, D. B. (1977). Efficient algorithms for shortest paths in sparse networks. Journal of the ACM (JACM), 24(1), 1–13.

Le, P. T. H., Truong, N. T. N., Kim, M., So, W., & Jung, J. H. (2018). Applying Theta* in Modern Game. J. Comput., 13(5), 527–536.

Moore, E. F. (1959). The shortest path through a maze. Proc. of the International Symposium on the Theory of Switching, 285–292.

Nicholson, T. A. J. (1966). Finding the shortest route between two points in a network. The Computer Journal, 9(3), 275–280.

Pohl, I. (1971). Bi-directional Search. Machine Intelligence, 6, 124–140.

Весь код в одном jupyter noteboook

Можно позапускать онлайн в биндере если не хочется локально настраивать питон

P. S.

Я не знаю, связано ли это с обширным использованием формул или чем-то ещё, но работа с формулами в WYSIWYG -- это адская боль. Уже недавно на хабре писали, что задача центрирования неразрешима, но есть еще более сложная задача: превращение -- в Надеюсь, что в будущем ситуация изменится.