網路流（2）：最大流
概述
我們有一張圖，要求從源點流向匯點的最大流量（可以有很多條路到達匯點），就是我們的最大流問題。
Ford-Fulkerson 增廣路演算法
該方法通過尋找增廣路來更新最大流，有 EK,dinic,SAP,ISAP 主流演算法。
求解最大流之前，我們先認識一些概念。
殘量網路
首先我們介紹一下一條邊的剩餘容量 \(c_f(u,v)\) （Residual Capacity），它表示的是這條邊的容量與流量之差，即 \(c_f(u,v)=c(u,v)-f(u,v)\) 。
對於流函式 \(f\) ，殘存網路 \(G_f\) （Residual Network）是網路 \(G\) 中所有結點 和剩餘容量大於 0
注意，剩餘容量大於 0 的邊可能不在原圖 \(G\) 中（根據容量、剩餘容量的定義以及流函式的斜對稱性得到）。可以理解為，殘量網路中包括了那些還剩了流量空間的邊構成的圖，也包括虛邊（即反向邊）。
増廣路
在原圖 \(G\) 中若一條從源點到匯點的路徑上所有邊的 剩餘容量都大於 0 ，這條路被稱為增廣路（Augmenting Path）。
或者說，在殘存網路 \(G_f\) 中，一條從源點到匯點的路徑被稱為增廣路。如圖：
我們從 \(4\) 到 \(3\)
-
\(4\rightarrow2\rightarrow3\) 這條 增廣路 的總流量為 \(20\) 。到 \(2\) 的時候還是 \(30\) ，到 \(3\) 了就只有 \(20\) 了。
-
\(4\rightarrow2\rightarrow1\rightarrow3\) 這樣子我們就很好的保留了 \(30\) 的流量。
所以我們這張圖的最大流就應該是 \(20+30=50\) 。
求最大流是很簡單的，接下來講解求最大流的 3 種方法。
Edmond-Karp 動能演算法（EK 演算法）
這個演算法很簡單，就是 BFS 找增廣路 ，然後對其進行 增廣 。你可能會問，怎麼找？怎麼增廣？
-
找？我們就從源點一直 BFS 走來走去，碰到匯點就停，然後增廣（每一條路都要增廣）。我們在 BFS 的時候就注意一下流量合不合法就可以了。
-
增廣？其實就是按照我們找的增廣路在重新走一遍。走的時候把這條路的能夠成的最大流量減一減，然後給答案加上最小流量就可以了。
再講一下 反向邊 。增廣的時候要注意建造反向邊，原因是這條路不一定是最優的，這樣子程式可以進行反悔。假如我們對這條路進行增廣了，那麼其中的每一條邊的反向邊的流量就是它的流量。
講一下一些小細節。如果你是用鄰接矩陣的話，反向邊直接就是從 \(table[x,y]\) 變成 \(table[y,x]\) 。如果是常用的鏈式前向星，那麼在加入邊的時候就要先加入反向邊。那麼在用的時候呢，我們直接 \(i\operatorname{xor}1\) 就可以了 ( \(i\) 為邊的編號）。為什麼呢？相信大家都是知道 \(\operatorname{xor}\) 的，那麼我們在加入正向邊後加入反向邊，就是靠近的，所以可以使用 \(\operatorname{xor}\) 。我們還要注意一開始的編號要設定為 \(tot=1\) ，因為邊要從編號 \(2\) 開始，這樣子 \(\operatorname{xor}\) 對編號 \(2,3\) 的邊才有效果。
EK 演算法的時間複雜度為 \(O(nm^2)\) （其中 \(n\) 為點數， \(m\) 為邊數）。效率還有很大提升空間。
#define maxn 250
#define INF 0x3f3f3f3f
struct Edge {
int from, to, cap, flow;
Edge(int u, int v, int c, int f) : from(u), to(v), cap(c), flow(f) {}
};
struct EK {
int n, m;
vector<Edge> edges;
vector<int> G[maxn];
int a[maxn], p[maxn];
void init(int n) {
for (int i = 0; i < n; i++) G[i].clear();
edges.clear();
}
void AddEdge(int from, int to, int cap) {
edges.push_back(Edge(from, to, cap, 0));
edges.push_back(Edge(to, from, 0, 0));
m = edges.size();
G[from].push_back(m - 2);
G[to].push_back(m - 1);
}
int Maxflow(int s, int t) {
int flow = 0;
for (;;) {
memset(a, 0, sizeof(a));
queue<int> Q;
Q.push(s);
a[s] = INF;
while (!Q.empty()) {
int x = Q.front();
Q.pop();
for (int i = 0; i < G[x].size(); i++) {
Edge& e = edges[G[x][i]];
if (!a[e.to] && e.cap > e.flow) {
p[e.to] = G[x][i];
a[e.to] = min(a[x], e.cap - e.flow);
Q.push(e.to);
}
}
if (a[t]) break;
}
if (!a[t]) break;
for (int u = t; u != s; u = edges[p[u]].from) {
edges[p[u]].flow += a[t];
edges[p[u] ^ 1].flow -= a[t];
}
flow += a[t];
}
return flow;
}
};
Dinic 演算法
Dinic 演算法 的過程是這樣的：每次增廣前，我們先用 BFS 來將圖分層。設源點的層數為 \(0\) ，那麼一個點的層數便是它離源點的最近距離。
通過分層，我們可以幹兩件事情：
- 如果不存在到匯點的增廣路（即匯點的層數不存在），我們即可停止增廣。
- 確保我們找到的增廣路是最短的。（原因見下文）
接下來是 DFS 找增廣路的過程。
我們每次找增廣路的時候，都只找比當前點層數多 \(1\) 的點進行增廣（這樣就可以確保我們找到的增廣路是最短的）。
Dinic 演算法有兩個優化：
- 多路增廣 ：每次找到一條增廣路的時候，如果殘餘流量沒有用完怎麼辦呢？我們可以利用殘餘部分流量，再找出一條增廣路。這樣就可以在一次 DFS 中找出多條增廣路，大大提高了演算法的效率。
- 當前弧優化 ：如果一條邊已經被增廣過，那麼它就沒有可能被增廣第二次。那麼，我們下一次進行增廣的時候，就可以不必再走那些已經被增廣過的邊。
設點數為 \(n\) ，邊數為 \(m\) ，那麼 Dinic 演算法的時間複雜度是 \(O(n^{2}m)\) ，在稀疏圖上效率和 EK 演算法相當，但在稠密圖上效率要比 EK 演算法高很多。
特別地，在求解二分圖最大匹配問題時，可以證明 Dinic 演算法的時間複雜度是 \(O(m\sqrt{n})\) 。
#define maxn 250
#define INF 0x3f3f3f3f
struct Edge {
int from, to, cap, flow;
Edge(int u, int v, int c, int f) : from(u), to(v), cap(c), flow(f) {}
};
struct Dinic {
int n, m, s, t;
vector<Edge> edges;
vector<int> G[maxn];
int d[maxn], cur[maxn];
bool vis[maxn];
void init(int n) {
for (int i = 0; i < n; i++) G[i].clear();
edges.clear();
}
void AddEdge(int from, int to, int cap) {
edges.push_back(Edge(from, to, cap, 0));
edges.push_back(Edge(to, from, 0, 0));
m = edges.size();
G[from].push_back(m - 2);
G[to].push_back(m - 1);
}
bool BFS() {
memset(vis, 0, sizeof(vis));
queue<int> Q;
Q.push(s);
d[s] = 0;
vis[s] = 1;
while (!Q.empty()) {
int x = Q.front();
Q.pop();
for (int i = 0; i < G[x].size(); i++) {
Edge& e = edges[G[x][i]];
if (!vis[e.to] && e.cap > e.flow) {
vis[e.to] = 1;
d[e.to] = d[x] + 1;
Q.push(e.to);
}
}
}
return vis[t];
}
int DFS(int x, int a) {
if (x == t || a == 0) return a;
int flow = 0, f;
for (int& i = cur[x]; i < G[x].size(); i++) {
Edge& e = edges[G[x][i]];
if (d[x] + 1 == d[e.to] && (f = DFS(e.to, min(a, e.cap - e.flow))) > 0) {
e.flow += f;
edges[G[x][i] ^ 1].flow -= f;
flow += f;
a -= f;
if (a == 0) break;
}
}
return flow;
}
int Maxflow(int s, int t) {
this->s = s;
this->t = t;
int flow = 0;
while (BFS()) {
memset(cur, 0, sizeof(cur));
flow += DFS(s, INF);
}
return flow;
}
};
ISAP
這個是 SAP 演算法的加強版 (Improved)。
struct Edge {
int from, to, cap, flow;
Edge(int u, int v, int c, int f) : from(u), to(v), cap(c), flow(f) {}
};
bool operator<(const Edge& a, const Edge& b) {
return a.from < b.from || (a.from == b.from && a.to < b.to);
}
struct ISAP {
int n, m, s, t;
vector<Edge> edges;
vector<int> G[maxn];
bool vis[maxn];
int d[maxn];
int cur[maxn];
int p[maxn];
int num[maxn];
void AddEdge(int from, int to, int cap) {
edges.push_back(Edge(from, to, cap, 0));
edges.push_back(Edge(to, from, 0, 0));
m = edges.size();
G[from].push_back(m - 2);
G[to].push_back(m - 1);
}
bool BFS() {
memset(vis, 0, sizeof(vis));
queue<int> Q;
Q.push(t);
vis[t] = 1;
d[t] = 0;
while (!Q.empty()) {
int x = Q.front();
Q.pop();
for (int i = 0; i < G[x].size(); i++) {
Edge& e = edges[G[x][i] ^ 1];
if (!vis[e.from] && e.cap > e.flow) {
vis[e.from] = 1;
d[e.from] = d[x] + 1;
Q.push(e.from);
}
}
}
return vis[s];
}
void init(int n) {
this->n = n;
for (int i = 0; i < n; i++) G[i].clear();
edges.clear();
}
int Augment() {
int x = t, a = INF;
while (x != s) {
Edge& e = edges[p[x]];
a = min(a, e.cap - e.flow);
x = edges[p[x]].from;
}
x = t;
while (x != s) {
edges[p[x]].flow += a;
edges[p[x] ^ 1].flow -= a;
x = edges[p[x]].from;
}
return a;
}
int Maxflow(int s, int t) {
this->s = s;
this->t = t;
int flow = 0;
BFS();
memset(num, 0, sizeof(num));
for (int i = 0; i < n; i++) num[d[i]]++;
int x = s;
memset(cur, 0, sizeof(cur));
while (d[s] < n) {
if (x == t) {
flow += Augment();
x = s;
}
int ok = 0;
for (int i = cur[x]; i < G[x].size(); i++) {
Edge& e = edges[G[x][i]];
if (e.cap > e.flow && d[x] == d[e.to] + 1) {
ok = 1;
p[e.to] = G[x][i];
cur[x] = i;
x = e.to;
break;
}
}
if (!ok) {
int m = n - 1;
for (int i = 0; i < G[x].size(); i++) {
Edge& e = edges[G[x][i]];
if (e.cap > e.flow) m = min(m, d[e.to]);
}
if (--num[d[x]] == 0) break;
num[d[x] = m + 1]++;
cur[x] = 0;
if (x != s) x = edges[p[x]].from;
}
}
return flow;
}
};
Push-Relabel 預流推進演算法
該方法在求解過程中忽略流守恆性，並每次對一個結點更新資訊，以求解最大流。
通用的預流推進演算法
首先我們介紹預流推進演算法的主要思想，以及一個可行的暴力實現演算法。
預流推進演算法通過對單個結點的更新操作，直到沒有結點需要更新來求解最大流。
演算法過程維護的流函式不一定保持流守恆性，對於一個結點，我們允許進入結點的流超過流出結點的流，超過的部分被稱為結點 \(u(u\in V-\{s,t\})\) 的 超額流 \(e(u)\) ：
\[e(u)=\sum_{(x,u)\in E}f(x,u)-\sum_{(u,y)\in E}f(u,y) \]
若 \(e(u)>0\) ，稱結點 \(u\) 溢位 。
預流推進演算法維護每個結點的高度 \(h(u)\) ，並且規定溢位的結點 \(u\) 如果要推送超額流，只能向高度小於 \(u\) 的結點推送；如果 \(u\) 沒有相鄰的高度小於 \(u\) 的結點，就修改 \(u\) 的高度（重貼標籤）。
高度函式
準確地說，預流推進維護以下的一個對映 \(h:V\to \mathbf{N}\) ：
- \(h(s)=|V|,h(t)=0\)
- \(\forall (u,v)\in E_f,h(u)\leq h(v)+1\)
稱 \(h\) 是殘存網路 \(G_f=(V_f,E_f)\) 的高度函式。
引理 1：設 \(G_f\) 上的高度函式為 \(h\) ，對於任意兩個結點 \(u,v\in V\) ，如果 \(h(u)>h(v)+1\) ，則 \((u,v)\) 不是 \(G_f\) 中的邊。
演算法只會在 \(h(u)=h(v)+1\) 的邊執行推送。
推送（Push）
適用條件：結點 \(u\) 溢位，且存在結點 \(v((u,v)\in E_f,c(u,v)-f(u,v)>0,h(u)=h(v)+1)\) ，則 push 操作適用於 \((u,v)\) 。
於是，我們儘可能將超額流從 \(u\) 推送到 \(v\) ，推送過程中我們只關心超額流和 \(c(u,v)-f(u,v)\) 的最小值，不關心 \(v\) 是否溢位。
如果 \((u,v)\) 在推送完之後滿流，將其從殘存網路中刪除。
重貼標籤（Relabel）
適用條件：如果結點 \(u\) 溢位，且 \(\forall (u,v)\in E_f,h(u)\leq h(v)\) ，則 relabel 操作適用於 \(u\) 。
則將 \(h(u)\) 更新為 \(min_{(u,v)\in E_f}h(v)+1\) 即可。
初始化
\[\begin{split} &\forall (u,v)\in E,&f(u,v)=\left\{\begin{split} &c(u,v)&,u=s\\ &0&,u\neq s\\ \end{split}\right. \\ &\forall u\in V,&h(u)=\left\{\begin{split} &|V|&,u=s\\ &0&,u\neq s\\ \end{split}\right. ,e(u)=\sum_{(x,u)\in E}f(x,u)-\sum_{(u,y)\in E}f(u,y) \end{split} \]
上述將 \((s,v)\in E\) 充滿流，並將 \(h(s)\) 抬高，使得 \((s,v)\notin E_f\) ，因為 \(h(s)>h(v)\) ，而且 \((s,v)\) 畢竟滿流，沒必要留在殘存網路中；上述還將 \(e(s)\) 初始化為 \(\sum_{(s,v)\in E}f(s,v)\) 的相反數。
通用演算法
我們每次掃描整個圖，只要存在結點 \(u\) 滿足 push 或 relabel 操作的條件，就執行對應的操作。
如圖，每個結點中間表示編號，左下表示高度值 \(h(u)\) ，右下表示超額流 \(e(u)\) ，結點顏色的深度也表示結點的高度；邊權表示 \(c(u,v)-f(u,v)\) ，綠色的邊表示滿足 \(h(u)=h(v)+1\) 的邊 \((u,v)\) （即殘存網路的邊 \(E_f\) ）：
整個演算法我們大致瀏覽一下過程，這裡筆者使用的是一個暴力演算法，即暴力掃描是否有溢位的結點，有就更新
最後的結果
可以發現，最後的超額流一部分回到了 \(s\) ，且除了源點匯點，其他結點都沒有溢位；這時的流函式 \(f\) 滿足流守恆性，為最大流，即 \(e(t)\) 。
核心程式碼：
const int N = 1e4 + 4, M = 1e5 + 5, INF = 0x3f3f3f3f;
int n, m, s, t, maxflow, tot;
int ht[N], ex[N];
void init() { // 初始化
for (int i = h[s]; i; i = e[i].nex) {
const int &v = e[i].t;
ex[v] = e[i].v, ex[s] -= ex[v], e[i ^ 1].v = e[i].v, e[i].v = 0;
}
ht[s] = n;
}
bool push(int ed) {
const int &u = e[ed ^ 1].t, &v = e[ed].t;
int flow = min(ex[u], e[ed].v);
ex[u] -= flow, ex[v] += flow, e[ed].v -= flow, e[ed ^ 1].v += flow;
return ex[u]; // 如果 u 仍溢位，返回 1
}
void relabel(int u) {
ht[u] = INF;
for (int i = h[u]; i; i = e[i].nex)
if (e[i].v) ht[u] = min(ht[u], ht[e[i].t]);
++ht[u];
}
HLPP 演算法
最高標號預流推進演算法（High Level Preflow Push）是基於預流推進演算法的優先佇列實現，該演算法優先推送高度高的溢位的結點，演算法演算法複雜度 \(O(n^2\sqrt m)\) 。
具體地說，HLPP 演算法過程如下：
- 初始化（基於預流推進演算法）；
- 選擇溢位結點（除 \(s,t\) ）中高度最高的結點 \(u\) ，並對它所有可以推送的邊進行推送；
- 如果 \(u\) 仍溢位，對它重貼標籤，回到步驟 2；
- 如果沒有溢位的結點，演算法結束。
BFS 優化
HLPP 的上界為 \(O(n^2\sqrt m)\) ，但在使用時卡得比較緊；我們可以在初始化高度的時候進行優化。具體來說，我們初始化 \(h(u)\) 為 \(u\) 到 \(t\) 的最短距離；特別地， \(h(s)=n\) 。
在 BFS 的同時我們順便檢查圖的連通性，排除無解的情況。
GAP 優化
HLPP 推送的條件是 \(h(u)=h(v)+1\) ，而如果在演算法的某一時刻， \(h(u)=t\) 的結點個數為 \(0\) ，那麼對於 \(h(u)>t\) 的結點就永遠無法推送超額流到 \(t\) ，因此只能送回 \(s\) ，那麼我們就在這時直接讓他們的高度變成 \(n+1\) ，以儘快推送回 \(s\) ，減少重貼標籤的操作。
LuoguP4722【模板】最大流 加強版/預流推進
#include <cstdio>
#include <cstring>
#include <queue>
using namespace std;
const int N = 1e4 + 4, M = 2e5 + 5, INF = 0x3f3f3f3f;
int n, m, s, t;
struct qxx {
int nex, t, v;
};
qxx e[M * 2];
int h[N], cnt = 1;
void add_path(int f, int t, int v) { e[++cnt] = (qxx){h[f], t, v}, h[f] = cnt; }
void add_flow(int f, int t, int v) {
add_path(f, t, v);
add_path(t, f, 0);
}
int ht[N], ex[N], gap[N]; // 高度；超額流；gap 優化
bool bfs_init() {
memset(ht, 0x3f, sizeof(ht));
queue<int> q;
q.push(t), ht[t] = 0;
while (q.size()) { // 反向 BFS, 遇到沒有訪問過的結點就入隊
int u = q.front();
q.pop();
for (int i = h[u]; i; i = e[i].nex) {
const int &v = e[i].t;
if (e[i ^ 1].v && ht[v] > ht[u] + 1) ht[v] = ht[u] + 1, q.push(v);
}
}
return ht[s] != INF; // 如果圖不連通，返回 0
}
struct cmp {
bool operator()(int a, int b) const { return ht[a] < ht[b]; }
}; // 偽裝排序函式
priority_queue<int, vector<int>, cmp> pq; // 將需要推送的結點以高度高的優先
bool vis[N]; // 是否在優先佇列中
int push(int u) { // 儘可能通過能夠推送的邊推送超額流
for (int i = h[u]; i; i = e[i].nex) {
const int &v = e[i].t, &w = e[i].v;
if (!w || ht[u] != ht[v] + 1) continue;
int k = min(w, ex[u]); // 取到剩餘容量和超額流的最小值
ex[u] -= k, ex[v] += k, e[i].v -= k, e[i ^ 1].v += k; // push
if (v != s && v != t && !vis[v])
pq.push(v), vis[v] = 1; // 推送之後，v 必然溢位，則入堆，等待被推送
if (!ex[u]) return 0; // 如果已經推送完就返回
}
return 1;
}
void relabel(int u) { // 重貼標籤（高度）
ht[u] = INF;
for (int i = h[u]; i; i = e[i].nex)
if (e[i].v) ht[u] = min(ht[u], ht[e[i].t]);
++ht[u];
}
int hlpp() { // 返回最大流
if (!bfs_init()) return 0; // 圖不連通
ht[s] = n;
memset(gap, 0, sizeof(gap));
for (int i = 1; i <= n; i++)
if (ht[i] != INF) gap[ht[i]]++; // 初始化 gap
for (int i = h[s]; i; i = e[i].nex) {
const int v = e[i].t, w = e[i].v; // 佇列初始化
if (!w) continue;
ex[s] -= w, ex[v] += w, e[i].v -= w, e[i ^ 1].v += w; // 注意取消 w 的引用
if (v != s && v != t && !vis[v]) pq.push(v), vis[v] = 1; // 入隊
}
while (pq.size()) {
int u = pq.top();
pq.pop(), vis[u] = 0;
while (push(u)) { // 仍然溢位
// 如果 u 結點原來所在的高度沒有結點了，相當於出現斷層
if (!--gap[ht[u]])
for (int i = 1; i <= n; i++)
if (i != s && i != t && ht[i] > ht[u] && ht[i] < n + 1) ht[i] = n + 1;
relabel(u);
++gap[ht[u]]; // 新的高度，更新 gap
}
}
return ex[t];
}
int main() {
scanf("%d%d%d%d", &n, &m, &s, &t);
for (int i = 1, u, v, w; i <= m; i++) {
scanf("%d%d%d", &u, &v, &w);
add_flow(u, v, w);
}
printf("%d", hlpp());
return 0;
}
感受一下執行過程