【網路流】最大流總結

2012-12-25

剛才把預流推進裡的隊列改成優先隊列【手寫堆，不會STL給跪了】，為啥速度變得更慢了QAQ測速用的POJ3469，普通預流推進2.7s，換成優先隊列變成3.2s了不科學，求解答。

2012-12-28

改成STL的優先隊列了，不過整體上的速度還是比SAP要慢，難道還有可以最佳化的地方麼~求解

2013-01-02

我的sap是單路增廣~今天看了一下多路增廣但是不是很明白，目前單路增廣的效率我已經很滿意了，多路增廣慢慢看看有沒有必要用吧

2013-01-05

補充一下dinic的模板~不過寫遞迴要小心棧溢出什麼的

這兩天把預流推進看了一下，今天完成了relabel to front的模板，不過還沒有用優先隊列實現最高標號推進什麼的。

這裡對求最大流的演算法做一個總結。

最大流演算法有增廣路演算法和預流推進演算法兩個常用的系列。理論時間複雜度一般是O（N*M^2）或O（N^2*M），但是最高標號預流推進的理論複雜度逆天到O（N^2*sqrt（M）），但是我沒有實現也沒有實測到底有多快。不過網路流的理論複雜度我完全不知道是怎麼算的，實際速度也比理論要快上好多的樣子。

先介紹增廣路系列：

增廣路系列的各種演算法在於找增廣路的方式不同，最早的FF據說是採用DFS找增廣路，然後導致理論複雜度無法計算什麼的。然後有了EK演算法把FF的DFS換成了BFS，複雜度變為O（N*M^2），EK演算法應該算的上是最簡單易懂易寫的最大流演算法了吧，前些天有智能專業的妹子說開了演算法課，然後說講了網路流TAT【我們學校竟然還有這麼神的專業】，然後說聽不懂讓我給講><當時就隨手寫了一個EK【我會說那是我第一次寫EK麼】。唔~扯遠了…………

然後就是可能大家最常用的dinic和sap了，這兩種演算法都是每次找最短增廣路，所謂最短增廣路，假設每條邊的長度都是1，我們使用一個dis數組，dis[i]表示點i到匯點的最短路，對於每一條邊u->v，只有滿足dis[u]==dis[v]+1的才是允許弧，也就是說每一次找到的增廣路上所有邊都滿足這個條件，當沒有增廣路時要重新計算dis，直到dis[s]>=n時演算法結束，注意到每次重新標號都會導致這個點的dis不變或者增大，因為最短增廣路裡包含的邊只會越來越多。黑書上說採用最短增廣路可以把每次找增廣路的平均時間複雜度從O（M）降到O（N），我表示沒有找過相關的證明，這裡就不對複雜度的計算做深入討論了。dinic與sap的區別在於，dinic每次重新標號是直接對整個圖進行BFS，sap則是每次對當前不能向下繼續增廣的點重新標號，因為我只用過沒加最佳化的dinic，所以不好對dinic和sap的速度做比較，不過只看實現的話sap比dinic還是好寫一些的。

因為sap是我最常用的，這裡介紹一下sap的幾個最佳化。

首先是gap最佳化，這個是對sap最佳化最大的，不過我沒有測過不加gap最佳化的sap的速度到底如何，所以我也說不清這個最佳化到底有多明顯。所謂gap最佳化是指gap[i]記錄dis為i有點的個數，如果出現gap[i]==0的話直接結束，這時的求出的就是最大流，因為如果出現gap[i]==0就說明層次網路裡出現了斷層，不存在最短增廣路。關於每個點的初始dis值，推薦在開始sap前先從匯點開始一次反向BFS，BFS過程的允許邊是流量為0的邊，也就是初始建圖時每條邊的反向邊，雖然也可以將所有點的dis都設為0然後直接跑，但是其實這個BFS在多數情況下還是很有用的，尤其是在點比較多的情況下，這一次bfs可以明顯減少sap過程中重新標號的次數，HDU4280那道平面圖網路流【正解不是sap】我後來就直接用加了BFS的sap過掉的，不加BFS就逾時……還有一個常數層級的當前弧最佳化，看很多人說這個最佳化的效果不明顯，不過我還是說明一下，當前弧最佳化的原理是記錄每一個點在當前情況下可以用於找增廣路的第一條弧，然後從這條弧開始增廣，當前弧之前的所有弧都不會出現在增廣路中，直到當前點被重新標號為止。

然後是預流推進：

relabel to front

上面這個連結是DD神寫的relabel to front相關。

預流推進也需要像sap一樣先計算出dis，注意dis[源點]=n。首先把與源點相連的連都置為滿流，對於每一個點用一個ef數組記錄這個點積攢的流量，對於非源匯點ef[i]最終應該為0，這樣用一個隊列維護ef[i]！=0的點，對於邊u->v，如果滿足dis[u]==dis[v]+1那麼將點u積攢的流量盡量壓入點v，壓入流量=min(ef[u],c[u][v])【c[u][v]是這條邊的容量】，如果ef[v]！=0且點v不是源點和匯點時將點v排入佇列，如果找不到這樣的邊就對點u重新標號，直到ef[u]==0，如果某個點的積攢流量無法向下壓時最終都會回到源點上，這樣當所有隊列為空白時ef[匯點]記錄的就是最大流。

預流推進同樣有gap最佳化，當gap[i]==0時將所有dis[u]>i 且n+1的點的dis置為n+1【n為圖中點的數量】，這個最佳化是非常明顯的，對於POJ3469那個題不加gap最佳化直接TLE，加了gap最佳化以後可以跑到2.7秒。如果某一層出現斷層，那麼比這一層更高的點的流量必然要都要被退回源點，不可能再繼續向下推進，直接標記為n+1等待回推。

那個理論複雜度只有O（N^2*sqrt（M））的最高標號預流推進就是用一個優先隊列每次取dis最大的點進行推進，完全不能理解為啥複雜度這麼低><不會算複雜度的給跪了。

高標預流我按網上流傳的一份吉大的模板寫了一下結果比普通的預流推進還慢，果然還是要寫優先隊列呀啊啊啊。

這裡先送上qinz的神SAP模板，跑起來相當快Qinz's SAP

接下來是我的模板：

先是SAP：【加滿了所有最佳化】

const int N=100010;const int M=400010;const int inf=1<<30;struct Edge{    int v,next,w;}e[M];int head[N],cnt;int n,m,ans;int s,t;int pre[N],cur[N],dis[N],gap[N];int q[N],open,tail;void addedge(int u,int v,int w){    e[cnt].v=v;    e[cnt].w=w;    e[cnt].next=head[u];    head[u]=cnt++;    e[cnt].v=u;    e[cnt].w=0;    e[cnt].next=head[v];    head[v]=cnt++;}void BFS(){int i,u,v;memset(gap,0,sizeof(gap));memset(dis,-1,sizeof(dis));open=tail=0;q[open]=t;dis[t]=0;while (open<=tail){u=q[open++];for(i=head[u];i!=-1;i=e[i].next){v=e[i].v;if(e[i].w!=0 || dis[v]!=-1) continue;q[++tail]=v;++gap[dis[v]=dis[u]+1];}}}int sap(int n){    int i,v,u,flow=0,aug=inf;    int flag;    BFS();    gap[0]=1;    for (i=1;i<=n;i++) cur[i]=head[i];    u=pre[s]=s;while (dis[s]<n){   flag=0;        for (int j=cur[u];j!=-1;j=e[j].next) {   v=e[j].v;            if (e[j].w>0&&dis[u]==dis[v]+1)   {                flag=1;    if (e[j].w<aug) aug=e[j].w;     pre[v]=u;      u=v;                if (u==t)   {                    flow+=aug;                    while (u!=s)   {  u=pre[u];   e[cur[u]].w-=aug;    e[cur[u]^1].w+=aug;    }                    aug=inf;                }                break;            }            cur[u]=e[j].next;        }        if (flag) continue;        int mindis=n;        for (int j=head[u];j!=-1;j=e[j].next)    {            v=e[j].v;            if (e[j].w>0&&mindis>dis[v])   {   mindis=dis[v];    cur[u]=j;}        }        if (--gap[dis[u]]==0) break;    ++gap[dis[u]=mindis+1];     u=pre[u];    }   return flow;}

然後是今天寫好的預流推進：

#define N 210#define M 40010#define INF (1<<30)int n,m,s,t;int head[N],cnt,cur[N];int dis[N],ef[N],gap[N];struct edge{int v,w,next;}e[M];queue<int> q;void addedge(int u,int v,int w){e[cnt]=(edge){v,w,head[u]}; head[u]=cnt++;e[cnt]=(edge){u,0,head[v]}; head[v]=cnt++;}void push(int u,int i){int aug=min(e[i].w,ef[u]),v=e[i].v;if(v!=s && v!=t && ef[v]==0) q.push(v);e[i].w-=aug;e[i^1].w+=aug;ef[u]-=aug;ef[v]+=aug;}void relabel(int u){int mindis=INF,nl=dis[u];--gap[dis[u]];for(int i=head[u];i!=-1;i=e[i].next)if(e[i].w!=0 && mindis>dis[e[i].v]){cur[u]=i;mindis=dis[e[i].v];}++gap[dis[u]=mindis+1];if(gap[nl]!=0 || nl>=n+1) return;for(int i=1;i<=n;i++)if(dis[i]>nl && dis[i]<=n && i!=s){--gap[dis[i]];++gap[n+1];dis[i]=n+1;}}void discharge(int u){while(ef[u]>0){if(cur[u]==-1) relabel(u);int v=e[cur[u]].v,w=e[cur[u]].w;if(w>0 && dis[u]==dis[v]+1) push(u,cur[u]);cur[u]=e[cur[u]].next;}}void BFS(){q.push(t); dis[t]=0;while(! q.empty()){int u=q.front();q.pop();for(int i=head[u];i!=-1;i=e[i].next)if(e[i].w==0 && dis[e[i].v]==-1 && e[i].v!=s){dis[e[i].v]=dis[u]+1;q.push(e[i].v);}}}int maxflow(){int i;for(i=1;i<=n;i++)cur[i]=head[i],dis[i]=-1,ef[i]=gap[i]=0;BFS();dis[s]=n;ef[s]=INF;for(i=1;i<=n;i++)if(dis[i]==-1) dis[i]=n;for(i=1;i<=n;i++)gap[dis[i]]++;for(i=head[s];i!=-1;i=e[i].next) push(s,i);while(! q.empty()){int u=q.front();q.pop();discharge(u);}return ef[t];}

dinic

#pragma comment(linker, "/STACK:32000000")#define N 100010#define M 200010#define INF (1<<30)int n,m,s,t,ans;int head[N],cnt,dis[N];struct edge{    int v,w,next;    }e[M];void addedge(int u,int v,int w){   e[cnt].v=v; e[cnt].w=w; e[cnt].next=head[u]; head[u]=cnt++;    e[cnt].v=u; e[cnt].w=w; e[cnt].next=head[v]; head[v]=cnt++;}queue<int> q;int bfs(){    memset(dis,0,sizeof(dis));        dis[s]=1; q.push(s);    while (! q.empty())    {        int u=q.front(),v;        q.pop();        for(int i=head[u];i!=-1;i=e[i].next)        if(e[i].w && dis[v=e[i].v]==0)  {    dis[v]=dis[u]+1;    q.push(v);    }    }    return dis[t];}int dfs(int s,int limit){    if(s==t) return limit;    int v,tmp,cost=0;    for(int i=head[s];i!=-1;i=e[i].next)    if(e[i].w && dis[s]==dis[v=e[i].v]-1)    {        tmp=dfs(v,min(limit-cost,e[i].w));        if(tmp>0)        {    e[i].w-=tmp;    e[i^1].w+=tmp;    cost+=tmp;    if(limit==cost) break;        }else dis[v]=-1;    }    return cost;}while (bfs()) ans0+=dfs(s,INF);