天天看點

最大流Dinic算法講解 && ISAP 算法解釋短路增值算法(MPLA)Dinic算法 網絡流-最大流問題 ISAP 算法解釋

轉自:http://blog.csdn.net/wall_f/article/details/8207595

為了更好的介紹Dinic算法,我們先來介紹最短增廣路算法。

最短增廣路算法

1、頂點的層次和層次網絡

頂點的層次:在殘留網絡中,把從源點到頂點u的最短路徑長度(該長度僅僅是值路徑上邊的數目,與容量無關),稱為頂點u的層次,記為level(u)。源點Vs的層次為0。

将殘留網絡中所有的頂點的層次标注出來的過程稱為分層。

注意:

(1)對殘留網路進行分層後,弧可能有3種可能的情況。

1、從第i層頂點指向第i+1層頂點。

2、從第i層頂點指向第i層頂點。

3、從第i層頂點指向第j層頂點(j < i)。

(2)不存在從第i層頂點指向第i+k層頂點的弧(k>=2)。

(3)并非所有的網絡都能分層。

層次網絡:對殘留網絡進行分層後,删去比彙點Vt層次更高的頂點和與彙點Vt同層的頂點(保留Vt),并删去這些頂點相關聯的弧,再删去從某層頂點指向同層頂點和低層頂點的弧,所剩餘的各條弧的容量與殘留網絡中的容量相同,這樣得到的網絡就是殘留網絡的子網絡,稱為層次網絡,記為G''(V'',E'')。

根據層次網絡定義,層次網絡中任意的一條弧<u,v>,有滿足level(u)+1 == level(v),這條弧也叫允許弧。直覺的說,層次網絡是建立在殘留網絡基礎之上的一張“最短路徑圖”。從源點開始,在層次網絡中沿着邊不管怎麼走,到達一個終點之後,經過的路徑一定是終點在殘留網絡中的最短路徑。

阻塞流:設容量網絡的一個可行流為f,當該網絡的層次G''中不存在增廣路(即從源點Vs到彙點Vt的路徑時),稱該可行流f為層次網絡G''的阻塞流。

2、最短路增廣路徑的算法思想

最短增廣路的算法思想是:每次在層次網絡中找一條含弧數最少的增廣路進行增廣。最短增廣路算法的具體步驟如下:

(1)初始化容量網絡和網絡流。

(2)構造殘留網絡和層次網絡,若彙點不在層次網絡中,則算法結束。

(3)在層次網絡中不斷用BFS增廣,直到層次網絡中沒有增廣路為止;每次增廣完畢,在層次網絡中要去掉因改進流量而導緻飽和的弧。

(4)轉步驟(2)。

在最短增廣路算法中,第(2)、(3)步被循環執行,将執行(2)、(3)步的一次循環稱為一個階段。在每個階段中,首先根據殘留網絡建立層次網絡,然後不斷用BFS在層次網絡中增廣,直到出現阻塞流。注意每次增廣後,在層次網絡中要去掉因改進流量而飽和的弧。該階段的增廣完畢後,進入下一階段。這樣的不斷重複,直到彙點不在層次網絡中為止。彙點不在層次網絡中意味着在殘留網絡中不存在一條從源點到彙點的路徑,即沒有增廣路。

在程式實作的時候,并不需要真正“構造”層次網絡,隻需要對每個頂點标記層次,增廣的時候,判斷邊是否滿足level(v) = level(u)+1這一限制條件即可。

3、最短增廣路算法複雜度分析

最短增廣路的複雜度包括建立層次網絡和尋找增廣路兩部分。

在最短增廣路中,最多建立n個層次網絡,每個層次網絡用BFS一次周遊即可得到。一次BFS的複雜度為O(m),是以建層次圖的總複雜度為O(n*m)。

每增廣一次,層次網絡中必定有一條邊會被删除。層次網絡中最多有m條邊,是以認為最多可以增廣m次。在最短增廣路算法中,用BFS來增廣,一次增廣的複雜度為O(n+m),其中O(m)為BFS的花費,O(n)為修改流量的花費。是以在每一階段尋找增廣路的複雜度為O(m*(m+n)) = O(m*m)。是以n個階段尋找增廣路的算法總複雜度為O(n*m*m)。

兩者之和為O(n*m*m)。

以上介紹最短增廣路算法隻是為下面介紹Dinic算法而提供給大家一個鋪墊,有了以上的基礎,接下來我們來介紹Dinic算法,Dinic其實是最短增廣路算法的優化版。

連續最短增廣路算法----Dinic算法

1、Dinic算法思路

Dinic算法的思想也是分階段地在層次網絡中增廣。它與最短增廣路算法不同之處是:最短增廣路每個階段執行完一次BFS增廣後,要重新啟動BFS從源點Vs開始尋找另一條增廣路;而在Dinic算法中,隻需一次DFS過程就可以實作多次增廣,這是Dinic算法的巧妙之處。Dinic算法具體步驟如下:

(1)初始化容量網絡和網絡流。

(2)構造殘留網絡和層次網絡,若彙點不再層次網絡中,則算法結束。

(3)在層次網絡中用一次DFS過程進行增廣,DFS執行完畢,該階段的增廣也執行完畢。

(4)轉步驟(2)。

在Dinic的算法步驟中,隻有第(3)步與最短增廣路相同。在下面執行個體中,将會發現DFS過程将會使算法的效率有非常大的提高。

Dinic算法複雜度分析

與最短增廣路算法一樣,Dinic算法最多被分為n個階段,每個階段包括建層次網絡和尋找增廣路兩部分,其中建立層次網絡的複雜度仍是O(n*m)。

現在來分析DFS過程的總複雜度。在每一階段,将DFS分成兩部分分析。

(1)修改增廣路的流量并後退的花費。在每一階段,最多增廣m次,每次修改流量的費用為O(n)。而一次增廣後在增廣路中後退的費用也為O(n)。是以在每一階段中,修改增廣路以及後退的複雜度為O(m*n)。

(2)DFS周遊時的前進與後退。在DFS周遊時,如果目前路徑的最後一個頂點能夠繼續擴充,則一定是沿着第i層的頂點指向第i+1層頂點的邊向彙點前進了一步。因為增廣路經長度最長為n,是以最壞的情況下前進n步就會遇到彙點。在前進的過程中,可能會遇到沒有邊能夠沿着繼續前進的情況,這時将路徑中的最後一個點在層次圖中删除。

注意到每後退一次必定會删除一個頂點,是以後退的次數最多為n次。在每一階段中,後退的複雜度為O(n)。

假設在最壞的情況下,所有的點最後均被退了回來,一共共後退了n次,這也就意味着,有n次的前進被“無情”地退了回來,這n次前進操作都沒有起到“尋找增廣路”的作用。除去這n次前進和n次後退,其餘的前進都對最後找到增廣路做了貢獻。增廣路最多找到m次。每次最多前進n個點。是以所有前進操作最多為n+m*n次,複雜度為O(n*m)。

于是得到,在每一階段中,DFS周遊時前進與後退的花費為O(m*n)。

綜合以上兩點,一次DFS的複雜度為O(n*m)。是以,Dinic算法的總複雜度即O(m*n*n)。

下面給出 Dinic具體實作的代碼。

預處理:

[cpp]  view plain  copy

  1. #include <iostream>  
  2. #include <cstdlib>  
  3. #include <cstdio>  
  4. #include <cstring>  
  5. #include <string>  
  6. #include <algorithm>  
  7. #include <queue>  
  8. using namespace std;  
  9. const int MAXN = 210;  
  10. const int MAXM = 210*210;  
  11. const int INF = 0x3f3f3f3f;  
  12. struct Edge  
  13. {  
  14.     int v, f;  
  15.     int next;  
  16. }edge[MAXM];  
  17. int n, m;  
  18. int cnt;  
  19. int first[MAXN], level[MAXN];  
  20. int q[MAXN];  
  21. void init()  
  22. {  
  23.     cnt = 0;  
  24.     memset(first, -1, sizeof(first));  
  25. }  
  26. void read_graph(int u, int v, int f)  
  27. {  
  28.     edge[cnt].v = v, edge[cnt].f = f;  
  29.     edge[cnt].next = first[u], first[u] = cnt++;  
  30.     edge[cnt].v = u, edge[cnt].f = 0;  //增加一條反向弧,容量為0   
  31.     edge[cnt].next = first[v], first[v] = cnt++;  
  32. }  

BFS建構層次網絡:

[cpp]  view plain  copy

  1. int bfs(int s, int t) //建構層次網絡   
  2. {  
  3.     memset(level, 0, sizeof(level));  
  4.     level[s] = 1;  
  5.     int front = 0, rear = 1;  
  6.     q[front] = s;  
  7.     while(front < rear)  
  8.     {  
  9.         int x = q[front++];  
  10.         if(x == t) return 1;  
  11.         for(int e = first[x]; e != -1; e = edge[e].next)  
  12.         {  
  13.             int v = edge[e].v, f = edge[e].f;  
  14.             if(!level[v] && f)  
  15.             {  
  16.                 level[v] = level[x] + 1;  
  17.                 q[rear++] = v;  
  18.             }  
  19.         }  
  20.     }  
  21.     return 0;  
  22. }  

DFS找尋增廣路:

[cpp]  view plain  copy

  1. int dfs(int u, int maxf, int t)  
  2. {  
  3.     if(u == t) return maxf;  
  4.     int ret = 0;  
  5.     for(int e = first[u]; e != -1; e = edge[e].next)  
  6.     {  
  7.         int v = edge[e].v, f = edge[e].f;  
  8.         if(level[u] + 1 == level[v] && f)  
  9.         {  
  10.             int Min = min(maxf-ret, f);  
  11.             f = dfs(v, Min, t);  
  12.             edge[e].f -= f;  
  13.             edge[e^1].f += f;  
  14.             ret += f;  
  15.             if(ret == maxf) return ret;  
  16.         }  
  17.     }  
  18.     return ret;  
  19. }  

Dinic:

[cpp]  view plain  copy

  1. int Dinic(int s, int t) //Dinic  
  2. {  
  3.     int ans = 0;  
  4.     while(bfs(s, t)) ans += dfs(s, INF, t);  
  5.     return ans;  
  6. }  

層次圖

最大流Dinic算法講解 &amp;&amp; ISAP 算法解釋短路增值算法(MPLA)Dinic算法 網絡流-最大流問題 ISAP 算法解釋
最大流Dinic算法講解 &amp;&amp; ISAP 算法解釋短路增值算法(MPLA)Dinic算法 網絡流-最大流問題 ISAP 算法解釋

短路增值算法(MPLA)

算法流程

最大流Dinic算法講解 &amp;&amp; ISAP 算法解釋短路增值算法(MPLA)Dinic算法 網絡流-最大流問題 ISAP 算法解釋

彙點不在層次圖内意味着在剩餘圖中不存在一條從源點到彙點的路徑,即沒有增廣路。

在程式實作的時候,層次圖并不用被“建”出來,我們隻需對每個頂點标記層次,增廣的時候,判斷邊是否滿足level(u) +1= level(v)這一限制即可。

Dinic算法

Dinic算法的思想也是分階段地在層次圖中增廣。

它與最短路徑增值算法不同之處是:在Dinic算法中,我們用一個dfs過程代替多次bfs來尋找阻塞流。下面給出其算法步驟:

算法流程

最大流Dinic算法講解 &amp;&amp; ISAP 算法解釋短路增值算法(MPLA)Dinic算法 網絡流-最大流問題 ISAP 算法解釋

增廣過程圖解

最大流Dinic算法講解 &amp;&amp; ISAP 算法解釋短路增值算法(MPLA)Dinic算法 網絡流-最大流問題 ISAP 算法解釋
最大流Dinic算法講解 &amp;&amp; ISAP 算法解釋短路增值算法(MPLA)Dinic算法 網絡流-最大流問題 ISAP 算法解釋

僞代碼描述

最大流Dinic算法講解 &amp;&amp; ISAP 算法解釋短路增值算法(MPLA)Dinic算法 網絡流-最大流問題 ISAP 算法解釋

在程式裡,p表示找到的增廣路徑,p.top為路徑中的最後一個頂點。一開始,p中隻有源點。

整個While循環分為2個操作。如果p的最後一個頂點為彙點,也就是說找到了增廣路,那麼對p增廣,注意到增廣後一定有一條或多條p中的邊被删除了。這時,我們使增廣路徑後退至p中從源點可到達的最後一個頂點。

如果p的最後一個頂點不為彙點,那麼觀察最後那個的頂點u 。若在層次圖中存在從u連出的一條邊,比如(u,v),我們就将頂點v放入路徑p中,繼續dfs周遊;否則,點u對之後的dfs周遊就沒有用了,我們将點u以及層次圖中連到u的所有邊删除,并且在p中後退一個點。

Dfs過程将會不斷重複這2個操作,直到從源點連出的邊全部被删除為止。

轉自:http://www.renfei.org/blog/isap.html

網絡流-最大流問題 ISAP 算法解釋

August 7, 2013 /  程式設計指南

ISAP 是圖論求最大流的算法之一,它很好的平衡了運作時間和程式複雜度之間的關系,是以非常常用。

約定

我們使用鄰接表來表示圖,表示方法可以見文章帶權最短路 Dijkstra, SPFA, Bellman-Ford, ASP, Floyd-Warshall 算法分析或二分圖的最大比對、完美比對和匈牙利算法的開頭(就不重複貼代碼了)。在下文中,圖的源點(source)表示為  s s ,彙點(sink)表示為  t t ,目前節點為  u u 。建圖時,需要建立雙向邊(設反向的邊容量為  0 0)才能保證算法正确。

引入

求解最大流問題的一個比較容易想到的方法就是,每次在殘量網絡(residual network)中任意尋找一條從  s s 到  t t 的路徑,然後增廣,直到不存在這樣的路徑為止。這就是一般增廣路算法(labeling algorithm)。可以證明這種不加改進的貪婪算法是正确的。假設最大流是  f f,那麼它的運作時間為 O( f⋅∣E∣) O( f⋅∣E∣) 。但是,這個運作時間并不好,因為它和最大流  f f 有關。

人們發現,如果每次都沿着殘量網絡中的最短增廣路增廣,則運作時間可以減為  O(∣E∣2⋅∣V∣)  O(∣E∣2⋅∣V∣)  。這就是最短增廣路算法。而 ISAP 算法則是最短增廣路算法的一個改進。其實,ISAP 的意思正是「改進的最短增廣路」 (Improved Shortest Augmenting Path)。

順便說一句,上面讨論的所有算法根本上都屬于增廣路方法(Ford-Fulkerson method)。和它對應的就是大名鼎鼎的預流推進方法(Preflow-push method)。其中最高标号預流推進算法(Highest-label preflow-push algorithm)的複雜度可以達到  O(∣V∣2∣E∣−−−√) O(∣V∣2∣E∣)。雖然在複雜度上比增廣路方法進步很多,但是預流推進算法複雜度的上界是比較緊的,是以有時差距并不會很大。

算法解釋

概括地說,ISAP 算法就是不停地找最短增廣路,找到之後增廣;如果遇到死路就 retreat,直到發現 s s,  t t不連通,算法結束。找最短路本質上就是無權最短路徑問題,是以采用 BFS 的思想。具體來說,使用一個數組 d d,記錄每個節點到彙點 t t的最短距離。搜尋的時候,隻沿着滿足 d[u]=d[v]+1 d[u]=d[v]+1的邊 u→v u→v(這樣的邊稱為允許弧)走。顯然,這樣走出來的一定是最短路。

原圖存在兩種子圖,一個是殘量網絡,一個是允許弧組成的圖。殘量網絡保證可增廣,允許弧保證最短路(時間界較優)。是以,在尋找增廣路的過程中,一直是在殘量網絡中沿着允許弧尋找。是以,允許弧應該是屬于殘量網絡的,而非原圖的。換句話說,我們沿着允許弧,走的是殘量網絡(而非原圖)中的最短路徑。當我們找到沿着殘量網絡找到一條增廣路,增廣後,殘量網絡肯定會變化(至少少了一條邊),是以決定允許弧的 d d數組要進行相應的更新(順便提一句,Dinic 的做法就是每次增廣都重新計算 d d數組)。然而,ISAP 「改進」的地方之一就是,其實沒有必要馬上更新 d d數組。這是因為,去掉一條邊隻可能令路徑變得更長,而如果增廣之前的殘量網絡存在另一條最短路,并且在增廣後的殘量網絡中仍存在,那麼這條路徑毫無疑問是最短的。是以,ISAP 的做法是繼續增廣,直到遇到死路,才執行 retreat 操作。

說到這裡,大家應該都猜到了,retreat 操作的主要任務就是更新 d d數組。那麼怎麼更新呢?非常簡單:假設是從節點 u u找遍了鄰接邊也沒找到允許弧的;再設一變量 m m,令 m m等于殘量網絡中 u u的所有鄰接點的 d d數組的最小值,然後令 d[u] d[u]等于 m+1 m+1即可。這是因為,進入 retreat 環節說明殘量網絡中 u u和  t t已經不能通過(已過時)的允許弧相連,那麼 u u和 t t實際上在殘量網絡中的最短路的長是多少呢?(這正是 d d的定義!)顯然是殘量網絡中 u u的所有鄰接點和 t t的距離加 1 1的最小情況。特殊情況是,殘量網絡中 u u根本沒有鄰接點。如果是這樣,隻需要把 d[u] d[u]設為一個比較大的數即可,這會導緻任何點到 u u的邊被排除到殘量網絡以外。(嚴格來說隻要大于等于 ∣V∣ ∣V∣即可。由于最短路一定是無環的,是以任意路徑長最大是 ∣V∣−1 ∣V∣−1)。修改之後,隻需要把正在研究的節點 u u沿着剛才走的路退一步,然後繼續搜尋即可。

講到這裡,ISAP 算法的架構内容就講完了。對于代碼本身,還有幾個優化和實作的技巧需要說明。

  1. 算法執行之前需要用 BFS 初始化 d d數組,方法是從 t t到 s s逆向進行。
  2. 算法主體需要維護一個「目前節點」 u u,執行這個節點的前進、retreat 等操作。
  3. 記錄路徑的方法非常簡單,聲明一個數組 p p,令 p[i] p[i]等于增廣路上到達節點 i i的邊的序号(這樣就可以找到從哪個頂點到的頂點 i i)。需要路徑的時候反向追蹤一下就可以了。
  4. 判斷殘量網絡中 s,t s,t不連通的條件,就是 d[s]≥∣V∣ d[s]≥∣V∣ 。這是因為當 s,t s,t不連通時,最終殘量網絡中 s s将沒有任何鄰接點,對 s s的 retreat 将導緻上面條件的成立。
  5. GAP 優化。GAP 優化可以提前結束程式,很多時候提速非常明顯(高達 100 倍以上)。GAP 優化是說,進入 retreat 環節後, u,t u,t之間的連通性消失,但如果 u u是最後一個和 t t距離 d[u] d[u](更新前)的點,說明此時 s,t s,t也不連通了。這是因為,雖然 u,t u,t已經不連通,但畢竟我們走的是最短路,其他點此時到 t t的距離一定大于 d[u] d[u](更新前),是以其他點要到 t t,必然要經過一個和 t t距離為 d[u] d[u](更新前)的點。GAP 優化的實作非常簡單,用一個數組記錄并在适當的時候判斷、跳出循環就可以了。
  6. 另一個優化,就是用一個數組儲存一個點已經嘗試過了哪個鄰接邊。尋找增廣的過程實際上類似于一個 BFS 過程,是以之前處理過的鄰接邊是不需要重新處理的(殘量網絡中的邊隻會越來越少)。具體實作方法直接看代碼就可以,非常容易了解。需要注意的一點是,下次應該從上次處理到的鄰接邊繼續處理,而非從上次處理到的鄰接邊的下一條開始。

最後說一下增廣過程。增廣過程非常簡單,尋找增廣路成功(目前節點處理到 t t)後,沿着你記錄的路徑走一遍,記錄一路上的最小殘量,然後從 s s到 t t更新流量即可。

實作

int source;         // 源點
int sink;           // 彙點
int p[max_nodes];   // 可增廣路上的上一條弧的編号
int num[max_nodes]; // 和 t 的最短距離等于 i 的節點數量
int cur[max_nodes]; // 目前弧下标
int d[max_nodes];   // 殘量網絡中節點 i 到彙點 t 的最短距離
bool visited[max_nodes];

// 預處理, 反向 BFS 構造 d 數組
bool bfs()
{
    memset(visited, 0, sizeof(visited));
    queue<int> Q;
    Q.push(sink);
    visited[sink] = 1;
    d[sink] = 0;
    while (!Q.empty()) {
        int u = Q.front();
        Q.pop();
        for (iterator_t ix = G[u].begin(); ix != G[u].end(); ++ix) {
            Edge &e = edges[(*ix)^1];
            if (!visited[e.from] && e.capacity > e.flow) {
                visited[e.from] = true;
                d[e.from] = d[u] + 1;
                Q.push(e.from);
            }
        }
    }
    return visited[source];
}

// 增廣
int augment()
{
    int u = sink, df = __inf;
    // 從彙點到源點通過 p 追蹤增廣路徑, df 為一路上最小的殘量
    while (u != source) {
        Edge &e = edges[p[u]];
        df = min(df, e.capacity - e.flow);
        u = edges[p[u]].from;
    }
    u = sink;
    // 從彙點到源點更新流量
    while (u != source) {
        edges[p[u]].flow += df;
        edges[p[u]^1].flow -= df;
        u = edges[p[u]].from;
    }
    return df;
}

int max_flow()
{
    int flow = 0;
    bfs();
    memset(num, 0, sizeof(num));
    for (int i = 0; i < num_nodes; i++) num[d[i]]++;
    int u = source;
    memset(cur, 0, sizeof(cur));
    while (d[source] < num_nodes) {
        if (u == sink) {
            flow += augment();
            u = source;
        }
        bool advanced = false;
        for (int i = cur[u]; i < G[u].size(); i++) { 
            Edge& e = edges[G[u][i]];
            if (e.capacity > e.flow && d[u] == d[e.to] + 1) {
                advanced = true;
                p[e.to] = G[u][i];
                cur[u] = i;
                u = e.to;
                break;
            }
        }
        if (!advanced) { // retreat
            int m = num_nodes - 1;
            for (iterator_t ix = G[u].begin(); ix != G[u].end(); ++ix)
                if (edges[*ix].capacity > edges[*ix].flow)
                    m = min(m, d[edges[*ix].to]);
            if (--num[d[u]] == 0) break; // gap 優化
            num[d[u] = m+1]++;
            cur[u] = 0;
            if (u != source)
                u = edges[p[u]].from;
        }
    }
    return flow;
}