网络流ISAP算法的简单介绍（zz）

http://www.zlinkin.com/?p=34

 

　这几天由于种种原因经常接触到网络流的题目，这一类型的题给人的感觉，就是要非常使劲的YY才能出来点比较正常的模型。尤其是看了Amber最小割应用的文章，里面的题目思路真是充满了绵绵不绝的YD思想。然而比赛中，当你YD到了这一层后，您不得不花比较多的时间去纠结于大量细节的实现，而冗长的代码难免会使敲错版后的调试显得异常悲伤，因此一些巧妙简短高效的网络流算法在此时便显得犹为重要了。本文力求以最简短的描述，对比较流行的网络流算法作一定的总结，并借之向读者强烈推荐一种效率与编程复杂度相适应的算法。

　　众所周知，在网络流的世界里，存在2类截然不同的求解思想，就是比较著名的预流推进与增广路，两者都需要反向边的小技巧。

　　其中预流推进的算法思想是以边为单元进行推流操作。具体流程如下:置初始点邻接边满流并用一次反向bfs对每个结点计算反向距离标号，定义除汇点外存量大于出量的结点为活动结点，每次对活动结点按允许边（u->v:d[u]=d[v]+1）进行推流操作，直到无法推流或者该点存量为0，若u点此时仍为活动结点，则进行重标号，使之等于原图中进行推操作后的邻接结点的最小标号+1，并将u点入队。当队列为空时，算法结束，只有s点和t点存量非0，网络中各顶点无存量，无法找到增广路继续增广，则t点存量为最大流。

　　而增广路的思想在于每次从源点搜索出一条前往汇点的增广路，并改变路上的边权，直到无法再进行增广，此时汇点的增广量即为最大流。两者最后的理论基础依然是增广路定理，而在理论复杂度上预流推进要显得比较优秀。其中的HLPP高标预流推进的理论复杂度已经达到了另人发指的O（sqrt(m)*n*n），但是其编程复杂度也是同样的令人发指- -

　　于是我们能否在编程复杂度和算法复杂度上找到一个平衡呢，答案是肯定的。我们使用增广路的思想，而且必须进行优化。因为原始的增广路算法（例如EK）是非常悲剧的。于是有人注意到了预流推进中的标号法，在增广路算法中引入允许弧概念，每次反搜残留网络得到结点标号，在正向增广中利用递归进行连续增广，于是产生了基于分层图的Dinic算法。一些人更不满足于常规Dinic所带来的提升，进而加入了多路分流增广的概念，即对同一顶点的流量，分多路同时推进，再加上比较复杂的手工递归，使得Dinic已经满足大部分题目的需要。

　　然而这样做就是增广路算法优化的极限么？答案永远是不。人们在Dinic中只类比了预流推进的标号技术，而重标号操作却没有发挥得淋漓尽致。于是人们在Dinic的基础上重新引入了重标号的概念，使得算法无须在每次增广后再进行BFS每个顶点进行距离标号，这种主动标号技术使得修正后算法的速度有了不少提高。但这点提高是不足称道的，人们又发现当某个标号的值没有对应的顶点后，即增广路被截断了，于是算法便可以提前结束，这种启发式的优化称为Gap优化。最后人们结合了连续增广，分层图，多路增广，Gap优化，主动标号等穷凶极恶的优化，更甚者在此之上狂搞个手动递归，于是产生了增广路算法的高效算法–ISAP算法。

　　虽然ISAP算法的理论复杂度仍然不可超越高标预流推进，但其编程复杂度已经简化到发指，如此优化，加上不逊于Dinic的速率（在效率上手工Dinic有时甚至不如递归ISAP），我们没有不选择它的理由。

　　因此本文强烈推荐ISAP作为网络流首选算法。

　　其实现方法见下文，除去例行的添边操作，不超过50行的代码，何乐而不为之，以下实现仍有优化的余地（在计算初始标号时，为减小代码量直接忽略之，其复杂度不变，但实现后效率有5%左右的下降，如果乐意修正的话可以进行改良，当然不修正不影响算法正确性）。

 

typedef  struct {int v,next,val;} edge;
const int MAXN=20010;
const int MAXM=500010;
 
edge e[MAXM];
int p[MAXN],eid;
 
inline void init(){memset(p,-1,sizeof(p));eid=0;}
 
//有向
inline void insert1(int from,int to,int val)
{
    e[eid].v=to;
    e[eid].val=val;
    e[eid].next=p[from];
    p[from]=eid++;
 
    swap(from,to);
 
    e[eid].v=to;
    e[eid].val=0;
    e[eid].next=p[from];
    p[from]=eid++;
}
 
//无向
inline void insert2(int from,int to,int val)
{
    e[eid].v=to;
    e[eid].val=val;
    e[eid].next=p[from];
    p[from]=eid++;
 
    swap(from,to);
 
    e[eid].v=to;
    e[eid].val=val;
    e[eid].next=p[from];
    p[from]=eid++;
}
 
int n,m;//n为点数 m为边数
int h[MAXN];
int gap[MAXN];
 
int source,sink;
inline int dfs(int pos,int cost)
{
    if (pos==sink)
    {
        return cost;
    }
 
    int j,minh=n-1,lv=cost,d;
 
    for (j=p[pos];j!=-1;j=e[j].next)
    {
        int v=e[j].v,val=e[j].val;
        if(val>0)
        {
            if (h[v]+1==h[pos])
            {
                if (lv<e[j].val) d=lv;
                else d=e[j].val;
 
                d=dfs(v,d);
                e[j].val-=d;
                e[j^1].val+=d;
                lv-=d;
                if (h[source]>=n) return cost-lv;
                if (lv==0) break;
            }
 
            if (h[v]<minh)    minh=h[v];
        }
    }
 
    if (lv==cost)
    {
        --gap[h[pos]];
        if (gap[h[pos]]==0) h[source]=n;
        h[pos]=minh+1;
        ++gap[h[pos]];
    }
 
    return cost-lv;
 
}
 
int sap(int st,int ed)
{
 
    source=st;
    sink=ed;
    int ret=0;
    memset(gap,0,sizeof(gap));
    memset(h,0,sizeof(h));
 
    gap[st]=n;
 
    while (h[st]<n)
    {
        ret+=dfs(st,INT_MAX);
    }
 
    return ret;
}