网络流费用流模板 ISAP+SPFA+ZKW

落日映苍穹つ 2021-11-02 14:42 323阅读 0赞

2020年4月20日重新发布。7年前的文章，几年前CSDN改版的时候变成了私密……重新发一下吧。

关于费用流ZKW算法的讲解：[从入门到精通: 最小费用流的“zkw算法”][_zkw]

关于Dinic的模板以及各模板的其他写法可以参考：[最大流（Dinic & Isap）+最小费用最大流（SPFAFlow） - 志当存高远][Dinic _ Isap_SPFAFlow_ -]

关于deque容器的学习：[C++ STL学习之三：容器deque深入学习][C_ STL_deque]

2013-8-21更新 ，修改了ISAP的模板；在最后贴一份别人用Dinic写的Poj2112。参考自：http://acm.hust.edu.cn/vjudge/problem/viewSource.action?id=1286342

Dinic我用的不多，暂时不仔细研究了。

#include <cstdio>
    #include <cstring>
    #include <iostream>
    using namespace std;
    
    const int INF = 0x3fffffff ;   //权值上限
    const int MAXPT = 100005 ;     //顶点数上限
    const int MAXEG = 120005 ;    //边数上限
    const int MAXQUE = 100005 ;	   // 队列长度
    
    /*
        s = 1 ; // 源点
    	t 根据初始化函数的不同改变
    */
    template<typename Type>
    class MNF_SAP
    {
    private:
        int s,t;
        int dis[MAXPT];   //距离标号
        int pre[MAXPT];   //前置顶点
        Type flow[MAXPT];  //到当前点为止，所有弧的最小值
        int curedge[MAXPT];   //当前弧cur
        int cnt[MAXPT];   //k标号出现次数
        int queue[MAXQUE],front,rear;
        bool vis[MAXPT];
    
        void BFS ()
        {
            int i,u;
            memset(vis,false,sizeof(vis));
            front=rear=0;
            dis[t]=0;
            vis[t]=true;
            queue[++rear]=t;
            while (front!=rear)
            {
                u=queue[(++front)%MAXQUE];
                for (i=head[u];i!=0;i=edges[i].next)
                    if (vis[edges[i].v]==false && !edges[i].cap)
                    {
                        dis[edges[i].v]=dis[u]+1;
                        vis[edges[i].v]=true;
                        queue[(++rear)%MAXQUE]=edges[i].v;
                    }
            }
    		for (i=1;i<=n;i++)
    			cnt[dis[i]]++;
        }
    public:
    
        struct Node
        {
            int v,next;
    		Type cap;
    		Node(){}
    		Node (int _v,Type _cap,int _next)
    		{
    			v=_v;
    			cap=_cap;
    			next=_next;
    		}
        }edges[MAXEG];
        int n;      //总节点数
    	int e;
        int head[MAXPT];
    
    	void Init (int _n)    //算法初始化
    	{
    		s=1,t=_n;   //源点1汇点n
    		n=_n;
    		e=2;
    		memset (head,0,sizeof(head));
    	}
    
    	void Init (int _s,int _t) //算法初始化,后续需要对n赋值
    	{
    		s=1;   //源点1,汇点指定
    		t=_t;
    		e=2;
    		memset (head,0,sizeof(head));
    	}
    
    	void Add (int u,int v,Type cap)   //始,终,量
    	{
    		edges[e]=Node(v,cap,head[u]);
    		head[u]=e++;
    		edges[e]=Node(u,0,head[v]);
    		head[v]=e++;
    	}
    
        Type SAP ()
        {
            int u,v,i;
    		Type maxflow=0;   //总最大流
            u=s;
            flow[s]=INF;
            for (i=1;i<=n;i++)
    			curedge[i]=head[i];     //当前弧初始化
            BFS ();
            cnt[0]=n;
            while (dis[s]<n)
            {
                for (i=curedge[u];i!=0;i=edges[i].next)        //找允许弧
                    if (edges[i].cap>0 && dis[edges[i].v]+1==dis[u])          //
    					break;
                if (i!=0)      //存在允许弧
                {
                    curedge[u]=i;         //设置当前弧
                    v=edges[i].v;
                    if (edges[i].cap<flow[u])
    					flow[v]=edges[i].cap;
                    else
    					flow[v]=flow[u];  //标记当前顶点为止经过的最小弧
                    u=v;
                    pre[v]=i;  //前置顶点边号
                    if (u==t)
                    {
                        do
                        {
                            edges[pre[u]].cap-=flow[t];          //正向弧减a[t]
                            edges[pre[u]^1].cap+=flow[t];        //通过异或操作找反向弧
                            u=edges[pre[u]^1].v;
                        }
                        while (u!=s);
                        maxflow+=flow[t];
                        //memset(flow,0,sizeof(flow));
                        flow[s]=INF;
                    }
                }
                else   //不存在允许弧
                {
                    if (--cnt[dis[u]]==0)
    					break;     //间隙优化
                    dis[u]=n;
                    curedge[u]=head[u];
                    for (i=head[u];i!=0;i=edges[i].next)
                        if (edges[i].cap && dis[edges[i].v]+1<dis[u])
                            dis[u]=dis[edges[i].v]+1;       //修改距离标号为 最小的非允许弧加1
                    cnt[dis[u]]++;
                    if (u!=s)
    					u=edges[pre[u]^1].v;
                }
            }
            return maxflow;
        }
    };
    
    MNF_SAP<int> ob;

#include <cmath>
    #include <queue>
    #include <cstdio>
    #include <cstring>
    #include <algorithm>
    using namespace std;
    #define max(a,b) ((a)>(b)?(a):(b))
    #define min(a,b) ((a)<(b)?(a):(b))
    
    const int INF = 0x3fffffff ;   //权值上限
    const int MAXPT = 5000*2+100 ;     //顶点数上限
    const int MAXEG = 120005 ;    //边数上限
    
    //template<typename Type>
    class MCMF_SPFA 
    {
    private:
    	int s,t,e;
    	int head[MAXPT],pre[MAXPT];
    	int flow[MAXPT],cost[MAXPT];
    	bool vis[MAXPT];
    
    	struct Node
    	{
    		int u,v,cost,flow,next;
    		Node(){}
    		Node (int _u,int _v,int _cost,int _flow,int _next)
    		{
    			u=_u;
    			v=_v;
    			cost=_cost;
    			flow=_flow;
    			next=_next;
    		}
    	}edges[MAXEG];
    
    	int SPFA ()
    	{
    		memset(pre,-1,sizeof(pre));
    		int i,u,v,f,c;
    		queue<int> Q;
    		for (i=0;i<=t;i++)
    			cost[i]=INF,vis[i]=false,flow[i]=0;
    		Q.push(s); 
    		cost[s]=0;
    		flow[s]=INF;
    		while (!Q.empty())
    		{
    			u=Q.front();
    			Q.pop();
    			vis[u]=false;
    			for (i=head[u];i!=-1;i=edges[i].next)
    			{
    				v=edges[i].v;
    				f=edges[i].flow;
    				c=edges[i].cost;
    				if (f>0 && cost[v]>cost[u]+c)
    				{
    					cost[v]=cost[u]+c;
    					flow[v]=min(flow[u],f);
    					pre[v]=i;
    					if (vis[v]==false)
    						vis[v]=true,Q.push(v);
    				}
    			}
    		}
    		return flow[t];
    	}
    public:	
    	void Init (int _s,int _t)    //算法初始化
    	{
    		s=_s;t=_t;
    		e=0;
    		memset (head,-1,sizeof(head));
    	}
    	
    	void Add (int u,int v,int flow,int cost)   //始,终,量
    	{
    		edges[e]=Node(u,v,cost,flow,head[u]);
    		head[u]=e++;
    		edges[e]=Node(v,u,-cost,0,head[v]);
    		head[v]=e++;
    	}
    
    	int MCMF ()
    	{
    		int ans=0,temp;
    		while (temp=SPFA ())
    			for (int i=pre[t];i!=-1;i=pre[edges[i].u])
    			{
    				ans+=temp*edges[i].cost;
    				edges[i].flow-=temp;
    				edges[i^1].flow+=temp;
    			}
    		return ans;
    	}
    }ob;

#include <cmath>
    #include <queue>
    #include <cstdio>
    #include <cstring>
    using namespace std;
    #define max(a,b) ((a)>(b)?(a):(b))
    #define min(a,b) ((a)<(b)?(a):(b))
    
    const int INF = 0x3fffffff ;   //权值上限
    const int MAXPT = 5000*2+100 ;     //顶点数上限
    const int MAXEG = 120005 ;    //边数上限
    
    //template<typename Type>
    class MCMF_ZKW
    {
    private:
    	int head[MAXPT],d[MAXPT],e;
    	bool vis[MAXPT];
    	int ans,cost,src,des,n;
    
    	struct Node
    	{
    		int v,flow,cost,next,re;
    		Node(){}
    		Node (int _v,int _flow,int _cost,int _next,int _re)
    		{
    			v=_v,flow=_flow,cost=_cost;
    			next=_next,re=_re;
    		}
    	}edges[MAXEG];
    
    	int aug (int u,int f)
    	{
    		if (u == des)
    		{
    			ans+=cost*f;
    			return f;
    		}
    		vis[u]=true;
    		int tmp=f;
    		for (int i=head[u]; i!=-1; i=edges[i].next)
    			if (edges[i].flow && !edges[i].cost && vis[edges[i].v]==false)
    			{
    				int delta = aug(edges[i].v, tmp < edges[i].flow ? tmp : edges[i].flow);
    				edges[i].flow -= delta;
    				edges[edges[i].re].flow += delta;
    				tmp -= delta;
    				if (tmp==0)
    					return f;
    			}
    		return f-tmp;
    	}
    
    	bool modlabel ()
    	{
    		for (int i=0;i<=n;i++)
    			d[i]=INF;
    		d[des] = 0;
    		deque<int>Q;
    		Q.push_back(des);
    		while (!Q.empty())
    		{
    			int u=Q.front(),tmp;
    			Q.pop_front();
    			for (int i=head[u]; i!=-1; i=edges[i].next)
    				if (edges[edges[i].re].flow && (tmp = d[u]-edges[i].cost) < d[edges[i].v])
    					(d[edges[i].v] = tmp) <= d[Q.empty() ? src : Q.front()] ? Q.push_front(edges[i].v) : Q.push_back(edges[i].v);
    		}
    		for (int u=1;u<=n;u++)
    			for (int i=head[u]; i!=-1; i=edges[i].next)
    				edges[i].cost += d[edges[i].v] - d[u];
    			cost+=d[src];
    			return d[src] < INF;
    	}
    public:
    	
    	void Init (int _s,int _t)    //算法初始化
    	{
    		src=_s;des=_t;
    		e=0,n=des;
    		memset (head,-1,sizeof(head));
    		ans=cost=0;
    	}
    
    	void Add (int u,int v,int flow,int cost)
    	{
    		edges[e]=Node(v,flow,cost,head[u],e+1);
    		head[u] = e++;
    		edges[e]=Node(u,0,-cost,head[v],e-1);
    		head[v] = e++;
    	}
    
    	int ZKW ()
    	{
    		while(modlabel())
    			do
    			{
    				memset(vis,false,sizeof(vis));
    			}while(aug(src,INF));
    		return ans;
    	}
    }ob;

#include <cstdio>
    #include <cstdlib>
    #include <cmath>
    #include <algorithm>
    #include <cstring>
    using namespace std;
    
    const int MAXPT = 255;
    const int MAXED = 100005;
    const int INF = 0x3fffffff;
    
    struct node{
        int u, v, flow;
        int opp;
        int next;
    };
    struct Dinic{
        node arc[MAXED];
        int vn, en, head[MAXPT];     //vn点个数(包括源点汇点),en边个数
        int cur[MAXPT];              //当前弧
        int q[MAXPT];                //bfs建层次图时的队列
        int path[MAXED], top;        //存dfs当前最短路径的栈
        int dep[MAXPT];              //各节点层次
        void init(int n){
            vn = n;
            en = 0;
            memset(head,-1,sizeof(head));
        }
        void insert_flow(int u, int v, int flow){
            arc[en].u = u;
            arc[en].v = v;
            arc[en].flow = flow;
            arc[en].opp = en + 1;
            arc[en].next = head[u];
            head[u] = en ++;
    
            arc[en].u = v;
            arc[en].v = u;
            arc[en].flow = 0;       //反向弧
            arc[en].opp = en - 1;
            arc[en].next = head[v];
            head[v] = en ++;
        }
        bool bfs(int s, int t){
    		memset(dep,-1,sizeof(dep));
            int lq = 0, rq = 1;
            dep[s] = 0;
            q[lq] = s;
            while(lq < rq){
                int u = q[lq ++];
                if (u == t){
                    return true;
                }
                for (int i = head[u]; i != -1; i = arc[i].next){
                    int v = arc[i].v;
                    if (dep[v] == -1 && arc[i].flow > 0){
                        dep[v] = dep[u] + 1;
                        q[rq ++] = v;
                    }
                }
            }
            return false;
        }
        int solve(int s, int t){
            int maxflow = 0;
            while(bfs(s, t)){
                int i,j,k;
                for (i = 1; i <= vn; i ++)  cur[i] = head[i];
                for (i = s, top = 0;;){
                    if (i == t){
                        int mink;
                        int minflow = 0x3fffffff;
                        for (k = 0; k < top; k ++)
                            if (minflow > arc[path[k]].flow){
                                minflow = arc[path[k]].flow;
                                mink = k;
                            }
                        for (k = 0; k < top; k ++)
                            arc[path[k]].flow -= minflow, arc[arc[path[k]].opp].flow += minflow;
                        maxflow += minflow;
                        top = mink;		//arc[mink]这条边流量变为0, 则直接回溯到该边的起点即可(这条边将不再包含在增广路内).
                        i = arc[path[top]].u;
                    }
                    for (j = cur[i]; j != -1; cur[i] = j = arc[j].next){
                        int v = arc[j].v;
                        if (arc[j].flow && dep[v] == dep[i] + 1)
                            break;
                    }
                    if (j != -1){
                        path[top ++] = j;
                        i = arc[j].v;
                    }
                    else{
                        if (top == 0)   break;
                        dep[i] = -1;
                        i = arc[path[-- top]].u;
                    }
                }
            }
            return maxflow;
        }
    }dinic;
    
    int map[250][250];
    void floyd(int n){
        for (int k = 0; k < n; k ++){
            for (int i = 0; i < n; i ++){
                if(map[i][k] == INF) continue;
                for (int j = 0; j < n; j ++){
                    if (map[j][k] == INF)    continue;
                    if (map[i][j] > map[i][k] + map[k][j])
                        map[i][j] = map[i][k] + map[k][j];
                }
            }
        }
        return ;
    }
    int go(int mid, int k, int c, int m)
    {
    	int i;
        dinic.init(k+c+2);
        for (i = 1; i <= k; i ++){
            dinic.insert_flow(k+c+1, i, m);
        }
        for (i = k+1; i <= k+c; i ++){
            dinic.insert_flow(i, k+c+2, 1);
        }
        for (i = 0; i < k; i ++){
            for (int j = k; j < k+c; j ++){
                if (map[i][j] <= mid){
                    dinic.insert_flow(i+1, j+1, 1);
                }
            }
        }
        return dinic.solve(k+c+1, k+c+2);
    }
    int BS(int k, int c, int m){
        int l = 0, r = 10000000;
        while(l < r){
            int mid = (l+r)>>1;
            if (go(mid, k, c, m) == c){
                r = mid;
            }
            else{
                l = mid + 1;
            }
        }
        return r;
    }
    int main ()
    {
        int k, c, m;
        scanf("%d %d %d", &k, &c, &m);
        for (int i = 0; i < k+c; i ++){
            for (int j = 0; j < k+c; j ++){
                scanf("%d", &map[i][j]);
                if (map[i][j] == 0)
                    map[i][j] = INF;
            }
        }
        floyd(k+c);
        printf("%d\n", BS(k, c, m));
    	return 0;
    }

[_zkw]: http://www.artofproblemsolving.com/blog/54262
[Dinic _ Isap_SPFAFlow_ -]: http://blog.csdn.net/fp_hzq/article/details/6675872
[C_ STL_deque]: http://www.2cto.com/kf/201204/126761.html