通過 USACO 4.2.1 Ditch 學習一下最大流算法 。可惜它給的測試數據幾乎沒有任何殺傷力，后面測試時我們采用 DD_engi 寫的程序生成的加強版數據。
總體上來說，最大流算法分為兩大類：增廣路 (Augmenting Path) 和預流推進重標號 (Push Relabel) 。也有算法同時借鑒了兩者的長處，如 Improved SAP 。本篇主要介紹增廣路類算法，思想、復雜度及實際運行效率比較，並試圖從中選擇一種兼顧代碼復雜度和運行效率的較好方案。以下我們將會看到，有時理論分析的時間復雜度並不能很好的反映一種算法的實際效率。
1. Ford - Fulkerson 方法
所有增廣路算法的基礎都是 Ford - Fulkerson 方法。稱之為方法而不是算法是因為 Ford - Fulkerson 只提供了一類思想，在此之上的具體操作可有不同的實現方案。
給定一個有向網絡 G(V,E) 以及源點 s 終點 t ，FF 方法描述如下：
Ford-Fulkerson 方法 (G,s,t) 1 將各邊上流量 f 初始化為 0 2 while 存在一條增廣路徑 p 3 do 沿路徑 p 增廣流量 f 4 return f
假設有向網絡 G 中邊 (i,j) 的容量為 c(i,j) ，當前流量為 f(i,j) ，則此邊的剩余流量即為 r(i,j) = c(i,j) - f(i,j) ，其反向邊的剩余流量為 r(j,i) = f(i,j) 。有向網中所有剩余流量 r(i,j) > 0 的邊構成殘量網絡 Gf ，增廣路徑p即是殘量網絡中從源點 s 到終點 t 的路徑。
沿路徑 p 增廣流量 f 的操作基本都是相同的，各算法的區別就在於尋找增廣路徑 p 的方法不同。例如可以尋找從 s 到 t 的最短路徑，或者流量最大的路徑。
2. Edmonds - Karp 算法
Shortest Augmenting Path (SAP) 是每次尋找最短增廣路的一類算法，Edmonds - Karp 算法以及后來著名的 Dinic 算法都屬於此。SAP 類算法可統一描述如下：
Shortest Augmenting Path 1 x <-- 0 2 while 在殘量網絡 Gx 中存在增廣路 s ~> t 3 do 找一條最短的增廣路徑 P 4 delta <-- min{rij:(i,j) 屬於 P} 5 沿 P 增廣 delta 大小的流量 6 更新殘量網絡 Gx 7 return x
在無權邊的有向圖中尋找最短路，最簡單的方法就是廣度優先搜索 (BFS)，E-K 算法就直接來源於此。每次用一遍 BFS 尋找從源點 s 到終點 t 的最短路作為增廣路徑，然后增廣流量 f 並修改殘量網絡，直到不存在新的增廣路徑。
E-K 算法的時間復雜度為 O(VE2)，由於 BFS 要搜索全部小於最短距離的分支路徑之后才能找到終點，因此可以想象頻繁的 BFS 效率是比較低的。實踐中此算法使用的機會較少。
3. Dinic 算法
BFS 尋找終點太慢，而 DFS 又不能保證找到最短路徑。1970年 Dinic 提出一種思想，結合了 BFS 與 DFS 的優勢，采用構造分層網絡的方法可以較快找到最短增廣路，此算法又稱為阻塞流算法 (Blocking Flow Algorithm)。
首先定義分層網絡 AN(f)。在殘量網絡中從源點 s 起始進行 BFS，這樣每個頂點在 BFS 樹中會得到一個距源點 s 的距離 d，如 d(s) = 0，直接從 s 出發可到達的點距離為 1，下一層距離為2 ... 。稱所有具有相同距離的頂點位於同一層，在分層網絡中，只保留滿足條件 d(i) + 1 = d(j) 的邊，這樣在分層網絡中的任意路徑就成為到達此頂點的最短路徑。
Dinic 算法每次用一遍 BFS 構建分層網絡 AN(f)，然后在 AN(f) 中一遍 DFS 找到所有到終點 t 的路徑增廣；之后重新構造 AN(f)，若終點 t 不在 AN(f) 中則算法結束。DFS 部分算法可描述如下：
1 p <-- s 2 while s 的出度 > 0 do 3 u <-- p.top 4 if u != t then 5 if u 的出度 > 0 then 6 設 (u,v) 為 AN(f) 中一條邊 7 p <-- p, v 8 else 9 從 p 和 AN(f) 中刪除點 u 以及和 u 連接的所有邊 10 else 11 沿 p 增廣 12 令 p.top 為從 s 沿 p 可到達的最后頂點 13 end while
實際代碼中不必真的用一個圖來存儲分層網絡，只需保存每個頂點的距離標號並在 DFS 時判斷 dist[i] + 1 = dist[j] 即可。Dinic 的時間復雜度為 O(V2E)。由於較少的代碼量和不錯的運行效率，Dinic 在實踐中比較常用。具體代碼可參考 DD_engi 網絡流算法評測包中的標程，這幾天 dinic 算法的實現一共看過和比較過將近 10 個版本了，DD 寫的那個在效率上是數一數二的，邏輯上也比較清晰。
4. Improved SAP 算法
本次介紹的重頭戲。通常的 SAP 類算法在尋找增廣路時總要先進行 BFS，BFS 的最壞情況下復雜度為 O(E)，這樣使得普通 SAP 類算法最壞情況下時間復雜度達到了 O(VE2)。為了避免這種情況，Ahuja 和 Orlin 在1987年提出了Improved SAP 算法，它充分利用了距離標號的作用，每次發現頂點無出弧時不是像 Dinic 算法那樣到最后進行 BFS，而是就地對頂點距離重標號，這樣相當於在遍歷的同時順便構建了新的分層網絡，每輪尋找之間不必再插入全圖的 BFS 操作，極大提高了運行效率。國內一般把這個算法稱為 SAP...顯然這是不准確的，畢竟從字面意思上來看 E-K 和 Dinic 都屬於 SAP，我還是習慣稱為 ISAP 或改進的 SAP 算法。
與 Dinic 算法不同，ISAP 中的距離標號是每個頂點到達終點 t 的距離。同樣也不需顯式構造分層網絡，只要保存每個頂點的距離標號即可。程序開始時用一個反向 BFS 初始化所有頂點的距離標號，之后從源點開始，進行如下三種操作：(1)當前頂點 i 為終點時增廣 (2) 當前頂點有滿足 dist[i] = dist[j] + 1 的出弧時前進 (3) 當前頂點無滿足條件的出弧時重標號並回退一步。整個循環當源點 s 的距離標號 dist[s] >= n 時結束。對 i 點的重標號操作可概括為 dist[i] = 1 + min{dist[j] : (i,j)屬於殘量網絡Gf}。具體算法描述如下：
algorithm Improved-Shortest-Augmenting-Path 1 f <-- 0 2 從終點 t 開始進行一遍反向 BFS 求得所有頂點的起始距離標號 d(i) 3 i <-- s 4 while d(s) < n do 5 if i = t then // 找到增廣路 6 Augment 7 i <-- s // 從源點 s 開始下次尋找 8 if Gf 包含從 i 出發的一條允許弧 (i,j) 9 then Advance(i) 10 else Retreat(i) // 沒有從 i 出發的允許弧則回退 11 return f procedure Advance(i) 1 設 (i,j) 為從 i 出發的一條允許弧 2 pi(j) <-- i // 保存一條反向路徑，為回退時准備 3 i <-- j // 前進一步，使 j 成為當前結點 procedure Retreat(i) 1 d(i) <-- 1 + min{d(j):(i,j)屬於殘量網絡Gf} // 稱為重標號的操作 2 if i != s 3 then i <-- pi(i) // 回退一步 procedure Augment 1 pi 中記錄為當前找到的增廣路 P 2 delta <-- min{rij:(i,j)屬於P} 3 沿路徑 P 增廣 delta 的流量 4 更新殘量網絡 Gf
算法中的允許弧是指在殘量網絡中滿足 dist[i] = dist[j] + 1 的弧。Retreat 過程中若從 i 出發沒有弧屬於殘量網絡 Gf 則把頂點距離重標號為 n 。
雖然 ISAP 算法時間復雜度與 Dinic 相同都是 O(V2E)，但在實際表現中要好得多。要提的一點是關於 ISAP 的一個所謂 GAP 優化。由於從 s 到 t 的一條最短路徑的頂點距離標號單調遞減，且相鄰頂點標號差嚴格等於1，因此可以預見如果在當前網絡中距離標號為 k (0 <= k < n) 的頂點數為 0，那么可以知道一定不存在一條從 s 到 t 的增廣路徑，此時可直接跳出主循環。在我的實測中，這個優化是絕對不能少的，一方面可以提高速度，另外可增強 ISAP 算法時間上的穩定性，不然某些情況下 ISAP 會出奇的費時，而且大大慢於 Dinic 算法。相對的，初始的一遍 BFS 卻是可有可無，因為 ISAP 可在循環中自動建立起分層網絡。實測加不加 BFS 運行時間差只有 5% 左右，代碼量可節省 15~20 行。
5. 最大容量路徑算法 (Maximum Capacity Path Algorithm)
1972年還是那個 E-K 組合提出的另一種最大流算法。每次尋找增廣路徑時不找最短路徑，而找容量最大的。可以預見，此方法與 SAP 類算法相比可更快逼近最大流，從而降低增廣操作的次數。實際算法也很簡單，只用把前面 E-K 算法的 BFS 部分替換為一個類 Dijkstra 算法即可。USACO 4.2 節的說明詳細介紹了此算法，這里就不詳述了。
時間復雜度方面。BFS 是 O(E)，簡單 Dijkstra 是 O(V2)，因此效果可想而知。但提到 Dijkstra 就不能不提那個 Heap 優化，雖然 USACO 的算法例子中沒有用 Heap ，我自己還是實現了一個加 Heap 的版本，畢竟 STL 的優先隊列太好用了不加白不加啊。效果也是非常明顯的，但比起 Dinic 或 ISAP 仍然存在海量差距，這里就不再詳細介紹了。
6. Capacity Scaling Algorithm
不知道怎么翻比較好，索性就這么放着吧。叫什么的都有，容量縮放算法、容量變尺度算法等，反正就那個意思。類似於二分查找的思想，尋找增廣路時不必非要局限於尋找最大容量，而是找到一個可接受的較大值即可，一方面有效降低尋找增廣路時的復雜度，另一方面增廣操作次數也不會增加太多。時間復雜度 O(E2logU) 實際效率嘛大約稍好於最前面 BFS 的 E-K 算法，稀疏圖時表現較優，但仍然不敵 Dinic 與 ISAP。
7. 算法效率實測！
重頭戲之二，雖然引用比較多，哎~
首先引用此篇強文 《Maximum Flow: Augmenting Path Algorithms Comparison》
對以上算法在稀疏圖、中等稠密圖及稠密圖上分別進行了對比測試。直接看結果吧：
稀疏圖：
ISAP 輕松拿下第一的位置，圖中最左邊的 SAP 應該指的是 E-K 算法，這里沒有比較 Dinic 算法是個小遺憾吧，他把 Dinic 與 SAP 歸為一類了。最大流量路徑的簡單 Dijkstra 實現實在是太失敗了 - -，好在 Heap 優化后還比較能接受……可以看到 Scaling 算法也有不錯的表現。
中等稠密圖：
ISAP 依然領先。最大流量算法依然不太好過……幾個 Scaling 類算法仍然可接受。
稠密圖：
ISAP……你無敵了！這次可以看出 BFS 的 Scaling 比 DFS 實現好得多，而且幾乎與 Improved Scaling 不相上下，代碼量卻差不少。看來除 ISAP 外 BFS Scaling 也是個不錯的選擇。
最后補個我自己實測的圖，比較算法有很多是 DD 網絡流算法評測包里的標程，評測系統用的 Cena，評測數據為 DD ditch 數據生成程序生成的加強版數據：
我一共寫了 7 個版本的 ISAP ，Dinic 算法也寫了幾個遞歸版的但效率都不高，只放上來一個算了。從上圖來看似乎 ISAP 全面超越了號稱最大流最快速算法的 HLPP，但在另外一台機器上測試結果與此卻不大相同，有時 ISAP 與 HLPP 基本持平，有時又稍稍慢一些。在這種差距非常小的情況下似乎編譯器的效果也比較明顯。那個 HLPP 是用 PASCAL 寫的，我不知在 Win32 平台下編譯代碼效率如何，至少我的幾個 ISAP 用 VC2008 + SP1 編譯比用 g++ 要強不少，也可能是參數設置的問題。
不過這些都是小事，關鍵是見證了 ISAP 的實際效率。從上面可以看出不加 GAP 優化的 ISAP 有幾個測試點干脆無法通過，而不加 BFS 卻無甚大影響，遞歸與非遞歸相差在 7% 左右的樣子。綜合以上表現，推薦采用 ISAP 不加 BFS，非遞歸 + GAP 優化的寫法，Ditch 這道題一共也才 80 行左右代碼。要提的一點是 GAP 優化用遞歸來表現的話不如 while 循環來得直接。另外，看起來 HLPP 表現確實很優秀，有機會也好好研究一下吧，預流推進重標號算法也是一大類呢……
最后附上一個 ISAP + GAP + BFS 的非遞歸版本代碼，網絡采用鄰接表 + 反向弧指針：
1 #include<cstdio> 2 #include<memory> 3 using namespace std; 4 5 const int maxnode = 1024; 6 const int infinity = 2100000000; 7 8 struct edge{ 9 int ver; // vertex 10 int cap; // capacity 11 int flow; // current flow in this arc 12 edge *next; // next arc 13 edge *rev; // reverse arc 14 edge(){} 15 edge(int Vertex, int Capacity, edge *Next) 16 :ver(Vertex), cap(Capacity), flow(0), next(Next), rev((edge*)NULL){} 17 void* operator new(size_t, void *p){ 18 return p; 19 } 20 }*Net[maxnode]; 21 int dist[maxnode]= {0}, numbs[maxnode] = {0}, src, des, n; 22 23 void rev_BFS(){ 24 int Q[maxnode], head = 0, tail = 0; 25 for(int i=1; i<=n; ++i){ 26 dist[i] = maxnode; 27 numbs[i] = 0; 28 } 29 30 Q[tail++] = des; 31 dist[des] = 0; 32 numbs[0] = 1; 33 while(head != tail){ 34 int v = Q[head++]; 35 for(edge *e = Net[v]; e; e = e->next){ 36 if(e->rev->cap == 0 || dist[e->ver] < maxnode)continue; 37 dist[e->ver] = dist[v] + 1; 38 ++numbs[dist[e->ver]]; 39 Q[tail++] = e->ver; 40 } 41 } 42 } 43 44 int maxflow(){ 45 int u, totalflow = 0; 46 edge *CurEdge[maxnode], *revpath[maxnode]; 47 for(int i=1; i<=n; ++i)CurEdge[i] = Net[i]; 48 u = src; 49 while(dist[src] < n){ 50 if(u == des){ // find an augmenting path 51 int augflow = infinity; 52 for(int i = src; i != des; i = CurEdge[i]->ver) 53 augflow = min(augflow, CurEdge[i]->cap); 54 for(int i = src; i != des; i = CurEdge[i]->ver){ 55 CurEdge[i]->cap -= augflow; 56 CurEdge[i]->rev->cap += augflow; 57 CurEdge[i]->flow += augflow; 58 CurEdge[i]->rev->flow -= augflow; 59 } 60 totalflow += augflow; 61 u = src; 62 } 63 64 edge *e; 65 for(e = CurEdge[u]; e; e = e->next) 66 if(e->cap > 0 && dist[u] == dist[e->ver] + 1)break; 67 if(e){ // find an admissible arc, then Advance 68 CurEdge[u] = e; 69 revpath[e->ver] = e->rev; 70 u = e->ver; 71 } else { // no admissible arc, then relabel this vertex 72 if(0 == (--numbs[dist[u]]))break; // GAP cut, Important! 73 CurEdge[u] = Net[u]; 74 int mindist = n; 75 for(edge *te = Net[u]; te; te = te->next) 76 if(te->cap > 0)mindist = min(mindist, dist[te->ver]); 77 dist[u] = mindist + 1; 78 ++numbs[dist[u]]; 79 if(u != src) 80 u = revpath[u]->ver; // Backtrack 81 } 82 } 83 return totalflow; 84 } 85 86 int main(){ 87 int m, u, v, w; 88 freopen("ditch.in", "r", stdin); 89 freopen("ditch.out", "w", stdout); 90 while(scanf("%d%d", &m, &n) != EOF){ // POJ 1273 need this while loop 91 edge *buffer = new edge[2*m]; 92 edge *data = buffer; 93 src = 1; des = n; 94 while(m--){ 95 scanf("%d%d%d", &u, &v, &w); 96 Net[u] = new((void*) data++) edge(v, w, Net[u]); 97 Net[v] = new((void*) data++) edge(u, 0, Net[v]); 98 Net[u]->rev = Net[v]; 99 Net[v]->rev = Net[u]; 100 } 101 rev_BFS(); 102 printf("%d\n", maxflow()); 103 delete [] buffer; 104 } 105 return 0; 106 }