Max flow, minimum cut

在微软苏州已经入职了一个月了，生活开始渐渐稳定下来。和在北京的各种不堪相比，我要感谢这次调职，得以让我从 easy 模式开始社畜生活。同时搁置好久的算法又可以得以继续研究下去了。

今天讲的算法是最大流最小割算法。实际上这个是我在《算法导论》上看来看去都看不明白的一个章节：术语太多，而且代码很少（我好菜啊）。幸运的是，hihoCoder 上面的一系列专题为我垫了不少砖头，让我得以掌握这一块知识。

对于一个最大流问题，其定义一般都是这样的：有一个有向图 G，源点为 s，汇点为 t，每条边都有其对应的流最大容量 c，源点 s 只有发出的流量，汇点 t 只有汇入的流量，求能从源点到汇点的最大可能流量。

那么，我先不想说证明过程，我直接给出结论：我们用一种反复查找的方法，试图在当前图 G 上找到一条从 s 到 t 的路径，其边上允许的流量非零。每次找到这一条路径之后，也就可以确定这条路径上的流量瓶颈，将路径上的边的可行流量减去这一流量瓶颈，在新图上进行下一次查找。我们一直持续这样的查找，直到无法从 s 到 t 走出一条流量瓶颈大于 0 的路径。这个算法叫做 Ford-Fulkerson 算法。这个查找方式可以使用 BFS 进行。

嘛，基本上就是这样。但是由于图中可能会存在反向边，对于我上面讲的这个情况，对于下面这种情况是不成立的：

An tricky example of flow network

如果我们用 BFS 找到了 A-B-D-F 这条可行路径，其瓶颈是 3，减去该瓶颈之后无法找到更新的路径。然而其最大流是 5：由 B 分出的流量有一份到 D，另外两份到 E，然后 A-C-D 有一个流量为 2 的路径。这时候，我们需要一个“残余网络”的概念帮助我们解决这个问题：每条边维护一条对应的反向边，其大小为当前在正向上面已经消耗掉的流量，而正向边的容量为当前正向剩余的可行流量。换句话说，正向边和反向边的存在，帮助我们维护了每条边还有多少剩余流量可用。如果有一条流量容量为 1 的边，其正向流量为 1000000，反向流量为 999999，这也是可行的！我们在残余网络上找到的一条可行路径，叫做增广路径。于是我上面描述的算法可以如下表达：

1
2
3
4
While ( findAugmentPath() ) // 判断是否有增广路
	maxFlow = maxFlow + delta // 最大流增加
	modifyGraph() // 对增广路进行修改
End While

与最小割的关系

我们通常说 max flow, minimun cut。那么什么被定义为 cut？在一个带有 s 和 t 的网络流里，有一种划分把点划分到两个不相交的 S 和 T 集合中，$s \in S, t \in T$. 净流 f(S, T) 被定义为穿过割 (S, T) 的流量之和（当一个割经过反向边的时候，反向边上的流量应为负值）。割的容量 C(S, T) 被定义为这条割上所有的容量之和（不包括反向边）。也就是说，$f(S, T) \le C(S, T)$。 可以证明，当前网络的流量总是等于任意一个割的净流。但是，任意的割不会有相同的容量。其中最小的那个割的容量对应的割，我们称之为最小割。

现在有一个结论：对于任意网络流图来说，最大流一定等于最小割的容量。这个结论是最大流最小割定理的直接推论，这个定理由三个等价的表达组成：

  1. f 是网络的最大流；
  2. 该网络的残余网络不包含任何增广路径。
  3. 流网络的某个切割 (S, T) 的容量等于 f。

因此，求出了最大流，也就求出了最小割的最大容量。使用上文所提及的 Ford-Fulkerson 算法，能够快速算出割、网络流量。

Production Code

使用 BFS 进行每次增广路径查找的 Ford-Fulkerson 实现叫做 Edmonds-Karp 算法。由于 BFS 时间复杂度（约）为 O(E)，流量递增操作的总次数为 O(VE)，该算法的时间复杂度为 $O(VE^2)$。这并不是最快的算法。现在有一种算法，通过记录每个点到汇点 t 的最短距离维持搜索顺序，使用 DFS 的方法进行增广路径的查找，能大大降低时间复杂度。这一算法的代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
// 本代码对应于 http://hihocoder.com/problemset/problem/1369
#include <cstdio>
#include <cstring>
#include <queue>

const int inf = 0x3f3f3f3f;
const int maxn = 505, maxm = 40005;
int N, M;

typedef struct {
	int u, v, c, prev;
}edge_t;
int head[maxn];
edge_t edges[maxm];

int source, sink;
int dep[maxn], depcnt[maxn], cur[maxn], stack[maxn];

void addedge(int u, int v, int c) {
	static int id = 0;
	edges[id].u = u; edges[id].v = v; edges[id].c = c; edges[id].prev = head[u]; head[u] = id++;
	edges[id].u = v; edges[id].v = u; edges[id].c = 0; edges[id].prev = head[v]; head[v] = id++;
}

void bfs() {
	memset(dep, 0xff, sizeof(dep));
	memset(depcnt, 0, sizeof(depcnt));
	dep[sink] = 0;
	depcnt[dep[sink]] = 1;
	std::queue<int> Q;
	Q.push(sink);
	int u, v;
	while (!Q.empty()) {
		u = Q.front();
		Q.pop();
		for (int e = head[u]; ~e; e = edges[e].prev) if (edges[e].c == 0 && dep[edges[e].v] == -1) {
			v = edges[e].v;
			Q.push(v);
			dep[v] = dep[u] + 1;
			depcnt[dep[v]]++;
		}
	}
}

int SAP() {
	bfs();
	memcpy(cur, head, sizeof(head));
	int u = source, maxflow = 0, top = 0, i;
	while (dep[source]<N) {
		if (u == sink) {
			//find the bottleneck
			int delta = inf, p;
			for (int t = top - 1; t >= 0; t--) {
				int e = stack[t];
				if (edges[e].c<delta) {
					delta = edges[e].c;
					p = t;
				}
			}
			for (int t = top - 1; t >= 0; t--) {
				int e = stack[t];
				edges[e].c -= delta;
				edges[e ^ 1].c += delta;
			}
			maxflow += delta;
			top = p;
			u = edges[stack[p]].u;
		}
		for (i = cur[u]; ~i; i = edges[i].prev) {
			if (edges[i].c>0 && dep[edges[i].v] == dep[u] - 1) {
				cur[u] = i;
				u = edges[i].v;
				stack[top++] = i;
				break;
			}
		}
		if (i == -1) {
			depcnt[dep[u]]--;
			if (depcnt[dep[u]] == 0) break;
			int mindep = N;
			for (int e = head[u]; ~e; e = edges[e].prev) {
				if (edges[e].c>0 && mindep > dep[edges[e].v]) {
					mindep = dep[edges[e].v];
					cur[u] = e;
				}
			}
			dep[u] = mindep + 1;
			depcnt[dep[u]]++;
			if (u != source) {
				u = edges[stack[--top]].u;
			}
		}
	}
	return maxflow;
}

int main() {
	freopen("testcase", "r", stdin);
	scanf("%d%d", &N, &M);
	source = 1, sink = N;
	memset(head, 0xff, sizeof(head));
	for (int i = 1; i <= M; i++) {
		int u, v, c;
		scanf("%d%d%d", &u, &v, &c);
		addedge(u, v, c);
	}

	printf("%d\n", SAP());
	return 0;
}

嘛，这个优化后的算法叫什么我也不知道，我只是从别人的解法里学习过来的。但是这个代码做了至少这么几个工作：

  1. 在进行 DFS 寻找增广路径之前用 BFS 从汇点搜索建立 dep 和 depcnt，即到汇点的最短距离，这样为在 DFS 的时候总是走最短路径提供了依据。
  2. 正向边和反向边成对存储，索引为偶数的总是正向边；正向边和反向边可以通过索引的异或切换；
  3. stack 里面存储通过 DFS 寻找的增广路径。
  4. 在每一次 DFS 寻找过后，回溯到瓶颈路径的起点重新搜索，节省时间；
  5. 维护一个 cur，其在物理意义上与 head 相同，也是为了直接从当前的可行边搜索，节省时间；
  6. 但是如果从某点找不到可行边，需要更新该点的 dep 与 depcnt。这也就意味着 cur 也要更新。

就是这样！但是有的文献指出，第一步建立 dep 与 depcnt 可以不进行。这主要的实现参考 http://hihocoder.com/contest/hiho117/solution/898175。

主要参考文献：

  1. hihocoder hiho一下 115、116、117、118、119 周
  2. 《算法导论》（第三版）第 26 章
  3. Max Flow-SAP-Improved Shortest Augmenting

留言

▲