スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Ford-Fulkerson のアルゴリズム

最大流問題を解くための Ford-Fulkerson のアルゴリズムを復習したので、簡単にまとめてみる。

Tags : グラフ理論 アルゴリズム プログラミング

※ 間違いがあるかもしれないので、鵜呑みにしないほうがいいです。

Ford-Fulkerson は、フローネットワークにおける最大フローを求めるためのアルゴリズム

以下、この記事を通して次の記号を使う。
・ G : フローネットワーク (ソース、シンクをもつ有向グラフで各辺に容量が与えられたもの)
・ n : G の頂点数
・ m : G の辺数
・ src : G のソース
・ snk : G のシンク



擬似コード (簡略版)

totalflow = 0
while src から snk への augment path がある
  augment path に沿って流せるだけフローを流す
  totalflow += 流したフロー
end while
return totalflow


一言でいえば、流せなくなるまでフローを流す Greedy アルゴリズムである。



いくつかの用語と重要そうな事項

G 上のフロー (flow) とは、G の辺から実数への関数である。
特に、値域は非負整数であることが多い。
正確にはフロー関数というべき。
厳密な定義は省略。たとえば Wikipedia - フローネットワーク にある。
特定の辺 e に割り当てられたフローの値についていいたいときは、e 上のフローなどという。
G の辺をパイプ、頂点をジャンクションとみたときに、フローはパイプ中を流れる水に対応する。
なので、グラフ中の有向道 R に沿ってフローの値を増やすときは、R に沿ってフローを流すなどという。

最大流問題とは、G 上に最大でどれだけのフローが流せるかを求める問題。

G の残余グラフ (residual graph) とは、フローを流す余地のある辺からなるグラフのこと。
「フローを流す余地のある辺」は次の 2 通りが考えられる。
・ G の辺で飽和していないもの (辺の容量 - 辺のフロー > 0 となっているもののこと; 空き容量があるイメージ)
・ G の辺 e と逆向きの辺で、e にフローが流れているもの (流れているフローを押し返すイメージ)
残余グラフを用いるとアルゴリズムの見通しがよくなる。

増加道 (augment path) とは、残余グラフの辺からなる、src から snk への有向道のこと。
augment path に沿ってフローを流すと、G 上のフローが増える。
このような操作をしても、更新後のフローは、フローの定義をみたしていることに注意。
特に整数フロー (辺容量がすべて整数) のときは、while ループ 1 回ごとにフローが 1 以上増えるので、アルゴリズムは必ず終了する。
このことから、整数フローではアルゴリズムの時間計算量は O(mC) で抑えられる。
ここで、C は G の最大フロー。

実数フローでは、Ford-Fulkerson は終了しない場合がある。
これは、フローの増加量がだんだん小さくなっていくような場合である。
何かの本で具体的に構成した例を見たことがあるが、書名を忘れた。

Ford-Fulkerson は augment path の探し方によって複数の実装が考えられ、それぞれ特性がある。

DFS で探索した場合、1 回あたりの探索時間は短くなりがちだが、良くない augment path を選んでしまうことがあり、最悪で O(mC) 時間かかる。
ただし、二部マッチングに応用するときなど、フローの増加量が 1 であることが分かっている場合には BFS より良いかもしれない。

BFS で探索したものは、特に Edmonds-Karp のアルゴリズムと呼ばれる。
Edmonds-Karp の時間計算量は O(m2n) であることがわかっている。
なお、Edmonds-Karp は実数フローでも必ず終了する。

一般に、フローの増加量が大きい augment path を優先的に選ぶことができれば、アルゴリズムの反復回数は少なく済むだろう。
このアイデアに基づいた、Capacity Scaling などと呼ばれる探索方法もある。
これは、パラメータ ⊿ を用いて、フローの増加量が ⊿ 以上になる augment path を優先的に見つける方法である。
所望の augment path が見つからなくなったら、⊿ の値を下げて、探索を繰り返す。
この方法の詳細は、『アルゴリズムデザイン』 7.3 節が詳しい。



C++ によるコード例

注意事項
ここのコードをテンプレとして使っている。
・ T は辺容量とフローの型。
・ AdjMatrix<T> は vector< vector<T> > と同じ。
・ Capacity Scaling による実装は整数フローのみ対応。

単純な DFS / BFS による実装
// DFS バージョン
template
<class T>
T findAugpath(int u,vi &aug,vb &visited,const AdjMatrix<T> &adj,const AdjMatrix<T> &flow,int snk){
aug.pb(u);
visited[u]=true;

if(u==snk) return INF;

int n=adj.size();
rep(v,n)if(!visited[v] && adj[u][v]-flow[u][v]>0){
T ans=findAugpath(v,aug,visited,adj,flow,snk);
if(ans>0) return min(ans,adj[u][v]-flow[u][v]);
}

aug.pop_back();
return 0;
}

// BFS バージョンで使うサブルーチン
vi makepath(const vi &parent,int s,int g){
vi path;
for(int u=g;;u=parent[u]){
path.pb(u);
if(u==s) break;
}
reverse(path.begin(),path.end());
return path;
}

// BFS バージョン
template<class T>
T findAugpath(int src,vi &aug,vb &visited,const AdjMatrix<T> &adj,const AdjMatrix<T> &flow,int snk){
int n=adj.size();
vi parent(n,-1);
visited[src]=true;

T ans=0;
queue< pair<T,int> > qu; qu.push(mp(INF,src));
while(!qu.empty()){
pair<T,int> a=qu.front(); qu.pop();
int u=a.second;
T water=a.first;
if(u==snk){ ans=water; break; }
rep(v,n)if(!visited[v] && adj[u][v]-flow[u][v]>0){
qu.push(mp(min(water,adj[u][v]-flow[u][v]),v));
parent[v]=u;
visited[v]=true;
}
}

if(ans>0) aug=makepath(parent,src,snk);

return ans;
}

template<class T>
T FordFulkerson(const AdjMatrix<T> &adj,int src,int snk){
int n=adj.size();
AdjMatrix<T> flow(n,vector<T>(n));

T ans=0;
vi aug;
vb visited(n);
while(1){
aug.clear();
visited.assign(n,false);
T water=findAugpath(src,aug,visited,adj,flow,snk);
if(water==0) break;

rep(i,aug.size()-1){
flow[aug[i]][aug[i+1]]+=water;
flow[aug[i+1]][aug[i]]-=water;
}
ans+=water;
}

return ans;
}


Capacity Scaling による実装
// DFS バージョン
template
<class T>
T findAugpath(int u,vi &aug,vb &visited,const AdjMatrix<T> &adj,const AdjMatrix<T> &flow,int snk,T delta){
aug.pb(u);
visited[u]=true;

if(u==snk) return INF;

int n=adj.size();
rep(v,n)if(!visited[v] && adj[u][v]-flow[u][v]>=delta){
T ans=findAugpath(v,aug,visited,adj,flow,snk,delta);
if(ans>0) return min(ans,adj[u][v]-flow[u][v]);
}

aug.pop_back();
return 0;
}

// makepath は上と同じ

// BFS バージョン
template<class T>
T findAugpath(int src,vi &aug,vb &visited,const AdjMatrix<T> &adj,const AdjMatrix<T> &flow,int snk,T delta){
int n=adj.size();
vi parent(n,-1);
visited[src]=true;

T ans=0;
queue< pair<T,int> > qu; qu.push(mp(INF,src));
while(!qu.empty()){
pair<T,int> a=qu.front(); qu.pop();
int u=a.second;
T water=a.first;
if(u==snk){ ans=water; break; }
rep(v,n)if(!visited[v] && adj[u][v]-flow[u][v]>=delta){
qu.push(mp(min(water,adj[u][v]-flow[u][v]),v));
parent[v]=u;
visited[v]=true;
}
}

if(ans>0) aug=makepath(parent,src,snk);

return ans;
}

template<class T>
T ScalingMaxFlow(const AdjMatrix<T> &adj,int src,int snk){
int n=adj.size();
AdjMatrix<T> flow(n,vector<T>(n));

T deltamax=0;
rep(u,n)rep(v,n) deltamax=max(deltamax,adj[u][v]);
for(T i=1;;i<<=1)if(deltamax<2*i){ deltamax=i; break; }

T ans=0;
vi aug;
vb visited(n);
for(T delta=deltamax;delta>=1;delta>>=1){
while(1){
aug.clear();
visited.assign(n,false);
T water=findAugpath(src,aug,visited,adj,flow,snk,delta);
if(water==0) break;

rep(i,aug.size()-1){
flow[aug[i]][aug[i+1]]+=water;
flow[aug[i+1]][aug[i]]-=water;
}
ans+=water;
}
}

return ans;
}


速度計測

POJ 1459 で速度を比較してみた。
問題のサイズは
0 ≦ n ≦ 100
0 ≦ m ≦ 10000
辺容量の最大値 ≦ 10000

制限時間は 2 秒

結果は
・ DFS : Time Limit Exceeded
・ BFS : 922ms
・ DFS by Capacity Scaling : 719 ms
・ BFS by Capacity Scaling : 750 ms

おおよそ予想通りの結果。
静的配列で実装すれば、あと 200 ms くらいは縮まると思う。



全然まとまらなかった。
スポンサーサイト

コメントの投稿

非公開コメント

プロフィール

fura2

Author : fura2
数学・コンピュータを中心に、考えたこと・やったことを書いていきます。

誤植等を含め、間違いはご指摘いただければ幸いです。

FC2カウンター
検索フォーム
最新記事
最新コメント
最新トラックバック
月別アーカイブ
リンク
RSSリンクの表示
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。