圖論演算法小結:網路流及其演算法
個人說明:最近學到了圖論演算法,但網路流這部分頗難理解,於是在網上找到了一片比較好的講解部落格。轉載之~
將每條有向邊想象成傳輸物質的管道。每個管道都有一個固定的容量,可以看作是物質能流經該管道的最大速度(譬如可以想象為水流和河槽)。頂點是管道間的連線點,除了源點(S,Source)和匯點(T,Target)以外,物質只流經這些頂點。而不聚集在頂點中。
注:下文提到的數字,基本都可加單位 “單位流量”來理解。
一、網路流基礎
網路流中的最大流研究的問題是:在不違背容量限制的條件下,把物質從源點傳輸到匯點的最大速率是多少?上邊的圖示中,還有些需要說明,設G = (V,E)是一個流網路,其容量函式為c(c(u, v),Capicity)。s 為源點,t 為匯點。G 的流是一個函式 f (Flow)。上圖中每條邊標記的是流網路的容量 c(u, v),右圖中,G中的流 |f| = 19。圖中只顯示正網路流。如果 f(u,v) > 0,則標記(u,v)為 f(u,v)/ c(u,v)(斜槓僅僅用來區分流和容量,不表示相除)。如果 f(u,v)<= 0,邊(u,v)只標記它的容量。譬如(s,v1
網路流演算法要基於三種思想:殘留網路(Residual Network),增廣路徑(Augmenting Path)和割(Cut)。
和開頭同一個例子,上圖,圖(b)就是殘留網路,仔細觀察一下,應該就可以明白,譬如邊(s,v1),從 s 到 v1 流過11,還剩5,所以殘餘容量是5,當然,流量是守恆的,所以反向從 v1 到 s 流過11。剩餘的容量 + 反向平衡的流量共同構成了殘留網路。對於名詞“殘留容量(Residual Capacity)”的定義:在不超過容量c(u,v)的條件下,從 u 到 v 之間可以壓入的額外網路流量,就是(u,v)的殘留容量(Residual Capacity),公式定義是:cf
圖(b),陰影覆蓋的邊為增廣路徑 p ,其殘留容量為 cf(p) = cf(v2,v3) = 4。看圖應該就可以大概理解什麼是增廣路了。看圖可以發現 s 到 v2還可以通過 5,v2到v3還可以通過 4,v3到 t 還可以流 5,按照常識,也可得出,這條路徑還可以再流過 4 的結論。而從 s 輸送 4 到 t 的這條路就是增廣路。"增廣路徑 p”的定義:p 為殘留網路Gf中從 s 到 t 的一條簡單路徑。在不違反邊的容量限制條件下,增廣路徑上的每條邊(u,v)可以容納從 u 到 v 的某額外正網路流。Starvae師兄在講解增廣路時,有一段通俗的描述:
假如有這麼一條路,這條路從源點開始一直一段一段的連到了匯點,並且,這條路上的每一段都滿足流量<容量,注意,是嚴格的<,而不是<=。那麼,我們一定能找到這條路上的每一段的(容量-流量)的值當中的最小值delta。我們把這條路上每一段的流量都加上這個delta,一定可以保證這個流依然是可行流。這樣我們就得到了一個更大的流,他的流量是之前的流量+delta,而這條路就叫做增廣路。
另外,這裡得出的,計算增廣路流量的公式非常重要:cf(p) = min { cf(u,v):(u,v)在 p 上},因為最大流演算法的核心基本也就是尋找增廣路了。
最大流最小割定理:一個流是最大流,當且僅當它的殘留網路不包含增廣路徑。
流網路的割(S,T)將V劃分為 S 和 T = V – S兩部分,使得 s ∈ S,t ∈ T。如果 f 是一個流,則穿過割 (S,T)的淨流被定義為 f (S,T)。割(S,T)的容量為 c(S,T)。一個網路的最小割就是網路中所有割中最小容量的割。
上圖中,流網路的割(S = {s,v1,v2},T = {v3,v4,t}),S中的頂點是黑色,T中的頂點是白色。穿越(S,T)的淨流量為:f(v1,v3) + f(v2,v3) + f(v2,v4) = 12 + (-4) + 11 = 19;容量為:c(v1,v3) + c(v2,v4) = 12 + 14 = 26。可以看出,割的淨流由雙向的網路流組成。而割的容量僅由 S 到 T 的邊計算而得。
還有一個很有重要的知識:反向邊。如圖(a),1是源點,4是匯點。很明顯如果第一次迭代找到的增廣路是1→2→4和1→3→4,則最大流是2,但是如果第一次迭代找到的增廣路是1→2→3→4,那麼流量只有1,如圖(b)是殘留網路,這時候,反向邊就起作用了,由於反向邊的原因,第二次迭代的時候,又找到一條增廣路1→2→3→4。這樣下來,總的流量還是2。還是Starvae師兄的解釋:
當我們第二次的增廣路走 3→2 這條反向邊的時候,就相當於把 2→3 這條正向邊已經是用了的流量給“退”了回去,不走 2→3 這條路,而改走從 2 點出發的其他的路也就是 2→4。(有人問如果這裡沒有 2→4 怎麼辦,這時假如沒有 2→4 這條路的話,最終這條增廣路也不會存在,因為他根本不能走到匯點)同時本來在 3→4 上的流量由 1→3→4 這條路來“接管”。而最終 2→3 這條路正向流量1,反向流量1,等於沒有流量。反向邊的作用就是給程式一個可以後悔的機會。—— 一語中的啊。(這句話我加的^^)
二、Ford-Fulkerson演算法
Ford-Fulkerson演算法在實際中並不常用,但是它提供了一種思想:先找到一條從源點到匯點的增廣路徑,這條路徑的容量是其中容量最小的邊的容量。然後,通過不斷找增廣路,一步步擴大流量,當找不到增廣路時,就得到最大流了(最大流最小割定理)。可以看看Ford-Fulkerson演算法的虛擬碼,
FORD-FULKERSON(G, s, t)
1 for each edge (u, v) ∈ E[G]
2 do f[u, v] ← 0
3 f[v, u] ← 0
4 while there exists a path p from s to t in the residual network Gf
5 do cf(p) ← min {cf(u, v) : (u, v) is in p}
6 for each edge (u, v) in p
7 do f[u, v] ← f[u, v] + cf(p)
8 f[v, u] ← -f[u, v]
從虛擬碼可以看出,Ford-Fulkerson的核心過程就是:第4~8行的 while 迴圈反覆找出 Gf 中的增廣路徑 p,並把沿 p 的流 f 加上其殘留容量cf(p)。當不再有增廣路徑時,流 f 就是一個最大流。完整圖示:
(f)最後在 while 迴圈測試的殘留網路,它沒有增光路徑,所以(e)中顯示的流就是最大流。
三、壓入重標記push_relabel演算法O(VE)
Push-relabel用到一個很有趣的概念一Preflow(前置流)他允許流進的量比流出的量還要多,有水來就先流進來,流不出去再說。
Push-relabel演算法的流程如下:
1 )我們先假定一個高度函式 h ( u ) ,他代表u點的高度,只有
h ( u )比較高的點才能夠將水流到h ( u )比較低的點;
2 ) 在程式一開始的時候,讓source node的高度是n ( node數) , 其它點的高度都是 0 ,這樣source node才有足夠的高度可以流往其它地方 ;
3 ) 然後, 讓source node往其它所有跟他直接相鄰的node ,都流水管寬度的水量 ( 流過去之後當然要計算剩下的網路情形,即計算 residual edges(殘留網路)
4 )對所有active node( 目前有水的 node )做 relabel(重標記)的動作: 在當某個node明明有水,但是他所連出去的所有物件的 h ( u ) 都比他還高, 則讓他的h ( u ) 增加為至少有一條水管可以流出去的量,也就是讓這個有水的 active node的高度變成比他連往的”高度最小的 n o d e ”+l , ( 流過去之後還是要計算剩下的網路情形
5 )對所有可以做push(壓入)動作的node做push的動作。 所謂的Push動作是指 : 當某個node有水,並且他有可以流出去的邊, 且他剛好比可以流出去的那個點高度高一點點 ( 高度恰好比他高 1 ) ,那就把某個node 的水流過去,要流多少呢?以下兩者取 min。 流出去的水管的量( 也就是說,這個active node的水量很多, 足夠把這條水管塞滿( 飽和) ,這個時候就叫做saturating Push)(飽和壓入)某個 n o d e 現在的水量( 這個 n o d e的水量不足以把流出去的這個水管填滿,稱作non saturating Push(不飽和壓入)
6 )重複Relabel和Push的工作,一直到沒有active node為止,此時從source node所流出的總流量( P r e f l o w) ,就是這個圖的最大流量。
Push-relabel algorithm 提供了最大流另一方向的思考,且就效率而言,Push -Relabel的複雜度為 o (vE )
#define MIN INT_MIN
#define MAX INT_MAX
#define N 110
int min(int a,int b){return a>b?b:a;}
int c[N][N];//殘留容量
int ef[N];//頂點餘流
int h[N];//頂點高度
int n;
int push_relabel(int s,int t){
int i,j;
int ans = 0;
memset(h,0,sizeof(h));
h[s] = t+1;//源點初始高度
memset(ef,0,sizeof(ef));
ef[s] = MAX;//源點初始餘流
queue<int> qq;
qq.push(s);
while(!qq.empty()){
int u = qq.front();
qq.pop();
for(i=0;i<=t;i++){
int p;
int v = i;
if(c[u][v]<ef[u])p = c[u][v];
else p = ef[u];
if(p>0 && (u==s || h[u] == h[v] +1)){
c[u][v] -= p;
c[v][u] += p;
if(v==t)ans+=p;//如果到達了匯點,就將流值加入到最大流中
ef[u] -= p;
ef[v] += p;
if(v!=s && v!=t)qq.push(v);//只有既不是源點也不是匯點才進隊
}
}
//如果不是源點且仍有餘流,則重標記高度再進隊。
//這裡只是簡單的將高度增加了一個單位,也可以像上面所說的一樣賦值為最低的相鄰頂點的高度高一個單位
if(u!= s && u!=t && ef[u]>0) {
h[u]++;
qq.push(u);
}
}
return ans;
}
int main(){
int np,nc,m;
while(scanf("%d%d%d%d",&n,&np,&nc,&m) != -1){
int s = n,t = n+1;
int i,j;
memset(c,0,sizeof(c));
char ss[30];
for(i=0;i<m;i++){
int u,v,w;
scanf("%s",ss);
sscanf(ss,"(%d,%d)%d",&u,&v,&w);
c[u][v] += w;
}
for(i=0;i<np;i++){
int u,w;
scanf("%s",ss);
sscanf(ss,"(%d)%d",&u,&w);
c[s][u] += w;
}
for(i=0;i<nc;i++){
int v,w;
scanf("%s",ss);
sscanf(ss,"(%d)%d",&v,&w);
c[v][t] += w;
}
printf("%d\n",push_relabel(s,t));
}
return 0;
}
四、Edmonds-Karp(EK)演算法 O(V*E*E)
EK演算法基於Ford-Fulkerson演算法,唯一的區別是將第 4 行用BFS(廣度優先搜尋)來實現對增廣路徑 p 的計算。EK演算法虛擬碼基本和上邊的Ford-Fulkerson演算法一樣。類似用DFS實現的還有Dinic演算法。它們都屬於SAP(Shortest Augmenting Path)演算法,從英文即可看出,它們每次都在尋找最短增廣路。對於EK演算法,每次用一遍 BFS 尋找從源點 s 到終點 t 的最短路作為增廣路徑,然後增廣流量 f 並修改殘量網路,直到不存在新的增廣路徑。E-K 演算法的時間複雜度為 O(VE^2),適用於稀疏邊,由於 BFS 要搜尋全部小於最短距離的分支路徑之後才能找到終點,因此頻繁的 BFS 效率是比較低的。實踐中此演算法使用的機會較少。
這裡以POJ 1273為例,這裡可以作為EK模板,也是我的第一道網路流,mark一下
#define MIN INT_MIN
#define MAX INT_MAX
#define N 204
int c[N][N];//邊容量
int f[N][N];//邊實際流量
int pre[N];//記錄增廣路徑
int res[N];//殘餘網路
queue<int> qq;
void init(){
while(!qq.empty())qq.pop();
memset(c,0,sizeof(c));
memset(f,0,sizeof(f));
}
int EK(int s,int t){
int i,j;
int ans=0;
while(1){
memset(res,0,sizeof(res));
res[s] = MAX;//源點的殘留網路要置為無限大!否則下面找增廣路出錯
pre[s] = -1;
qq.push(s);
//bfs找增廣路徑
while(!qq.empty()){
int x = qq.front();
qq.pop();
for(i=1;i<=t;i++){
if(!res[i] && f[x][i] < c[x][i]){
qq.push(i);
pre[i] = x;
res[i] = min(c[x][i] - f[x][i], res[x]);//這裡類似dp,如果有增廣路,那麼res[t]就是增廣路的最小權
}
}
}
if(res[t]==0)break;//找不到增廣路就退出
int k = t;
while(pre[k]!=-1){
f[pre[k]][k] += res[t];//正向邊加上新的流量
f[k][pre[k]] -= res[t];//反向邊要減去新的流量,反向邊的作用是給程式一個後悔的機會
k = pre[k];
}
ans += res[t];
}
return ans;
}
int main(){
int n,m;
while(scanf("%d%d",&n,&m) != -1){
int i,j;
init();
while(n--){
int a,b,v;
scanf("%d%d%d",&a,&b,&v);
c[a][b]+=v;
}
printf("%d\n",EK(1,m));
}
return 0;
}
四、Improved SAP(ISAP)演算法
ISAP字面意思是改良的最短增廣路演算法。關於ISAP,一位叫 DD_engi 的神牛講非常清楚,引用一下:
SAP演算法(by dd_engi):求最大流有一種經典的演算法,就是每次找增廣路時用BFS找,保證找到的增廣路是弧數最少的,也就是所謂的 Edmonds-Karp 演算法。可以證明的是在使用最短路增廣時增廣過程不超過 V * E次,每次 BFS 的時間都是O(E),所以 Edmonds-Karp 的時間複雜度就是O(V * E^2)。
如果能讓每次尋找增廣路時的時間複雜度降下來,那麼就能提高演算法效率了,使用距離標號的最短增廣路演算法就是這樣的。所謂距離標號,就是某個點到匯點的最少的弧的數量(另外一種距離標號是從源點到該點的最少的弧的數量,本質上沒什麼區別)。設點 i 的標號為D[i],那麼如果將滿足D[i] = D[j] + 1的弧(i,j))叫做允許弧,且增廣時只走允許弧,那麼就可以達到“怎麼走都是最短路”的效果。每個點的初始標號可以在一開始用一次從匯點沿所有反向邊的BFS求出,實踐中可以初始設全部點的距離標號為0,問題就是如何在增廣過程中維護這個距離標號。
維護距離標號的方法是這樣的:當找增廣路過程中發現某點出發沒有允許弧時,將這個點的距離標號設為由它出發的所有弧的終點的距離標號的最小值加一。這種維護距離標號的方法的正確性我就不證了。由於距離標號的存在,由於“怎麼走都是最短路”,所以就可以採用DFS找增廣路,用一個棧儲存當前路徑的弧即可。當某個點的距離標號被改變時,棧中指向它的那條弧肯定已經不是允許弧了,所以就讓它出棧,並繼續用棧頂的弧的端點增廣。為了使每次找增廣路的時間變成均攤O(V),還有一個重要的優化是對於每個點儲存“當前弧”:初始時當前弧是鄰接表的第一條弧;在鄰接表中查詢時從當前弧開始查詢,找到了一條允許弧,就把這條弧設為當前弧;改變距離標號時,把當前弧重新設為鄰接表的第一條弧,還有一種在常數上有所優化的寫法是改變距離標號時把當前弧設為那條提供了最小標號的弧。當前弧的寫法之所以正確就在於任何時候我們都能保證在鄰接表中當前弧的前面肯定不存在允許弧。
還有一個常數優化是在每次找到路徑並增廣完畢之後不要將路徑中所有的頂點退棧,而是隻將瓶頸邊以及之後的邊退棧,這是借鑑了Dinic演算法的思想。注意任何時候待增廣的“當前點”都應該是棧頂的點的終點。這的確只是一個常數優化,由於當前邊結構的存在,我們肯定可以在O(n)的時間內復原路徑中瓶頸邊之前的所有邊。
優化:
1.鄰接表優化:
如果頂點多的話,往往N^2存不下,這時候就要存邊:
存每條邊的出發點,終止點和價值,然後排序一下,再記錄每個出發點的位置。以後要呼叫從出發點出發的邊時候,只需要從記錄的位置開始找即可(其實可以用連結串列)。優點是時間加快空間節省,缺點是程式設計複雜度將變大,所以在題目允許的情況下,建議使用鄰接矩陣。
2.GAP優化:
如果一次重標號時,出現距離斷層,則可以證明ST無可行流,此時則可以直接退出演算法。
3.當前弧優化:
為了使每次找增廣路的時間變成均攤O(V),還有一個重要的優化是對於每個點儲存“當前弧”:初始時當前弧是鄰接表的第一條弧;在鄰接表中查詢時從當前弧開始查詢,找到了一條允許弧,就把這條弧設為當前弧;改變距離標號時,把當前弧重新設為鄰接表的第一條弧。
另外,ISAP簡化的描述是:程式開始時用一個反向 BFS 初始化所有頂點的距離標號,之後從源點開始,進行如下三種操作:(1)當前頂點 i 為終點時增廣 (2) 當前頂點有滿足 dist[i] = dist[j] + 1 的出弧時前進 (3) 當前頂點無滿足條件的出弧時重標號並回退一步。整個迴圈當源點 s 的距離標號 dist[s] >= n 時結束。對 i 點的重標號操作可概括為 dist[i] = 1 + min{dist[j] : (i,j)屬於殘量網路Gf}。
const int MAXN=20010;
const int MAXM=500010;
int n,m;//n為點數 m為邊數
int h[MAXN];
int gap[MAXN];
int p[MAXN],ecnt;
int source,sink;
struct edge{
int v;//邊的下一點
int next;//下一條邊的編號
int val;//邊權值
}e[MAXM];
inline void init(){memset(p,-1,sizeof(p));eid=0;}
//有向
inline void insert1(int from,int to,int val){
e[ecnt].v=to;
e[ecnt].val=val;
e[ecnt].next=p[from];
p[from]=eid++;
swap(from,to);
e[ecnt].v=to;
e[ecnt].val=0;
e[ecnt].next=p[from];
p[from]=eid++;
}
//無向
inline void insert2(int from,int to,int val){
e[ecnt].v=to;
e[ecnt].val=val;
e[ecnt].next=p[from];
p[from]=eid++;
swap(from,to);
e[ecnt].v=to;
e[ecnt].val=val;
e[ecnt].next=p[from];
p[from]=eid++;
}
inline int dfs(int pos,int cost){
if (pos==sink){
return cost;
}
int j,minh=n-1,lv=cost,d;
for (j=p[pos];j!=-1;j=e[j].next){
int v=e[j].v,val=e[j].val;
if(val>0){
if (h[v]+1==h[pos]){
if (lv<e[j].val) d=lv;
else d=e[j].val;
d=dfs(v,d);
e[j].val-=d;
e[j^1].val+=d;
lv-=d;
if (h[source]>=n) return cost-lv;
if (lv==0) break;
}
if (h[v]<minh) minh=h[v];
}
}
if (lv==cost){
--gap[h[pos]];
if (gap[h[pos]]==0) h[source]=n;
h[pos]=minh+1;
++gap[h[pos]];
}
return cost-lv;
}
int sap(int st,int ed){
source=st;
sink=ed;
int ans=0;
memset(gap,0,sizeof(gap));
memset(h,0,sizeof(h));
gap[st]=n;
while (h[st]<n){
ans+=dfs(st,INT_MAX);
}
return ans;
}
對於EK演算法與ISAP演算法的區別:
EK演算法每次都要重新尋找增廣路,尋找過程只受殘餘網路的影響,如果改變殘餘網路,則增廣路的尋找也會隨之改變;SAP演算法預處理出了增廣路的尋找大致路徑,若中途改變殘餘網路,則此演算法將重新進行。EK處理在運算過程中需要不斷加邊的最大流比SAP更有優勢。
*****************************************************************************************************************************************************
本文只介紹基礎的網路流知識,網路流非常強大,許多問題都可以轉化為網路流模型,進一步的,可以去看看Starfall大神的這篇《【網路流】總結》,這篇只是一個目錄性質,裡邊很多超連結,後面資源更豐富。其中提到了6篇國家集訓隊論文,正好電腦裡都有,就打包傳到CSDN了,點選跳轉到下載頁面。我只看過其中的兩三篇,這幫高中生寫的實在牛叉。
記得當年看時,覺得這篇《最大流在資訊學競賽中應用的一個模型》當做入門非常好,看完後會發現,原來組合數學都可以用最大流來解,還有什麼不可以的。另外,最經典,最全面的要數胡波濤這篇的《最小割模型在資訊學競賽中的應用》,如果想深入學習網路流,這幫傢伙的論文絕對不能錯過。