[理性愉悦]最大流算法学习

本文总结了一些OI中常用和不常用的最大流算法,并对运行效率进行简单测试。主要是闲着无聊没事干

Tip:除了Dinic以外其他算法我都是现学现写的,理解/代码上可能有各种错误,欢迎拍砖qwq

0. Edmonds-Karp

由于这种方法实在太慢了,这里不提。

不过这个做法告诉了我们一个重要道理:网络流找增广路要加反向边。

1. Dinic

这个经典网络流算法相信大家都耳熟能详。

我们在残余网络上bfs,把访问到的点按层标号,接着将跨层的边进行dfs增广。(当然一句话讲的也不是很清楚,不会的自行百度吧)

正确的dinic分层姿势是从T开始倒着分层,从S开始正着分层复杂度似乎不变,但是从T开始跑的快一点。

Dinic主要有两个优化,一个是当前弧优化(即不重复访问一条边),另一个是当一个点无法增广时就放弃这个点。还有一个优化比较玄学,因为bfs很慢,所以我们随便搜出一条路就返回,复杂度全看脸,详见代码中的注释。

#define SZ 666666
typedef long long ll;
int fst[SZ],nxt[SZ],vb[SZ],cap[SZ],ff[SZ],M=1,N=0;
#define ad_de ad_dl
void _ad_dl(int a,int b,int c) {++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;cap[M]=c;}
void ad_dl(int a,int b,int c)
{_ad_dl(a,b,c); _ad_dl(b,a,0);}
int S,T,q[SZ],d[SZ];
bool bfs()
{
    memset(d,-1,4*(N+1));
    d[T]=0; q[1]=T; int h=1,t=2;
    while(h!=t)
    {
        int cur=q[h++];
        for(int e=fst[cur];e;e=nxt[e])
        {
            int b=vb[e];
            if(d[b]==-1&&cap[e^1])
            {
                d[b]=d[cur]+1, q[t++]=b;
                if(b==S) return 1; //这是玄学 
            }
        }
    }
    return d[S]!=-1;
}
int dfs(int x,int f)
{
    if(f<=0) return 0;
    if(x==T) return f;
    int ca=0;
    for(int& e=fst[x];e;e=nxt[e]) //当前弧优化 
    {
        int b=vb[e];
        if(d[b]+1!=d[x]) continue;
        int w=dfs(b,(cap[e]<f-ca)?cap[e]:(f-ca));
        cap[e]-=w; cap[e^1]+=w; ca+=w;
        if(ca==f) break;
    }
    if(!ca) d[x]=-1; //另一个优化 
    return ca;
}
#define inf 2e9+3
int dinic()
{
    int ans=0;
    memcpy(ff,fst,4*(N+1));
    while(bfs())
    {
        //当前弧优化
        ans+=dfs(S,inf);
        memcpy(fst,ff,4*(N+1));
    }
    return ans;
}

理论时间复杂度为O(n^2m),据说寻找增广路时使用LCT可以优化到O(nmlog(n))image)。

2. Sap

考虑仍然像dinic一样把从T开始反着把图分层,这回我们增广时采用dfs,只进行单路增广,只走相邻层的边,找到一条增广路就增广返回起点。

这个做法听上去非常naive,但是可以加入两个优化。

Dinic的时候我们每次增广完一堆路径,这时候我们会重新跑一遍bfs,然后再增广一堆路径,而Sap不需要重新跑一遍bfs,而是当某个点找不到增广路时,把它到T的层数用它所有有流量的后继层数最小值+1来更新再返回,这样就避免了重新跑bfs。

此外还有一个叫gap的玄学优化,就是说如果从S到T中间有一个层一个点也没有,那么这个层就是一个“gap”,阻止了流量的传输,那么我们就可以直接退出bfs节省时间。

当然当前弧优化也是可以使用的。实际写的时候并不需要真的dfs,而是记录每个点的pre(反正增广也要用),直接回溯。

#define SZ 666666
int fst[SZ],nxt[SZ],vb[SZ],cap[SZ],ff[SZ],M=1,N=0;
#define ad_de ad_dl
void _ad_dl(int a,int b,int c) {++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;cap[M]=c;}
void ad_dl(int a,int b,int c)
{_ad_dl(a,b,c); _ad_dl(b,a,0);}
int S,T,q[SZ],d[SZ],gap[SZ];
int pre[SZ],low[SZ],pe[SZ];
int sap()
{
    for(int i=0;i<=N;++i) d[i]=N,ff[i]=fst[i];
    d[T]=0; q[1]=T; int h=1,t=2;
    while(h^t)
    {
        int x=q[h++];
        for(int e=fst[x];e;e=nxt[e])
        {
            int b=vb[e];
            if(d[b]==N&&cap[e^1])
                d[b]=d[x]+1, q[t++]=b;
        }
    }
    for(int i=1;i<=N;++i) ++gap[d[i]];
    int x=S,ans=0;low[S]=2e9;
    while(d[S]<N)
    {
        int ce=-1;
        for(int& e=fst[x];e;e=nxt[e]) //当前弧优化 
        {
            if(!cap[e]||d[vb[e]]+1!=d[x]) continue;
            ce=e; break;
        }
        if(ce==-1)
        {
            --gap[d[x]];
            if(!gap[d[x]]) break; //gap优化 
            d[x]=N; fst[x]=ff[x]; //回退当前弧 
            for(int e=fst[x];e;e=nxt[e])
                (cap[e]&&d[vb[e]]+1<d[x])?(d[x]=d[vb[e]]+1):0; //sap
            ++gap[d[x]];
            if(x!=S) x=pre[x]; //回溯 
        }
        else
        {
            int b=vb[ce]; pre[b]=x;pe[b]=ce;
            low[b]=(low[x]<cap[ce])?low[x]:cap[ce];
            x=b; if(x!=T) continue;
            //找到了增广路
            ans+=low[x];
            for(int p=x;p!=S;p=pre[p])
                cap[pe[p]]-=low[x],cap[pe[p]^1]+=low[x];
            x=S;
        }
    }
    return ans;
}

它的理论时间复杂度仍然为O(n^2m)

3. HLPP

http://blog.csdn.net/hechenghai/article/details/42024683

https://en.wikipedia.org/wiki/Push%E2%80%93relabel_maximum_flow_algorithm

这个算法全称叫最高标号预流推进,和Sap类似,仍然从T开始给每个点一个距离标号,无法增广时就更新标号。开始每个点的距离标号可以赋为0,也可以正常bfs。

既然叫“最高标号”预流推进算法,这个算法按标号从大到小维护了一个优先队列,表示要增广的点。

那么什么叫“预流推进”呢?“预流”就表示想往外流,每次考虑要增广一条边(a,b),那么a的预流就减掉这条边的流量,b的预流加上这条边的流量。

这个算法主体用bfs增广,增广一个点时,考虑每条通往下层的边,按上面的做法推流,推到了没进队的点就加入优先队列,当这个点的预流为0时它的历史使命就完成了。如果一个点想要推流却啥也没推出去,那么就更新它的标号。最后检查这个点的预流是否为0,没0继续加入队列。

注意S和T除了一开始S加入队列外不应该加入队列。开始将S的预流赋为∞,该推的流都推完以后T的预流就是答案。

当前弧优化在这个算法中似乎不适用。

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <string>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <algorithm>
#include <sstream>
#include <stack>
#include <iomanip>
using namespace std;
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
typedef double ld;
typedef vector<int> vi;
#define fi first
#define se second
#define fe first
#define FO(x) {freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);}
#define Edg int M=0,fst[SZ],vb[SZ],nxt[SZ];void ad_de(int a,int b){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;}void adde(int a,int b){ad_de(a,b);ad_de(b,a);}
#define Edgc int M=0,fst[SZ],vb[SZ],nxt[SZ],vc[SZ];void ad_de(int a,int b,int c){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;vc[M]=c;}void adde(int a,int b,int c){ad_de(a,b,c);ad_de(b,a,c);}
#define es(x,e) (int e=fst[x];e;e=nxt[e])
#define esb(x,e,b) (int e=fst[x],b=vb[e];e;e=nxt[e],b=vb[e])
#define VIZ {printf("digraph G{\n"); for(int i=1;i<=n;i++) for es(i,e) printf("%d->%d;\n",i,vb[e]); puts("}");}
#define VIZ2 {printf("graph G{\n"); for(int i=1;i<=n;i++) for es(i,e) if(vb[e]>=i)printf("%d--%d;\n",i,vb[e]); puts("}");}
#define SZ 666666
int fst[SZ],nxt[SZ],vb[SZ],cap[SZ],M=1,N=0;
#define ad_de ad_dl
void _ad_dl(int a,int b,int c) {++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;cap[M]=c;}
void ad_dl(int a,int b,int c)
{_ad_dl(a,b,c); _ad_dl(b,a,0);}
int S,T,q[SZ],d[SZ],ex[SZ];
bool inq[SZ]; priority_queue<pii> pq;
#define gq(x) ((x)!=S&&(x)!=T&&ex[x]&&!inq[x])
int hlpp()
{
    for(int i=0;i<=N;++i) d[i]=N+2;
    d[T]=0; q[1]=T; int h=1,t=2;
    while(h^t)
    {
        int x=q[h++];
        for(int e=fst[x];e;e=nxt[e])
        {
            int b=vb[e];
            if(d[b]==N+2&&cap[e^1])
                d[b]=d[x]+1, q[t++]=b;
        }
    }
    d[S]=N+2;
    pq.push(pii(d[S],S)); inq[S]=1; ex[S]=2e9;
    while(!pq.empty())
    {
        int x=pq.top().se; pq.pop(); inq[x]=0;
        bool p=0;
        for(int e=fst[x];e;e=nxt[e])
        {
            int b=vb[e];
            if(!cap[e]||(x!=S&&d[b]+1!=d[x])) continue;
            p=1; int f=(cap[e]<ex[x])?cap[e]:ex[x];
            ex[x]-=f; ex[b]+=f;
            cap[e]-=f; cap[e^1]+=f;
            if(gq(b)) inq[b]=1,pq.push(pii(d[b],b));
            if(!ex[x]) break;
        }
        if(!p)
        {
            d[x]=N+2;
            for(int e=fst[x];e;e=nxt[e])
                (cap[e]&&d[vb[e]]+1<d[x])?(d[x]=d[vb[e]]+1):0;
        }
        if(gq(x)) inq[x]=1,pq.push(pii(d[x],x));
    }
    return ex[T];
}
namespace FF
{
char ch,B[1<<15],*S=B,*T=B;
#define getc() (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++)
#define isd(c) (c>='0'&&c<='9')
int aa,bb;int F(){
    while(ch=getc(),!isd(ch)&&ch!='-');ch=='-'?aa=bb=0:(aa=ch-'0',bb=1);
    while(ch=getc(),isd(ch))aa=aa*10+ch-'0';return bb?aa:-aa;
}
}
#define gi FF::F()
int main()
{
    N=gi; int m=gi;
    S=1; T=N;
    for(int i=1;i<=m;i++)
    {
        int x=gi,y=gi,z=gi;
        ad_dl(x,y,z);
    }
    printf("%d\n",hlpp());
}

UPD:注意!这份代码是错误的,原因我现在还不清楚,反正是错的,下面有关HLPP的测试结果请忽略。

理论时间复杂度比上面两个优秀!为O(n^2\sqrt{m}),这个界是紧的。

4. Goldberg-Rao

http://www.cs.princeton.edu/courses/archive/fall06/cos528/handouts/Goldberg-Rao.pdf

http://www.cs.tau.ac.il/~haimk/adv-alg-2013/goldbergrao.pdf

https://people.orie.cornell.edu/dpw/orie633/LectureNotes/lecture10.pdf

首先我们需要知道,dinic的复杂度十分玄学,在容量为1的图上跑的飞快,经过一番严谨的证明可以证出上界是O(min(n^{2/3},m^{1/2})m)(详见维基)。

我们考虑将dinic算法改造一番将它在容量大的时候也跑的很快。

在开始之前,我们设\Lambda=min(n^{2/3},m^{1/2}),最多能增广的量不超过F,\Delta=\lceil{F}/\Lambda\rceil

考虑形式化地定义dinic算法,我们定义一条边u->v是好的,当且仅当d(v)=d(u)+l,其中l表示这条边的长度。在原始dinic算法中,l恒为1,但在这个算法中,l可能为0或1。d的定义同理,形式化地,d(T)=0,对于每个点v,d(v)\leq{d(w)+l(v~to~w)},并且等号能取到。

我们定义一条边为大的,当且仅当它的残余容量\geq3\Delta;一条边为特别的,这个条件比较多,有三个条件:残余容量在[2\Delta,3\Delta),两端点距离标号相等,它的反向边残余容量\geq3\Delta

那么我们现在可以定义一条边的l了,定义十分简单,如果这条边是大的或者特别的,l为0,否则为1。

仔细一想你会觉得有些不对,l依赖“特别”来定义,“特别”依赖距离标号来定义,而距离标号却又依赖l来计算...但是事实上我们可以先忽略特别的边,即只考虑小的边将其l设为0,这样计算出距离标号之后再重新计算一遍l,因为特别的边显然不会改变距离标号。

然后我们开始先定一个F的上界,为了简化我们设U为每条边的最大容量,那么就开始设F=mU

算法流程大致如下:

初始化F
while F>=1:
  计算△
  重复5Λ次操作:
    初始化距离标号和边的l值
    在好的边上找一个堵塞的或者流量为△的流
    增广
  F/=2

“堵塞”的流就表示是增广出来满的一个流。是不是看起来十分简单呢?

我们还遗漏了一个重要细节。在正常dinic中,l都为1,增广轻松愉快。但在这个算法中l可能等于0,那么有可能出现流量循环的情况,因为同一层的点间可能会有流量。

流量循环?Tarjan!我们在好的边上跑一遍tarjan,把每个强连通分量缩成一个点,那么就不会有什么循环的情况出现了!

那么缩完点之后的流量怎么办?我们在强连通分量里随便找个点r当做根,同样在好边上dfs两遍搞出一棵In Tree,一棵Out Tree。盗图一张:

image

如果有了流量,我们通过In Tree上的边流给r,再从Out Tree上的边流出去,这样就搞完啦!

不过还有一个小优化,考虑每次F/=2显得非常蠢,我们考虑想办法乱搞求出残余流量的一个上界,进行估价,每次F=min(F/2,计算出的上界)。

我们考虑最小割=最大流,那么考虑算出一个比较小的割。考虑在原图上bfs求出距离标号,然后每次考虑对于每个k,把距离标号为k和k+1的点之间的边全部割掉,这显然是一个合法的割,虽然不是很小,我们就用这个进行估价。同时,开始的F也可以设为这个值。

当前弧优化还是可以继续使用的。另外如果找到的堵塞流流量为0,可以直接break。

正确性和复杂度证明留作习题

注意下面这份代码可能有各种各样的错误,数据中的边权都是随机的,所以增广一次就够了...更强的数据我并不会造,就先这样好了...

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include <stdlib.h>
#include <string>
#include <bitset>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <algorithm>
#include <sstream>
#include <stack>
#include <iomanip>
using namespace std;
#define pb push_back
#define mp make_pair
typedef pair<int,int> pii;
typedef long long ll;
typedef double ld;
typedef vector<int> vi;
#define fi first
#define se second
#define fe first
#define FO(x) {freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);}
#define Edg int M=0,fst[SZ],vb[SZ],nxt[SZ];void ad_de(int a,int b){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;}void adde(int a,int b){ad_de(a,b);ad_de(b,a);}
#define Edgc int M=0,fst[SZ],vb[SZ],nxt[SZ],vc[SZ];void ad_de(int a,int b,int c){++M;nxt[M]=fst[a];fst[a]=M;vb[M]=b;vc[M]=c;}void adde(int a,int b,int c){ad_de(a,b,c);ad_de(b,a,c);}
#define es(x,e) (int e=fst[x];e;e=nxt[e])
#define esb(x,e,b) (int e=fst[x],b=vb[e];e;e=nxt[e],b=vb[e])
#define VIZ {printf("digraph G{\n"); for(int i=1;i<=n;i++) for es(i,e) printf("%d->%d;\n",i,vb[e]); puts("}");}
#define VIZ2 {printf("graph G{\n"); for(int i=1;i<=n;i++) for es(i,e) if(vb[e]>=i)printf("%d--%d;\n",i,vb[e]); puts("}");}
#define SZ 666666
#define FEdge \
int M=1,fst[SZ],vb[SZ],cap[SZ],nxt[SZ],va[SZ];\
void ad_de_(int a,int b,int c)\
{++M; nxt[M]=fst[a]; vb[M]=b; cap[M]=c; fst[a]=M; va[M]=a;}\
void ad_de(int a,int b,int c)\
{ad_de_(a,b,c); ad_de_(b,a,0);}
#define GetD(l) \
int qs[SZ];\
bool inq[SZ];\
void spfa()\
{\
    int h=0,t=1;qs[0]=T;\
    for(int i=1;i<=n;i++) d[i]=1e9;\
    d[T]=0; inq[T]=1;\
    while(h^t)\
    {\
        int x=qs[h++];\
        for(int e=fst[x];e;e=nxt[e])\
        {\
            if(cap[e^1]);else continue;\
            int y=d[x]+l(e^1),b=vb[e];\
            if(d[b]<=y) continue;\
            d[b]=y; if(inq[b]) continue;\
            inq[b]=1; qs[t++]=b;\
        }\
        inq[x]=0;\
    }\
}
int n,m,S,T;
namespace dn
{
int T,d[SZ],M=1,fst[SZ],vb[SZ],cap[SZ],nxt[SZ],l[SZ];
int addedge(int a,int b,int c,int g)
{++M; nxt[M]=fst[a]; vb[M]=b; cap[M]=c; fst[a]=M; l[M]=g; return M;}
int dfs(int x,int f)
{
    if(f<=0) return 0;
    if(x==T) return f;
    int ca=0;
    for(int& e=fst[x];e;e=nxt[e])
    {
        int b=vb[e];
        if(d[b]+l[e]!=d[x]) continue;
        int w=dfs(b,min(cap[e],f-ca));
        cap[e]-=w; cap[e^1]+=w; ca+=w;
        if(ca==f) break;
    }
    if(!ca) d[x]=-233;
    return ca;
}
}
namespace ng
{
FEdge int d[SZ];
namespace og
{
int d[SZ],sm[SZ]; GetD(1+0*)
int gj()
{
    spfa();
    int mx=d[S];
    if(mx==1e9) return 0;
    for(int i=0;i<=mx;++i) sm[i]=0;
    for(int i=1;i<=n;i++)
    {
        for(int e=fst[i];e;e=nxt[e])
        {
            int b=vb[e];
            if(!cap[e]||d[b]+1!=d[i]) continue;
            sm[d[i]]+=cap[e];
        }
    }
    int ans=2147483647;
    for(int i=1;i<=mx;++i)
        ans=min(ans,sm[i]);
    return ans;
}
}
ll F; int L,D,vs[SZ],ans=0;
bool bg(int e) {return cap[e]>=D+D+D;}
bool fb(int e) {return cap[e]<D+D+D;}
int val(int e)
{return !(bg(e)||(cap[e]>=D+D&&cap[e]<D+D+D
&&d[va[e]]==d[vb[e]]&&cap[e^1]>=D+D+D));}
GetD(fb)
int low[SZ],dfn[SZ],cur=0;
bool ins[SZ];
int ss[SZ],sn=0,bl[SZ],C=0;
void tarjan(int x)
{
    dfn[x]=low[x]=++cur; ss[++sn]=x; ins[x]=1;
    for esb(x,e,b)
    {
        if(vs[e]) continue;
        if(ins[b]) low[x]=min(low[x],dfn[b]);
        else if(!dfn[b])
            tarjan(b),low[x]=min(low[x],low[b]);
    }
    if(dfn[x]!=low[x]) return;
    int cs=0; ++C;
    while(1)
    {
        ++cs; int t=ss[sn--];
        bl[t]=C; ins[t]=0;
        if(t==x) break;
    }
}
int fa[SZ],tmp[SZ],oi[SZ]/*in-out*/;
void dfs(int x,int g)
{
    for esb(x,e,b)
    {
        if(bl[b]!=g||fa[b]) continue;
        fa[b]=x; dfs(b,g);
        cap[e]+=oi[b]; cap[e^1]-=oi[b];
        oi[x]+=oi[b]; oi[b]=0;
    }
}
//大概有错 
void doit()
{
    L=min(pow(n,2.0/3),pow(m,1.0/2));
    if(L<1) L=1; F=1e9;
    int tm=0;
    while(F>=1)
    {
        D=ceil(F/(ld)L);
        for(int i=1;i<=5*L;++i)
        {
            spfa();
            for(int i=2;i<=M;++i) vs[i]=val(i);
            sn=cur=C=0;
            for(int i=1;i<=n;i++)
                dfn[i]=low[i]=ins[i]=oi[i]=0;
            for(int i=1;i<=n;i++)
                if(!low[i]) tarjan(i);
            if(bl[S]==bl[T]) break;
            dn::M=1; dn::T=bl[T];
            for(int i=1;i<=C;i++) dn::fst[i]=0;
            for(int i=1;i<=n;i++) dn::d[bl[i]]=d[i];
            for(int i=2;i<=M;i++)
            {
                int a=va[i],b=vb[i];
                if(bl[a]!=bl[b]) tmp[i]=dn::addedge(bl[a],bl[b],cap[i],vs[i]);
            }
            ll cur=dn::dfs(bl[S],D); ans+=cur;
            if(!cur) break;
            for(int i=1;i<=n;i++)
                for esb(i,e,b)
                    if(bl[i]!=bl[b])
                    {
                        ll l=cap[e]-dn::cap[tmp[e]];
                        oi[i]-=l; cap[e]-=l;
                    }
            for(int i=1;i<=n;i++)
                if(!low[i]) fa[i]=-1,dfs(i,bl[i]);
        }
        F=min(int(F/2),og::gj());
    }
    printf("%d\n",ans);
}
}
namespace FF
{
char ch,B[1<<15],*S=B,*T=B;
#define getc() (S==T&&(T=(S=B)+fread(B,1,1<<15,stdin),S==T)?0:*S++)
#define isd(c) (c>='0'&&c<='9')
int aa,bb;int F(){
    while(ch=getc(),!isd(ch)&&ch!='-');ch=='-'?aa=bb=0:(aa=ch-'0',bb=1);
    while(ch=getc(),isd(ch))aa=aa*10+ch-'0';return bb?aa:-aa;
}
}
#define gi FF::F()
int main()
{
    n=gi,m=gi; S=1; T=n;
    for(int i=1;i<=m;++i)
    {
        int a,b,c;
        a=gi,b=gi,c=gi;
        ng::ad_de(a,b,c);
    }
    ng::doit();
}

理论复杂度为O(min(n^{2/3},m^{1/2})mlog(n)log(mU)),其中U表示边的最大容量。

测试数据中n=3000,m=200000,11组数据,时限2s,采用各种姿势生成,目测挺强的。如果按理论复杂度这些算法一个也过不了

sap:100分/492ms,dinic:100分/508ms,假的goldberg-rao:100分/1812ms,HLPP:45~54分(写挂了)

终于填上了这个大坑好高兴啊,过几天开新坑。

posted @ 2017-03-29 22:59 fjzzq2002 阅读(...) 评论(...) 编辑 收藏