Как и большинство других алгоритмов нахождения максимального потока, алгоритм Форда-Фалкерсона работает на остаточной сети, построенной на основе исходного графа. Сам алгоритм заключается в нахождении раз за разом (увеличивающего) пути от источника к стоку и проталкивании максимального потока по этому пути с ограничением на пропускную способность ребер. Как только путь от источника к стоку перестает существовать, алгоритм завершается.
Запишем все вышесказанное в виде кода, где неуточненные места вынесены в функции, с которыми мы будем разбираться позднее. Итак, наше решение заключено в функцию max_flow, которая принимает две вершины s (источник) и t (сток). Очевидно, что если источник и сток совпадают, то алгоритм можно не запускать. Сам алгоритм работает на уже построенной остаточной сети и описан в функции ford_fulkerson. В ней мы в цикле пытаемся построить какой-нибудь путь от вершины s. Если при этом вершина t даже не была посещена, то пути не существует, и алгоритм завершится. Если же вершина t была посещена, то путь от s к t существует, и мы его восстанавливаем. Затем мы находим максимальную величину потока, который еще можно протолкнуть по найденному пути, и увеличиваем поток на эту величину по всем ребрам пути.
def ford_fulkerson(s, t):
result_flow = 0
while True:
for v in range(n):
visited[v] = False
pred[v] = None
find_path(s)
if not visited[t]:
break
path = restore_path(t)
flow = path_capacity(path)
push_path(path, flow)
result_flow += flow
return result_flow
def max_flow(s, t):
if s == t:
return None
build_network()
return ford_fulkerson(s, t)При эффективной реализации на списках смежности и с целочисленными пропускными способностями алгоритм может достигать времени работы O((n + m)f), где f -- величина максимального потока. Каждый найденный путь увеличивает поток хотя бы на 1, а поиск пути может быть выполнен за O(m). Обнуление массивов pred и visited выполняется за O(n) на каждой итерации.
Некоторые части алгоритма мы вынесли в отдельные функции. Теперь пора их реализовать.
Начнем с поиска пути. Воспользуемся для этого поиском в глубину, при этом будем запоминать ребро, по которому мы пришли в вершину, в массиве pred и помечать посещенные вершины в массиве visited. Одним отличием от обычного поиска в глубину является то, что мы не можем перейти по ребру, если оно заблокировано потоком (по ребру пущен поток величиной равной пропускной способности ребра).
Так как мы работаем не с исходным графом, а с остаточной сетью, то введем дополнительные функции для работы с ней: edges (получение исходящих ребер в остаточной сети), target (получение конечной вершины дуги), available (вычисление доступной пропускной способности).
def find_path(v):
visited[v] = True
for e in edges(v):
u = target(e)
if not visited[u] and available(e) > 0:
pred[u] = e
find_path(u)Следующим этапом является восстановление пути по массиву pred. Начиная со стока t, будем по ребрам перемещаться все ближе к источнику s. В этой цепочке вершин s -- единственная вершина, у которой нет ребра, по которому в нее пришли, так как с нее был начат обход. Именно по этому критерию мы завершим восстановление пути.
Введем еще одну функцию для работы с остаточной сетью source, которая возвращает начальную вершину дуги.
def restore_path(v):
path = []
while pred[v] is not None:
e = pred[v]
path.append(e)
v = source(e)
return pathМаксимальная величина потока, который можно протолкнуть по найденному пути, равна наименьшей остаточной пропускной способности среди ребер этого пути, поэтому пройдем по всем ребрам пути и найдем минимум. Воспользуемся фактом, что путь содержит хотя бы одно ребро (следует из условия, что s != t), взяв за базовое значение остаточную пропускную способность первого ребра в пути.
def path_capacity(path):
cap = available(path[0])
for e in path:
cap = min(cap, available(e))
return capЭта функция может быть также упрощена до одной строчки:
def path_capacity(path):
return min(map(available, path))Проталкивание потока вдоль пути заключается увеличении потока вдоль каждого из ребер пути на указанную величину. Для этого введем функцию push, которая увеличивает поток вдоль одного ребра.
def push_path(path, flow):
for e in path:
push(e, flow)Остались не описаны следующие функции:
build_networkдля построения остаточной сети на основе исходного графаsourceдля получения начальной вершины ребра в остаточной сетиtargetдля получения конечной вершины ребра в остаточной сетиavailableдля получения остаточной пропускной способности ребраflowдля получения величины потока, пропущенного по ребруedgesдля получения исходящих из вершины ребер в остаточной сетиpushдля увеличения потока вдоль ребра в остаточной сети
Существуют разные подходы к предоставлению подобного интерфейса, и основная проблема при этом заключается в том, что в остаточной сети все ребра дублируются и изменение одного ребра влечет за собой изменение и второго ребра. При этом находить второе ребро нужно быстро и желательно за константное время. Мы разберем один из подходов при организации остаточной сети на списках смежности.
Для начала определим структуру ребра. Ребро характеризуется начальной и конечной вершинами, пропускной способностью и потоком, пропущенным по этому ребру.
@dataclass
class Edge:
source: int
target: int
capacity: int
flow: intКаждому ребру (v, cu, u) в графе ставится в соответствие два ребра в остаточной сети: (v, cu, u) и (u, 0, v). При этом изначально по ним пропущен поток величины 0. Будем хранить ребра остаточной сети в виде списка, причем каждая пара ребер будет занимать позиции 2i и 2i + 1. Таким образом, по четности индекса легко определить индекс обратного ребра. Более того, эту операцию легко выразить через ^ 1 (xor 1), так как (2i) ^ 1 = 2i + 1 и, соответственно, (2i + 1) ^ 1 = 2i.
Саму остаточную сеть будем, как и исходный граф, представлять в виде списков смежности, но с одним отличием: вместо пары (пропускная способность, вершина) будем хранить индекс ребра в списке ребер.
network = [[] for v in range(n)]
flow_edges = []
def build_network():
for v in range(n):
for cu, u in g[v]:
network[v].append(len(flow_edges))
flow_edges.append(Edge(source=v, target=u, capacity=cu, flow=0))
network[u].append(len(flow_edges))
flow_edges.append(Edge(source=u, target=v, capacity=0, flow=0))Таким образом, достаточно легко реализовать работающие за константу недостающие функции.
def edges(v):
return network[v]
def available(e):
edge = flow_edges[e]
return edge.capacity - edge.flow
def flow(e):
edge = flow_edges[e]
return edge.flow
def target(e):
edge = flow_edges[e]
return edge.target
def source(e):
edge = flow_edges[e]
return edge.sourceПри увеличении потока по ребру в остаточной сети необходимо уменьшить поток в обратном ребре на ту же величину.
def push(e, flow):
edge = flow_edges[e]
edge.flow += flow
edge = flow_edges[e ^ 1]
edge.flow -= flowИсходя из кода можно заметить, что для любого ребра (прямого или обратного) всегда выполняется соотношение flow <= capacity. По прямым рёбрам течёт неотрицательный поток, по обратным -- неположительный. Сумма потока по прямому и обратному ребру всегда равна нулю.
Предыдущий алгоритм не ограничивал выбор алгоритма для поиска увеличивающего пути. Алгоритм Эдмондса-Карпа, в свою очередь, предлагает использовать поиск в ширину, чтобы увеличивающий путь был всегда кратчайшим. Таким образом, необходимо изменить лишь одну функцию:
def find_path(s):
q = queue()
visited[s] = True
q.enqueue(s)
while not q.empty():
v = q.dequeue()
for e in edges(v):
u = target(e)
if not visited[u] and available(e) > 0:
visited[u] = True
pred[u] = e
q.enqueue(u)Сам каркас алгоритма при этом останется неизменным.
def edmonds_karp(s, t): # <-- здесь!
result_flow = 0
while True:
for v in range(n):
visited[v] = False
pred[v] = None
find_path(s)
if not visited[t]:
break
path = restore_path(t)
flow = path_capacity(path)
push_path(path, flow)
result_flow += flow
return result_flow
def max_flow(s, t):
if s == t:
return None
build_network()
return edmonds_karp(s, t) # <-- здесь!Каждая итерация алгоритма выполняется за O(n + m). Каждый следующий найденный данным алгоритмом увеличивающий путь будет не короче предыдущего. При нахождении увеличивающего пути хотя бы одно ребро блокируется потоком и может стать вновь доступно лишь при проталкивании потока по обратному ребру, что может случиться лишь при более длинном увеличивающем пути. Так как увеличивающий путь кратчайший, его длина не превосходит n - 1. Таким образом, число итераций составляет O(nm). Получаем время работы алгоритма O((n+m)nm).
В этой задаче ребрам графа ставится в соответствие не только пропускная способность, но и удельная стоимость потока (неотрицательная). Соответственно, наша цель найти такой поток среди максимальных потоков в этом графе, что его стоимость минимальна. При этом стоимость потока рассчитывается как сумма по всем ребрам произведений удельной стоимости потока ребра на величину потока в ребре (учитываются лишь ребра исходного графа).
Добавим удельную стоимость к представлению ребра, а также определим новую функцию cost к нашей абстракции, которая возвращает удельную стоимость для переданного ребра:
@dataclass
class Edge:
source: int
target: int
capacity: int
cost: int # <-- здесь!
flow: int
def cost(e):
edge = flow_edges[e]
return edge.costВ связи с этим достроим нашу остаточную сеть, чтобы она поддерживала также удельные стоимости, при этом удельная стоимость для обратного ребра будет равна удельной стоимости усходного ребра со знаком минуса. Будем считать, что списки смежности g содержат тройки (пропускная способность, удельная стоимость, вершина).
def build_network():
for v in range(n):
for cu, pu, u in g[v]: # <-- здесь!
network[v].append(len(flow_edges))
flow_edges.append(Edge(source=v, target=u, capacity=cu, cost=pu, flow=0)) # <-- здесь!
network[u].append(len(flow_edges))
flow_edges.append(Edge(source=u, target=v, capacity=0, cost=-pu, flow=0)) # <-- и здесь!Существует несколько алгоритмов нахождения максимального потока минимальной стоимости. Мы рассмотрим два из них.
Первый подход заключается в изначальном построении потока минимальной стоимости, с постепенным увеличением его величины, пока не получим максимальный поток. При этом на каждой итерации мы получаем поток минимальной стоимости для данной величины. Этот алгоритм называется алгоритмом Басакера-Гоуэна, и его тоже можно рассматривать как некоторую модификацию алгоритма Форда-Фалкерсона. Отличие опять же заключается в поиске увеличивающего пути. Алгоритм предлагает искать кратчайшие по удельной стоимости пути в остаточной сети.
Приведем сначала каркас алгоритма. Заметим, что теперь мы считаем не только величину потока, но и его стоимость. Также мы ввели новую функцию path_cost, которая вычисляет удельную стоимость потока для найденного пути.
def busacker_gowen(s, t):
result_flow = 0
result_cost = 0 # <-- здесь!
while True:
for v in range(n):
pred[v] = None
dist[v] = None # <-- заменили `visited` на `dist`
find_path(s)
if dist[t] is None:
break
path = restore_path(t)
flow = path_capacity(path)
pcost = path_cost(path) # <-- здесь!
push_path(path, flow)
result_flow += flow
result_cost += flow * pcost # <-- еще здесь!
return result_flow, result_cost # <-- и здесь!
def min_cost_max_flow(s, t):
if s == t:
return None
build_network()
return busacker_gowen(s, t) # <-- и здесь!Функция path_cost достаточно проста и очевидна, так что начнем с нее. Удельной стоимостью потока для пути является сумма удельных стоимостей ребер, входящих в этот путь.
def path_cost(path):
return sum(map(cost, path))Теперь вернемся к функции find_path. Как уже говорилось, мы должны огранизовать поиск кратчайшего по удельной стоимости пути. Так как обратные ребра имеют отрицательную стоимость, то мы можем использовать алгоритм Форда-Беллмана. Так как на каждом этапе мы находим поток наименьшей стоимости, то отрицательных циклов в остаточной сети не имеется, и алгоритм отработает корректно. Также нельзя забывать про то, что ребро в сети должно обрабатываться только при условии, если его остаточная пропускная способность положительна.
def find_path(s):
dist[s] = 0
for i in range(n - 1): # <-- не делаем последнюю итерацию
for e in network_edges():
v = source(e)
u = target(e)
c = cost(e)
if dist[v] is None or available(e) == 0:
continue
if dist[u] is None or dist[u] > dist[v] + c:
dist[u] = dist[v] + c
pred[u] = eТак как Форд-Беллман работает со списком ребер, то добавим функцию его получения к абстракции:
def network_edges():
return range(len(flow_edges))Еще один подход заключается в построении любого максимального потока с дальнейшим уменьшением его стоимости. Первая часть выполняется любым из алгоритмов нахождения максимального потока, например алгоритмом Эдмондса-Карпа. Во второй части мы добавляем в сеть удельные стоимости. Если поток не обладает минимальной стоимостью, то в остаточной сети существуют отрицательные по удельной стоимости циклы. Соответственно, алгоритм заключается в нахождении и устранении отрицательных циклов в сети.
Сначала представим каркас нашего алгоритма:
def min_cost_max_flow(s, t):
if s == t:
return None
build_network()
result_flow = edmonds_karp(s, t)
remove_cycles()
result_cost = flow_cost()
return result_flow, result_costНачнем с простой функции flow_cost, которая посчитает стоимость потока в сети. Напомним, что она равна сумме произведений удельных стоимостей ребер на велечину потока в них. При этом считать необходимо только ребра из исходного графа. Но при этом, эта же величина для обратных ребер будет одинакова за счет зеркальности величины потока и удельной стоимости, поэтому можно посчитать значение для всех ребер остаточной сети и разделить на два.
def flow_cost():
edge_cost = lambda e: cost(e) * flow(e)
return sum(map(edge_cost, network_edges())) // 2Теперь перейдем к функции remove_cycles. В ней мы будем в цикле искать отрицательные циклы, пропускать по ним поток, устраняя их и тем самым уменьшая стоимость потока. Как только отрицательный цикл не найден, алгоритм можно завершать. Также нельзя начинать поиск цикла из вершины s, так как в остаточной сети, по которой уже пропущен максимальный поток, некоторые вершины и, соответственно, циклы могут быть недостижимы из вершины s. Поэтому начнем алгоритм "одновременно" из всех вершин сети, присвоив им значение 0 в массиве dist (при этом физический смысл значений в массиве теряется).
def remove_cycles():
while True:
for v in range(n):
pred[v] = None
dist[v] = 0 # <-- здесь!
v = find_cycle()
if v is None:
break
cycle = restore_cycle(v)
flow = path_capacity(cycle)
push_path(cycle, flow)Находить цикл отрицательной стоимости будем с помощью алгоритма Форда-Беллмана. Если на последней итерации алгоритма была найдена вершина, для которой проведена релаксация, то эта вершина достижима из искомого цикла, и именно ее мы вернем.
def find_cycle():
for i in range(n): # <-- делаем последнюю итерацию
for e in network_edges():
v = source(e)
u = target(e)
c = cost(e)
if available(e) == 0:
continue
if dist[u] > dist[v] + c:
dist[u] = dist[v] + c
pred[u] = e
if i + 1 == n:
return u
return NoneТеперь, когда у нас есть вершина u, достижимая из цикла, мы можем добраться до вершины v, принадлежащей циклу, следуя по массиву pred. Так как на n-й итерации алгоритма Форда-Беллмана мы построили кратчайшие пути длины не более n, то достаточно пройти назад по ребрам n раз, чтобы гарантированно попасть в цикл. Затем, начиная от вершины v, мы можем восстановить все ребра цикла, следуя по массиву pred, пока вновь не встретим вершину v.
def restore_cycle(start):
for _ in range(n):
start = source(pred[start])
cycle = []
v = start
while True:
e = pred[v]
cycle.append(e)
v = source(e)
if v == start:
break
return cycle