网络流——Dinic

上一篇写ISAP时说过要补上Dinic的博客。

如果为避免有些朋友不懂网络流，说Dinic之前还需要先介绍网络流以及增广路。这些我就不想介绍太多了（曾经想写网络流的博客也是因为不想写这些而放弃了），因此就只简单说。

网络流问题是一类图论问题，求解在网络流模型下的一些问题。网络流模型指一个有向图，包含源点和汇点，每条边有一个容量，且任意状态下每条边有一个流量。流量沿边的方向流动，总的方向为源点到汇点。直观地理解网络流模型的话，可以想象一个水管系统，源点就是无限流出水的点，汇点就是可以吸收无限水的点，边就是一些有流量上限的水管。有时候会把其他问题转化为网络流问题求解。

网络流里的状态有一些称为可行流，需要满足以下条件：1）每条边的流量不超过容量；2）除源点和汇点，对于任意一点，流入的流量和等于流出的流量和。

一般用S代表源点，T代表汇点。

增广路算法求解网络最大流（求一个可行流，使得流入T的流量和最大）。增广路算法求解最大流时首先从一个可行流构造残量网络：残量网络中，两点之间的有向边E(u,v)的边权f表示原图中从u到v还有f单位流量可以增加。举个例子：假设原图中有一条边E(u,v)，容量为c，此时流量为f，那么对应残量网络中就有两条边：E'(u,v,c-f)和E''(v,y,f)，反向边E''表示原图中有f单位流量可以减少（之后解释减少）。

同时算法构造一个层次网络：在残量网络上bfs，S的层为0，每个点的层是其到源点的最短路（不计边权）。

增广路算法构造好残量网络后，每进行一次bfs，就在层次网络上按照最短路进行增广。一条增广路是层次网络上一条从S到T的路径，且每条边都从一层指向下一层。只要满足层次网络上有增广路，就表明从S到T还可以增加流量。每次找到一条增广路后，使残量网络上路径中每条边的边权（表示可增加的流量）减去路径上最小的边权，同时反向边加上这个值。当一条边权值变为0，就从残量网络上暂时去掉这条边（bfs时不考虑这条边）。由于这条增广路不一定包含在最终方案里，所以反向边是用来反悔的，可以让这条增广路在需要的时候被消去。

由于增广一次后一定有一条边被删除了，每次增广到没有增广路后需要再次分层。

网络流的算法复杂度受到图的影响较大，因此分析出上界很松，若用邻接表，复杂度不超过O(|V|*|E|²)。看起来什么数据都跑不过，然而一般情况下达不到这么大。

接下来介绍Dinic。Dinic是增广路算法的优化版本。由于增广路算法每分一次层用bfs找增广路，每次只找到1条，效率比较低。Dinic把bfs换成dfs。dfs的原理是：从点u出发dfs，每条边寻找是否有增广路，全部加入到点u的增广路中。每次出发寻找都会分配一个最大可用流量，表示下一个点最多可能增广这些流量。S分配的最大可用流量是+∞（由于S需要流出无限水）；若点u出发可以到达点v，那么为点v分配的最大可用流量是点u的现存可用流量和E(u,v)的残量之间的最小值。出发寻找之前，现存可用流量是最大可用流量；每次找到一条增广路，现存可用流量减去这条增广路的最小残量。同时，E(u,v)的残量减去同样的值，反向边则加上这个值。

Dinic同样进行分层，每次dfs严格按照从一层到下一层的边增广。当dfs返回0时，重新分层。

Dinic的复杂度上界是O(|V|²*|E|)。看起来还是跑不过什么数据，然而大多数情况下比这个快到不知哪里去了。顺便一提，ISAP也是这个上界，不过加了gap优化后比Dinic还快。所以说网络流的复杂度都是玄学。

Dinic的小技巧：若只跑一次Dinic，可以不记录原图，只记录残量网络，加边时把原图中这条边的流量暂时定为0，这样可以保证原图是可行流。加边时从0号开始加边，正反两条边一起加，可以用位运算瞎搞，边i的反向边为i^1。

还有一点吐槽的：Dinic有一个当前弧优化，不过不让算法变慢的写法对我来讲依旧是个谜，因此没有写。

Dinic代码：

#include<cstdio>
#include<cstring>
#include<iostream>
#define min(a,b) (a<b?a:b)
#define MXN 10000+2
#define MXM 100000+10
int v[2*MXM],f[2*MXM],fst[MXN],nxt[2*MXM];
int layer[MXN];
int etop,S,T;
int n,m,x,y,z,w,ans;
int cur[MXM];
void add_edge(int x,int y,int z){
    v[etop] = y;
    f[etop] = z;
    nxt[etop] = fst[x];
    fst[x] = etop++;
    
    v[etop] = x;
    f[etop] = 0;
    nxt[etop] = fst[y];
    fst[y] = etop++;
}
bool DinicBfs(){
    int queue[MXN+100],head,tail;
    head=tail=0;
    memset(layer,-1,sizeof(layer));
    queue[tail++]=S;
    layer[S]=0;
    while(head<tail){
        int x=queue[head++];
        for(int y=fst[x]; y!=-1; y=nxt[y]){
            if(layer[v[y]]==-1&&f[y]>0){
                layer[v[y]]=layer[x]+1;
                queue[tail++]=v[y];
            }
        }
    }
    if(layer[T]==-1) return false;
    else return true;
}
int DinicDfs(int now,int curflow){
    if(now==T||curflow<=0) return curflow;
    int sum=0,temp;
    for(int &y=cur[now]; y!=-1; y=nxt[y]){
        if(layer[now]+1==layer[v[y]]&&f[y]>0){
            temp=DinicDfs(v[y],min(f[y],curflow));
            if(temp!=0){
                sum+=temp;
                f[y]-=temp;
                curflow-=temp;
                f[y^1]+=temp;
            }
        }
    }
    return sum;
}
int Dinic(){
    int sum=0,temp;
    while(DinicBfs()){
        memcpy(cur,fst,sizeof(cur));
        while(temp=DinicDfs(S,0x6f6f6f6f)) sum+=temp;
    }
    return sum;
}
int main(){
    scanf("%d%d%d%d",&n,&m,&S,&T);
    memset(fst,-1,sizeof(fst));
    memset(nxt,-1,sizeof(nxt));
    etop=0;
    for(int i=0; i<m; i++){
        scanf("%d%d%d",&x,&y,&z);
        add_edge(x,y,z);
    }
    ans=Dinic();
    printf("%d",ans);
    return 0;
}