Алгоритмы поиска максимального потока[TBD]

1. Краткая сводка теории

Теорема. Максимальный поток равен минимальному разрезу.

Величина разреза — сумма пропускных способностей только попутных дуг (из части с истоком в часть со стоком). Встречные дуги не считаем.
Когда максимальный поток построен, на границе минимального разреза все попутные дуги заполнены ( $f = c$ ), а все встречные пусты ( $f = 0$ ). Из $s$ через такую границу не пройти.

Как найти сам разрез. После окончания алгоритма запускаем ещё один DFS из $s$ по открытым дугам. Вершины, до которых он дошёл, образуют левую часть разреза $S$ , остальные вершины составляют $T$ . Это лишь один из минимальных разрезов, других может быть несколько.

Три алгоритма отличаются способом поиска увеличивающего пути:

Алгоритм	Как ищем путь	Асимптотика
Форд — Фалкерсон	DFS	$O (E \cdot F_{max})$
Эдмондс — Карп	BFS, то есть путь с наименьшим числом рёбер	$O (V \cdot E^{2})$
Диниц	DFS в слоистой сети, построенной BFS	$O (V^{2} E)$

Не стоит думать, что Эдмондс — Карп «лучше» Форда — Фалкерсона. Асимптотика у него просто другая: она не зависит от величины потока, но зависит от размеров графа. Какой алгоритм выбрать, решают ограничения задачи.

2. Как выбирать алгоритм: оценка по ограничениям

В учебной задаче: $V = 100$ , $E = 1000$ , пропускные способности до $10^{5}$ .

Поток не больше суммы пропускных способностей всех рёбер: $F_{max} \leq 10^{3} \cdot 10^{5} = 10^{8}$ .
Форд — Фалкерсон: $E \cdot F_{max} \leq 10^{11}$ . Оценка очень грубая и на пределе, но тесты слабые, поэтому пройдёт.
Эдмондс — Карп: $V \cdot E^{2} = 10^{2} \cdot 10^{6} = 10^{8}$ . Заметно лучше.
Диниц: $V^{2} E = 10^{4} \cdot 10^{3} = 10^{7}$ . Ещё лучше.

Предупреждение про вещественные пропускные способности: Форд — Фалкерсон хорошо работает именно для целых чисел. С вещественными его DFS может «улететь в бесконечность» (см. конспект прошлой лекции).

3. Хранение графа

Во всех задачах на потоки удобно хранить граф матрицей смежности. Это странно, но вершин обычно немного, а обращаться к паре $(v, u)$ приходится постоянно.

Нужны две матрицы $N \times N$ :

capacity[v][u] — пропускная способность дуги $v \to u$ (не меняется в ходе алгоритма);
flow[v][u] — текущий поток по дуге (меняется).

Остаточная пропускная способность вдоль дуги $v \to u$ вычисляется единым выражением:

d = capacity [v] [u] - flow [v] [u] .

Почему это работает и для движения против ребра

Когда мы пускаем поток $Δ$ по дуге $v \to u$ , делаем две вещи:

cpp

flow[v][u] += delta;
flow[u][v] -= delta;   // отрицательное число, это нормально

Теперь рассмотрим, что покажет формула $d = capacity - flow$ для разных ситуаций.

Пример. Сеть из 4 вершин (после перевода в нумерацию с нуля): дуги $0 \to 1$ (пропускная способность 1), $0 \to 2$ (2), $2 \to 1$ (1), $1 \to 3$ (2), $2 \to 3$ (1). Максимальный поток равен 3: минимальный разрез можно провести хоть вокруг истока, хоть вокруг стока, и величина будет $3$ , а, например, разрез с попутной дугой $2 \to 1$ даёт $4$ , он не минимальный.

Что проверяем	`capacity`	`flow`	$d = c - f$	Смысл
Дуга $0 \to 1$ после того, как пустили 1	$1$	$1$	$0$	дуга заполнена, закрыта
Дуга $0 \to 2$ , потока нет	$2$	$0$	$2$	можно протолкнуть 2
Обратный проход $1 \to 0$ (дуги нет)	$0$	$- 1$	$0 - (- 1) = 1$	можно вернуть 1 единицу потока
Пары, между которыми дуги нет вообще	$0$	$0$	$0$	дуга закрыта

Видно, что единственная формула покрывает все случаи. Значение $d = 0$ означает ровно одно из трёх:

попутная дуга полностью заполнена;
встречная дуга полностью пуста (по ней нечего «возвращать»);
такой дуги в графе просто нет.

Ещё деталь: capacity мы не трогаем, а flow меняем «хитро». При добавлении потока по $v \to u$ увеличиваем flow[v][u], но уменьшаем flow[u][v], чтобы обратная дуга оказалась открытой и по ней можно было пойти и уменьшить поток.

Чтение входных данных

cpp

int n, m;
cin >> n >> m;
vector<vector<int>> capacity(n, vector<int>(n));
vector<vector<int>> flow(n, vector<int>(n));

for (int i = 0; i < m; ++i) {
    int v, u, c;               // c — capacity, а не вес
    cin >> v >> u >> c;
    --v; --u;                  // в задаче нумерация с 1
    capacity[v][u] += c;       // += на случай кратных рёбер
}

4. Форд — Фалкерсон: реализация

Общая схема очень проста:

text

ответ = 0
пока true:
    сбросить used
    delta = dfs(исток, INF)
    если delta == 0: выйти
    ответ += delta

INF — «сколько мы вообще могли бы протолкнуть»: по условию пропускные способности не больше $10^{5}$ , поэтому берём, например, $10^{6}$ ;
delta тащим вдоль пути и берём минимум остаточных пропускных способностей, в итоге получаем узкое место пути.

Массив used сбрасывается перед каждым поиском пути. Делаем его vector<int>, а не vector<bool>: специализированный vector<bool> работает медленно, и это может испортить асимптотику.

cpp

int main() {
    // ... чтение, как выше ...
    const int INF = 1e6;
    vector<int> used(n);

    function<int(int, int)> dfs = [&](int v, int delta) -> int {
        if (v == n - 1) return delta;           // дошли до стока: ура
        used[v] = 1;
        for (int u = 0; u < n; ++u) {
            int d = capacity[v][u] - flow[v][u];
            if (used[u] || d <= 0) continue;    // были там или дуга закрыта
            int add = dfs(u, min(delta, d));    // запоминаем узкое место
            if (add > 0) {                      // путь найден, откатываемся
                flow[v][u] += add;
                flow[u][v] -= add;
                return add;
            }
        }
        return 0;                               // увеличивающего пути здесь нет
    };

    int answer = 0;
    while (true) {
        fill(used.begin(), used.end(), 0);
        int delta = dfs(0, INF);
        if (delta == 0) break;
        answer += delta;
    }
    cout << answer << "\n";
}

Как это работает:

зашли в вершину — пометили used;
перебираем всех соседей $u$ по матрице смежности;
вычисляем $d = c - f$ ; если вершина уже посещена или $d = 0$ (дуга закрыта), пропускаем;
рекурсивно идём в $u$ , передавая $min (delta, d)$ ;
если вернулось положительное число, значит, путь найден и узкое место известно. На откате по рекурсии обновляем поток на каждой дуге пути и возвращаем это же число наверх.

Лямбда с захватом [&] здесь выбрана только ради краткости, чтобы не таскать used, capacity, flow параметрами и не загромождать суть алгоритма. «По жизни» так писать не стоит: функции лучше передавать данные явно.

5. Эдмондс — Карп: реализация

Внешний цикл тот же самый, но вместо DFS запускаем BFS, который ищет кратчайший по числу рёбер путь. Затем по списку предков восстанавливаем путь и пускаем по нему поток.

Как устроен BFS:

массив dist заполнен «бесконечностью». Бесконечностью можно взять n: BFS строит дерево обхода, в нём путь не длиннее $n - 1$ ребра, поэтому $n$ недостижимо;
массив parent (список предков) заполнен значением -1; в нём храним, из какой вершины пришли;
dist[0] = 0, кладём исток в очередь;
берём вершину из очереди и просматриваем соседей. Если $c - f = 0$ (дуга закрыта или её нет), пропускаем. Иначе обычная релаксация: если dist[v] + 1 < dist[u], то обновляем dist[u], записываем parent[u] = v и кладём $u$ в очередь.

Здесь мы смотрим только на то, открыта ли дуга, а не на величину пропускной способности: нужен просто путь с минимальным числом рёбер. Вершины естественным образом раскладываются по слоям (расстоянию от истока).

cpp

const int INF_CAP = 1e6;
vector<int> dist, parent;

auto bfs = [&]() -> bool {
    dist.assign(n, n);          // «бесконечность»
    parent.assign(n, -1);       // тип int, а не size_t!
    queue<int> q;
    dist[0] = 0;
    q.push(0);
    while (!q.empty()) {
        int v = q.front(); q.pop();
        for (int u = 0; u < n; ++u) {
            if (capacity[v][u] - flow[v][u] <= 0) continue;
            if (dist[v] + 1 < dist[u]) {
                dist[u] = dist[v] + 1;
                parent[u] = v;
                q.push(u);
            }
        }
    }
    return parent[n - 1] != -1;     // дошли ли до стока
};

int answer = 0;
while (bfs()) {
    // 1-й проход по предкам: считаем узкое место
    int delta = INF_CAP;
    for (int u = n - 1; u != 0; u = parent[u]) {
        int v = parent[u];
        delta = min(delta, capacity[v][u] - flow[v][u]);
    }
    // 2-й проход: обновляем поток
    for (int u = n - 1; u != 0; u = parent[u]) {
        int v = parent[u];
        flow[v][u] += delta;
        flow[u][v] -= delta;
    }
    answer += delta;
}

Замечания:

Путь в parent восстанавливается от стока к истоку: рёбра идут вперёд, а список предков указывает назад.
Приходится пройти по списку предков два раза: сначала узнать минимум, потом изменить поток (пока минимум не известен, потоки менять нельзя). Это можно написать лучше, например, считать минимум по пути прямо в BFS отдельным массивом, но здесь пишем «как обсуждали».
Для DFS-версии минимум считался «по дороге туда», для BFS приходится восстанавливать путь «по дороге обратно».

6. Алгоритм Диница

Идея: два алгоритма в один

Диниц — это «смесь BFS и DFS»:

слоистую сеть строим так же, как в Эдмондсе — Карпе (BFS);
по слоистой сети проталкиваем поток как в Форде — Фалкерсоне (DFS), но так, чтобы никогда не ходить дважды одним и тем же тупиковым путём.

Слоистая сеть

Запускаем BFS из истока и получаем для каждой вершины расстояние dist (число рёбер от истока). Вершины раскладываются по слоям: слой 0, слой 1, слой 2, …

В слоистую сеть входят только рёбра, которые ведут из слоя $k$ в слой $k + 1$ (и только открытые, то есть с $c - f > 0$ ). Остальные рёбра (внутри слоя и назад) для DFS неинтересны, их игнорируем.

Блокирующий поток

Эдмондс — Карп после построения слоистой сети толкает поток вдоль одного пути и строит BFS заново. Это расточительно. Вместо этого:

строим слоистую сеть;
находим в ней увеличивающий путь DFS-ом и проталкиваем поток;
после этого хотя бы одно ребро пути становится заполненным (узкое горлышко) и выпадает из сети;
запускаем DFS снова и снова, пока он находит путь в той же слоистой сети;
когда путь в слоистой сети найти нельзя («поток заблокирован»), строим сеть заново.

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

Указатель `ptr` для амортизации

Как запомнить, что ребро выкинуто? Для каждой вершины $v$ храним указатель ptr[v]: индекс соседа, с которого нужно продолжать перебор.

Индекс только растёт (монотонно). Если из $v$ по соседу ничего не нашли (или ребро заполнилось), сдвигаем ptr[v] дальше и никогда не возвращаемся.
Так мы никогда не заходим повторно туда, где уже были и ничего не нашли.

Оценка сложности

Каждый DFS идёт вдоль слоистой сети, а путь в ней не длиннее $V - 1$ ребра. За счёт ptr проталкивание всех путей в одной слоистой сети обходится дёшево, суммарно порядка $V \cdot E$ .
Слоистых сетей строить придётся не более $V$ штук: после блокировки кратчайший путь до стока в остаточной сети строго удлиняется (рёбра, по которым пустили поток, становятся обратными и открывают новые пути, поэтому следующая сеть «глубже»), а длина пути не превосходит $V$ .
Один BFS стоит $O (V + E)$ .

Итого — квадратичная по $V$ оценка $O (V^{2} E)$ .

Код

cpp

vector<int> dist, ptr;

// BFS строит слоистую сеть; список предков больше не нужен,
// важно лишь расстояние от истока
auto bfs = [&]() -> bool {
    dist.assign(n, n);
    queue<int> q;
    dist[0] = 0;
    q.push(0);
    while (!q.empty()) {
        int v = q.front(); q.pop();
        for (int u = 0; u < n; ++u) {
            if (capacity[v][u] - flow[v][u] <= 0) continue;
            if (dist[v] + 1 < dist[u]) {
                dist[u] = dist[v] + 1;
                q.push(u);
            }
        }
    }
    return dist[n - 1] != n;       // достижим ли сток
};

// DFS толкает поток вдоль слоистой сети с запоминанием ptr
function<int(int, int)> dfs = [&](int v, int pushed) -> int {
    if (v == n - 1) return pushed;
    for (int& u = ptr[v]; u < n; ++u) {            // продолжаем с того места, где остановились
        if (dist[u] != dist[v] + 1) continue;      // не слоистая сеть
        int d = capacity[v][u] - flow[v][u];
        if (d <= 0) continue;                      // дуга закрыта
        int add = dfs(u, min(pushed, d));
        if (add > 0) {
            flow[v][u] += add;
            flow[u][v] -= add;
            return add;
        }
    }
    return 0;
};

int answer = 0;
while (bfs()) {                                    // пока строится слоистая сеть
    fill(ptr.begin(), ptr.end(), 0);
    while (int pushed = dfs(0, INF)) answer += pushed;   // блокирующий поток
}

Что нужно заметить:

Внутренняя часть DFS — это ровно Форд — Фалкерсон, с той разницей, что ходим только по рёбрам слоистой сети и используем ptr[v] вместо used.
Внешняя часть — Эдмондс — Карп (BFS и цикл «пока сток достижим»).
Пока можно построить слоистую сеть, максимальный поток не достигнут; как только сеть построить нельзя, поток максимален.
Диниц пишется легко и работает быстро. Проще всего его объяснять как «написали Эдмондса — Карпа, написали Форда — Фалкерсона, а потом объединили».

Алгоритмы поиска максимального потока[TBD] ​

1. Краткая сводка теории ​

2. Как выбирать алгоритм: оценка по ограничениям ​

3. Хранение графа ​

Почему это работает и для движения против ребра ​

Чтение входных данных ​

4. Форд — Фалкерсон: реализация ​

5. Эдмондс — Карп: реализация ​

6. Алгоритм Диница ​

Идея: два алгоритма в один ​

Слоистая сеть ​

Блокирующий поток ​

Указатель ptr для амортизации ​

Оценка сложности ​

Код ​