最小费用最大流

最小费用最大流（Minimum Cost Maximum Flow），简称费用流（MCMF），是最大流算法的一个分支。

现实生活中，如运输物资，管道输水等等不仅要考虑最大流量，有时还要考虑运输费用问题，费用流问题由此而生。在一个图中，每条边不仅有一个流量值，还有一个费用值，一般表示是在这条边上流一个单位的流量所需的费用，所以对于每一条边上的流量flow[i]，实际花费为flow[i]*cost[i]。另外有些特殊情况一般很难遇见，以后在进行讨论。

求费用流的算法一般可以分为两类，一是每次都找最大流来更新最小费用直到无法更新，二是每次在剩余图中找以费用为权值的最短路来更新流量使其不断接近最大流直到没有增广路。在此笔者推荐第二种方法，并且介绍一种基于SPFA的算法。

Q:为什么不能用Dijkstra来找最短路呢？A：反向边的费用=正向边的费用的相反数，所以会出现负权边，用Dijkstra这不老脸一黑-_-||

以每条边的费用作为路径长度，每次从S到T找一条最短路（当然要满足路上所有边的流量>0），最后这条路上的流量就是流量最小的一条边的流量，不断重复这个步骤，直到找不到增广路为止，这时候得到的已有流量就是最大流。

那么问题来了，怎么存储增广路呢？掌握SPFA的同学都知道，结点的最短路值是会不断更新并加入队列的，此时我们不能简单的存储路径，而是采用逆序的思想。pre[i]表示i号结点是由第pre[i]条边走过来，那么i在路径中的上一个点就是from[pre[i]]，这样我们在找到一条增广路以后，从T开始不断沿着from[pre[i]]向前找，直到走到S，采用这种覆盖式的存储而不是简单地添加，就有效地避免了存进错误的结点或者重复存储等情况。那为什么不能直接存储i号结点的前趋结点，而是存它的前趋边再由此找到前趋结点呢？我们在处理增广路的时候，需要更新边上的流量值。由于采用了数组模拟链表的形式，很轻松的可以知道每条边对应的那些结点，但是只知道结点很难知道它对应的是哪一条边。那为什么不直接使用邻接矩阵存储呢？很多网络流图都是稀疏图，用邻接矩阵的话在DFS时复杂度会几何级增长，所以非常不建议使用，而且数组模拟链表稍微多写几次也很容易上手。具体参看代码。

 

好消息是，基于这种思想的MCMF算法代码非常简洁，不需要加入SAP或者另外的什么东西，主过程由两部分组成，一是SPFA，二是处理增广路，其中SPFA也只需加入两行代码即可（已在程序中标示出）。代码如下：

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
using namespace std;

const int maxn=10005;
const int maxm=1000005;
int node[maxm],next[maxm],cost[maxm],flow[maxm],head[maxn];
int n,m,ans,pre[maxn],tot,minf,dist[maxn];
bool vis[maxn];

void add(int x,int y,int z,int c){
	node[++tot]=y;
	next[tot]=head[x];head[x]=tot;
	flow[tot]=z;cost[tot]=c;
}

void init(){
	tot=-1;
	int x,y,z,c;			//z是流量 c是费用 
	cin>>n>>m;			
	memset(head,-1,sizeof head);
	for (int i=1;i<=m;i++) {
		cin>>x>>y>>z>>c;
		add(x,y,z,c);
		add(x,y,0,-c);		//为什么是-c这个就不用再说了吧 
	}
	ans=0;	
}

bool spfa(){
	int q[maxn*10],l,r,now,nex;		//由于spfa中可能会重复加进结点 所以队列要开大一点（开queue的当我没说） 
	memset(vis,1,sizeof vis);vis[1]=0;
	memset(dist,10,sizeof dist);dist[1]=0;
	q[1]=l=r=1;
	while (l<=r){
		now=q[l];vis[now]=1;
		for (int i=head[now];i!=-1;i=next[i])
			if (flow[i]>0){				//首先要保证可以走 
				nex=node[i];
				if (dist[nex]>dist[now]+cost[i]){
					dist[nex]=dist[now]+cost[i];
					pre[nex]=i;		//记录路径 
					if (vis[nex]){
						vis[nex]=0;
					 	q[++r]=nex;
					 	
					}
				}
			}
		l++;
	}
	return dist[n]!=168430090;
	//可以通过判断dist[n]有没有被更新过来判断图中是否存在增广路 memset中10等于168430090 见line35
}

int main(){
	init();
	
	while (spfa()){
		minf=0x7fffff;
		for (int i=n;i!=1;i=node[pre[i]^1])	//这里用了个小技巧使得少开了一个数组 
			minf=min(minf,flow[pre[i]]);	//i的出发结点就是它的反向边指向的结点 即node[pre[i]^1] 
	  						//这样我们就可以省去一个用于存储出发结点的数组 
		for (int i=n;i!=1;i=node[pre[i]^1])
			flow[pre[i]]-=minf,
			flow[pre[i]^1]+=minf;
		ans+=minf;
	}
	
	cout<<ans<<endl;
	
	return 0;
}