ネットワークフロー-最大フロー問題ISAPアルゴリズム解釈(Renfei Song's Blogから転送)

ネットワークフロー-最大フロー問題ISAPアルゴリズム解釈
August 7, 2013 / 
プログラミングガイド
ISAPは最大ストリームを求めるアルゴリズムの1つであり,実行時間とプログラムの複雑さの関係をよくバランスさせるため,非常によく用いられる.
に約束
隣接表を用いて図を表し,表現方法は文章の重み付き最短Dijkstra,SPFA,Bellman−Ford,ASP,Floyd−Warshallアルゴリズム解析または二分図の最大整合,完全整合,ハンガリーアルゴリズムの先頭を見ることができる(コードを重複しない).以下では、図のソースポイントをsと表す
、シンクタンク(sink)は
t,現在のノードは
u .図面を作成する場合は、双方向エッジを作成する必要があります(逆方向のエッジ容量を
0
)を使用して、アルゴリズムが正しいことを保証できます.
導入
最大ストリーム問題を解く比較的容易な方法は,残量ネットワーク(residual network)のたびにsから
tの経路は,その後,このような経路が存在しないまで拡張される.これが一般的な拡張アルゴリズム(labeling algorithm)である.この改良されていない貪欲なアルゴリズムが正しいことを証明することができる.最大ストリームが
f,それではその運転時間は
O( f⋅∣E∣) .しかし、この運転時間は、最大ストリームと
f
関連する.
毎回残量ネットワークの最短拡張路に沿って拡張すれば,実行時間をOに減らすことができることが分かった(∣E∣2⋅V∣)
.これが最短拡張アルゴリズムです.ISAPアルゴリズムは最短拡張アルゴリズムの改良である.実は、ISAPの意味は「改良された最短拡張路」(Improved Shortest Augmenting Path)です.
ちなみに,上述したすべてのアルゴリズムは根本的に拡張路法(Ford−Fulkerson method)に属する.それに対応するのが有名なプレフロー推進法(Preflow-push method).ここで、最高符号プリフロー推進アルゴリズム(Highest-label preflow-push algorithm)の複雑さはO(
.複雑度は拡張路法よりはるかに進歩しているが,プリフロー推進アルゴリズムの複雑度の上界は比較的きついため,差が大きくない場合がある.
アルゴリズム解釈
要約すると、ISAPアルゴリズムは最短の拡張路を探し続け、見つけてから拡張することである.死の道に出会ったらsを発見するまでretreatし、
tは連通せず,アルゴリズムは終了する.最短ルートを探すことは本質的に最短ルートの問題であるため,BFSの考え方を採用する.具体的には、配列を使用します.
d,各ノードから合流点までを記録する
tの最短距離.検索するときは、満足に沿ってのみ
d[u]=d[v]+1のエッジ
u→v
(このようなエッジを許容アークと呼ぶ)歩く.明らかに、このように出てきたのは一番短いに違いない.
原図には2つのサブマップが存在し,1つは残量ネットワークであり,1つはアーク組成を許容する図である.残量ネットワーク保証は広がり、アーク保証が最も短絡することを可能にする(時間界が優れている).従って,拡張路を探す過程で,残量ネットワークにおいて許容アークに沿って探し続けた.したがって、許容アークは、原図ではなく残量ネットワークに属するべきである.すなわち,許容アークに沿って,原図ではなく残量ネットワークの最短経路を歩いた.残量ネットワークに沿って拡張路を見つけると、拡張後、残量ネットワークは必ず変化する(少なくとも1つのエッジが少ない)ため、アークのd配列を対応する更新を許可することにした(ちなみにDinicのやり方は、拡張するたびに再計算することである
d配列).しかし、ISAPの「改善」の一つは、すぐに更新する必要がないことです.
d
配列.これは、1つのエッジを削除すると、経路がより長くなる可能性があり、拡張前の残量ネットワークに別の最短回路が存在し、拡張後の残量ネットワークに依然として存在する場合、この経路が最も短いことは間違いないからである.したがって、ISAPのアプローチは、デッドラインに遭遇するまで拡大し続け、retreat操作を実行することです.
ここまで言うと、retreat操作の主な任務はd配列を更新することであると推測される.では、どのように更新しますか?非常に簡単:ノードから
u隣接するエッジを探しても許容アークは見つからなかった.変数の再設定
m,令
mイコール残量ネットワークにおける
uのすべての隣接点の
d配列の最小値、そして
d[u]は
m+1でいいです.これは、retreatの一環に入って残量ネットワークを説明するためです.
u和
tは(時代遅れ)の許容アークで接続できなくなり、
u和
tは実際に残量ネットワークの中で最も短絡的な長さはいくらですか?(これはまさに
dの定義!)明らかに残量ネットワークでは
uのすべての隣接点と
tの距離加算
1の最小ケース.特に、残量ネットワークでは
u隣接点がまったくない.もしそうなら、ただ
d[u]は比較的大きな数に設定すればよいので、任意の点から
uのエッジは残量ネットワークから除外される.(厳密には等しい以上であれば
ゞ∣V∣即可.ゞ最短路は必ず無環であるため、任意の経路長が最大
∣V∣−1).修正後、検討中のノードを
u
さっき歩いた道に沿って一歩下がって、検索を続けてください.
ここまで言うと、ISAPアルゴリズムのフレームワークの内容は終わります.コード自体については、いくつかの最適化と実装のテクニックについて説明する必要があります.

アルゴリズムの実行前にBFSでd

を初期化する必要がある.
配列、メソッドは
tから
s
逆方向に進む.アルゴリズム本体は「現在のノード」を維持する必要があります
u,このノードの前進,retreatなどの操作を実行する.パスを記録する方法は非常に簡単で、配列を宣言します.
p,令
p[i]拡大路到達ノードに等しい
iのエッジのシーケンス番号(頂点から頂点までの頂点を見つけることができます)
i).経路が必要なときは逆方向に追跡すればいいです.判定残量ネットワークにおける
s,tが連通しない条件は,
d[s]≥∣V∣ .なぜなら
s,tが連通しない場合,最終残量ネットワークにおいて
sには隣接点はありません
sのretreatは上記の条件の成立を招く.GAP最適化.GAP最適化はプログラムを早期に終了させることができ,多くの場合,スピードアップが非常に顕著である(100倍以上に達する).GAP最適化とは、retreatコーナーに入ると、
u,t間の連通性は消失するが,
uは最後の和
t距離
d[u](更新前)の点、この場合を説明
s,tもつながっていない.なぜなら
u,tはもうつながっていませんが、結局私たちが歩いているのは最短で、他の点はこの時着きます.
tの距離はきっと大きい
d[u](更新前)なので、他のポイントは
t,必ず1つの和を通らなければならない
t距離は
d[u]

(更新前)のポイントです.GAP最適化の実現は非常に簡単で,1つの配列で記録し,適切なタイミングで判断し,ループから飛び出せばよい.

のもう一つの最適化は、1つの配列で1つの点がどの隣接エッジを試みたかを保存することである.拡張を探すプロセスは実際にはBFSプロセスに似ているため,以前に処理された隣接エッジは再処理する必要はない(残量ネットワーク内のエッジはますます少なくなるだけである).具体的な実装方法はコードを直接見ればよく,非常に理解しやすい.注意すべき点は、前回処理された隣接辺の次の条からではなく、次回は前回処理された隣接辺から処理を継続すべきである.

最後に拡張過程についてお話しします.拡張プロセスは非常に簡単で、拡張パスが成功した(現在のノードがtに処理されている)ことを探した後、あなたが記録したパスに沿って歩いて、道中の最小残量を記録して、それから
sから
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 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;
}

私のTemplate:

const int maxn=1000;
const int INF=0x3f3f3f3f;
int s,t,first[maxn],nxt[maxn];
struct Edge{
	int u,v,cap,flow;
}e[maxn];
bool vis[maxn];
int q[maxn],d[maxn],p[maxn],num[maxn],cur[maxn];
void bfs()//  bfs     。  d[]  ;
{
    memset(vis,false,sizeof(vis));
    int head=0,tail=1;
    q[0]=t;
    d[t]=0;
    vis[t]=true;
    while(head!=tail){
        int now=q[head];head++;
        for(int i=first[now];i;i=nxt[i])
            if(!vis[e[i].u]&&e[i].cap>e[i].flow){
                vis[e[i].u]=true;
                d[e[i].u]=d[now]+1;
                q[tail++]=e[i].u;
            }
    }
}
int Agument()
{
    int x=t,a=INF;//        p        ;
    while(x!=s){
        a=min(a,e[p[x]].cap-e[p[x]].flow);//      ;
        x=e[p[x]].u;
    }
    x=t;//    ;
    while(x!=s){
        e[p[x]].flow+=a;
        e[p[x]^1].flow-=a;
        x=e[p[x]].u;
    }
    return a;
}
int ISAP()
{
    int flow=0;//flow      int   ;
    bfs();//    d[]  ;
    memset(num,0,sizeof(num));//num     ，   i         ，  gap  ;
    for(int i=1;i<=n;i++)num[d[i]]++;
    int x=s;
    memset(cur,0,sizeof(cur));//currently，         ;
    while(d[s]e[i].flow&&d[x]==d[e[i].v]+1){
                ok=true;
                p[e[i].v]=i;
                cur[x]=i;
                x=e[i].v;
                break;
            }
        if(!ok){//retreat  ：  d[]  
            int mn=n-1;
            for(int i=first[x];i;i=nxt[i])
                if(e[i].cap>e[i].flow)mn=min(mn,d[e[i].v]);
            if(--num[d[x]]==0)break;//gap   ：  x         distance d[x]  ，     s、t      ，     ;
            num[d[x]=mn+1]++;//      ;
            cur[x]=0;
            if(x!=s)x=e[p[x]].u;
        }
    }
    return flow;
}
void Link()//     ;
{
    memset(nxt,0,sizeof(nxt));
    memset(first,0,sizeof(first));
    for(int a,b,c;m;m--){
        scanf("%d%d%d",&a,&b,&c);
        e[++ecnt].u=a,e[ecnt].v=b,e[ecnt].cap=c,e[ecnt].flow=0;
        nxt[ecnt]=first[a],first[a]=ecnt;
        e[++ecnt].u=b,e[ecnt].v=a,e[ecnt].cap=0,e[ecnt].flow=0;
        nxt[ecnt]=first[b],first[b]=ecnt;
    }
}