转载于 http://www.cppblog.com/guyuecanhui/articles/88393.html
最近在看网络流,把几个常用的算法总结下,正确性的证明和一些理论的东西就不写了,参看算法导论和神牛们的论文,我只写算法的理解和实现模板。
Ford-Fulkerson方法
每次找增广路,把这条路上的所有点的流量加上这条路上的残余容量,再找新的增广路,直到找不到为止,它有很多种实现方法,下面给出算法导论上的伪代码
Ford_Fulkerson( G, s, t )
{
for each edge( u, v )∈E[G]
do f[u,v]= 0
f[v,u]= 0
while there exists a path p from s to t in the residual network Gf
do Cf(p)= min{ Cf(u,v) | (u,v) is in p }
for each edge(u,v) in p
do f[u,v]+= Cf(p)
f[v,u]= -f[u,v]
Edmonds-Karp算法
就是用广度优先搜索来实现Ford-Fulkerson方法中对增广路径的计算,时间复杂度为O(VE2) (代码参考NOCOW),算法导论上说Edmonds-Karp选择源点s到汇点t的最短路径作为增广路径,但下个代码只是通过广度优先顺序找到一条增广路径,怀疑此算法为Ford-Fulkerson,而不是Edmonds-Karp,先放在这,以后再来考虑这个问题
#define VMAX 201
int n, m; //分别表示图的顶点数和边数
//书上的例子,检验上述算法
int c[VMAX][VMAX]={{0,16,13},{0,0,10,12},{0,4,0,0,14},{0,0,9,0,0,20},{0,0,0,7,0,4}};
int Edmonds_Karp( int s, int t ){ //输入源点和汇点
int p, q, queue[VMAX], u, v, pre[VMAX], flow= 0, aug;
while(true){
memset(pre,-1,sizeof(pre)); //记录父节点
for( queue[p=q=0]=s; p<=q; p++ ){ //广度优先搜索
u= queue[p];
for( v=0; v<n&&pre[t]<0; v++ )
if( c[u][v]>0 && pre[v]<0 )
pre[v]=u, queue[++q]=v;
if( pre[t]>=0 ) break;
}
if( pre[t]<0 ) break; //不存在增广路
aug= 0x7fff; //记录最小残留容量
for( u=pre[v=t]; v!=s; v=u,u=pre[u] )
if(c[u][v]<aug) aug=c[u][v];
for( u=pre[v=t]; v!=s; v=u,u=pre[u] )
c[u][v]-=aug, c[v][u]+=aug;
flow+= aug;
}
return flow;
}int main()
{
n=6;
m=10;
int f=Edmonds_Karp(0,5);
printf("%d\n",f);
return 0;
}
PoJ1273
#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#define _DEBUG 1
#define INF 0x3c3c3c3c
#define VMAX 201
int n, m; //分别表示图的顶点数和边数,注意与题目中意思相反
//书上的例子,检验上述算法
int c[VMAX][VMAX];
int parent[VMAX];
int mincost[VMAX];
bool isVisited[VMAX];
bool singleSource(int s,int t){//单源最短路径,找到增广路径返回true,否则返回false
memset(parent,-1,sizeof(parent));
memset(mincost,INF,sizeof(mincost));
memset(isVisited,false,sizeof(isVisited));
mincost[s]=0;
int i;
while(1){
//找到最小点下标
int u=-1;
for(i=1;i<=n;i++){
if(isVisited[i])//已访问过
continue;
if(u==-1 && mincost[i]<INF){
u=i;
continue;
}
if(mincost[i] < mincost[u]){
u = i;
}
}
if(u==-1){//没有增广路经到t
return false;
}
isVisited[u]=true;
if(u == t)//结束
return true;
for(i=1;i<=n;++i){
if(!isVisited[i] && c[u][i]>0 && mincost[u]+c[u][i]<mincost[i]){
mincost[i] = mincost[u]+c[u][i];
parent[i] = u;
}
}
}
}
int Edmonds_Karp( int s, int t ){ //输入源点和汇点
int p, q, queue[VMAX], u, v,flow= 0, aug;
while(true){
if(!singleSource(s,t))
break;
aug= 0x7fff; //记录最小残留容量
for( u=parent[v=t]; v!=s; v=u,u=parent[u] )
if(c[u][v]<aug) aug=c[u][v];
for( u=parent[v=t]; v!=s; v=u,u=parent[u] )
c[u][v]-=aug, c[v][u]+=aug;
flow+= aug;
}
return flow;
}
int main()
{
int i,j;
#if _DEBUG == 1
freopen("POJ1273.in","r",stdin);
#endif
while(scanf("%d",&m)!=EOF){
memset(c,0,sizeof(c));
scanf("%d",&n);
for(i=0;i<m;i++){
int u,v,cap;
scanf("%d %d %d",&u,&v,&cap);
c[u][v] += cap;//存在多重边
}
int f=Edmonds_Karp(1,n);
printf("%d\n",f);
}
return 0;
}
注意1、此题存在多重边,如:
3 2
1 2 4
1 2 5
1 2 6
输出应该为15
2、边是有向的,如
2 2
1 2 150
2 1 150
输出应该为150
3、要检查残余网络的计算是否正确,如
10 8
1 3 2
3 4 2
4 5 2
5 8 2
4 6 4
1 2 4
2 6 2
6 7 2
7 8 4
2 7 2
输出应该为6
(这个对一些童鞋的算法不一定是存在残余网络,那就要自己再想一个咯)
The Push-Relabel Algorithm 最大流的压入与重标记算法
#include <algorithm>
#include <iostream>
#include <list>
#include <queue>
#define _DEBUG 1
using namespace std;
const int N = 201;
bool flag[N]; //顶点是否在队列 l 中
int c[N][N], //有向边的容量
e[N], //余流
f[N][N], //流
h[N], //高度
n; //顶点数
list<int> l, //待排除顶点
near[N]; //邻接表
void Push(int x, int y)
{
int df = min(e[x], c[x][y]-f[x][y]);
f[x][y] += df;
f[y][x] = -f[x][y];
e[x] -= df;
e[y] += df;
}
void Relabel(int x)
{
h[x] = n*2-1;
for (list<int>::iterator iter = near[x].begin(); iter != near[x].end(); ++iter)
if (h[*iter] < h[x] && c[x][*iter] > f[x][*iter])
h[x] = h[*iter];
++h[x];
}
void Discharge(int x)
{
list<int>::iterator iter = near[x].begin();
while (e[x] > 0) //有余流
{
if (iter == near[x].end())
{
Relabel(x);
iter = near[x].begin();
}
if (h[x] == h[*iter]+1 && c[x][*iter] > f[x][*iter])
{
Push(x, *iter);
if (e[*iter] > 0 && !flag[*iter])
l.push_back(*iter);
}
++iter;
}
}
int Push_Relabel()
{
l.clear();
//初始化前置流
h[0] = n; e[0] = 0;
memset(flag+1, 0, n-2); //1~n-2 不在 l 中
flag[0] = true; flag[n-1] = true; //防止源、汇进入 l
for (int i = 1; i < n; ++i)
{
h[i] = 0; //初始化各顶点高度为 h(i)=0
f[0][i] = c[0][i]; //初始化源与其他与之相连的 边流量f(0,i)=边容量c(0,i)
f[i][0] = -c[0][i]; //反向边容量
e[0] -= c[0][i]; //初始化源的 点余量e(0)=-c(0,i)
e[i] = c[0][i]; //初始化其他 点余量e(i)=c(0,i)
if (c[0][i] > 0 && i != n-1)
{
l.push_back(i); //待排除顶点
flag[i] = true;
}
}
//构造邻接表,每个点i的列表near[i]中存储与点i在图上相邻的点
for (int i = 0; i < n-1; ++i)
for (int j = i; ++j < n; )
if (c[j][i] > 0 || c[i][j] > 0)
{
near[i].push_back(j);
near[j].push_back(i);
};
//排除队列中的顶点
while (!l.empty())
{
int x = l.front();
Discharge(x);
l.pop_front();
flag[x] = false;
}
return e[n-1];
}
int main()
{
#if _DEBUG == 1
freopen("POJ1273.in","r",stdin);
#endif
int m;
while(scanf("%d",&m)!=EOF){
memset(c,0,sizeof(c));
memset(e,0,sizeof(e));
memset(f,0,sizeof(f));
memset(h,0,sizeof(h));
scanf("%d",&n);
for (; m; --m)
{
int x, y, w;
scanf("%d%d%d", &x, &y, &w);
c[x-1][y-1] += w;
}
printf("%d\n",Push_Relabel()); //计算从0(源)到 n-1(汇)的最大流
}
}
奇怪的是Push-Relabel Algorithm运行时间为47ms,而Edmonds_Karp算法只要16ms,从理论上讲Push-Relabel Algorithm应该快呀!
求最大流的Relabel-to-Front算法
转于:http://cuitianyi.com/blog/%E6%B1%82%E6%9C%80%E5%A4%A7%E6%B5%81%E7%9A%84relabel-to-front%E7%AE%97%E6%B3%95/
除了用各种方法在剩余网络中不断找增广路(augmenting)的Ford-Fulkerson系的算法外,还有一种求最大流的算法被称为压入与重标记(Push-Relabel)算法。它的基本操作有:压入,作用于一条边,将边的始点的预流尽可能多的压向终点;重标记,作用于一个点,将它的高度(也就是label)设为所有邻接点的高度的最小值加一。Push-Relabel系的算法普遍要比Ford-Fulkerson系的算法快,但是缺点是相对难以理解。
Relabel-to-Front使用一个链表保存溢出顶点,用Discharge操作不断使溢出顶点不再溢出。Discharge的操作过程是:若找不到可被压入的临边,则重标记,否则对临边压入,直至点不再溢出。算法的主过程是:首先将源点出发的所有边充满,然后将除源和汇外的所有顶点保存在一个链表里,从链表头开始进行Discharge,如果完成后顶点的高度有所增加,则将这个顶点置于链表的头部,对下一个顶点开始Discharge。
Relabel-to-Front算法的时间复杂度是O(V^3),还有一个叫Highest Label Preflow Push的算法复杂度据说是O(V^2*E^0.5)。我研究了一下HLPP,感觉它和Relabel-to-Front本质上没有区别,因为Relabel-to-Front每次前移的都是高度最高的顶点,所以也相当于每次选择最高的标号进行更新。还有一个感觉也会很好实现的算法是使用队列维护溢出顶点,每次对pop出来的顶点discharge,出现了新的溢出顶点时入队。
Push-Relabel类的算法有一个名为gap heuristic的优化,就是当存在一个整数0<k<V,没有任何顶点满足h[v]=k时,对所有h[v]>k的顶点v做更新,若它小于V+1就置为V+1。不过这个我还没实现,因为现在算法的效率我已经很满意了,以后有空再实现吧。
今天实现的Relabel-to-Front比昨天的距离标号最短增广路又要高效一倍左右。最后感叹一下,网络流真是有趣啊……甚至都想买那本专门讲网络流的书来看看。
示例程序(还是USACO Training Ditch,不过都是自己出数据测的):
/*
ID: dd.ener1
PROG: ditch
LANG: C++
*/
#include <cstdio>
#include <cstring>
using namespace std;
const int maxn=1024,OO=1000000000;
struct edge{
int x,y;//两个顶点
int w;//容量
int f;
edge *next,*back;
edge(){}
edge(int x,int y,int w,edge* next):x(x),y(y),w(w),f(0),next(next),back(0){}
void* operator new(size_t, void* p){
return p;
}
}*E[maxn],*cur[maxn];
int H[maxn];//高度函数
int P[maxn];//preflow
int N,S,T;
void input(){
int m;
scanf("%d%d",&m,&N);
memset(E,0,sizeof(E));
edge* data=new edge[2*m];
while(m--){
int x,y,w;
scanf("%d%d%d",&x,&y,&w);
E[x]=new((void*)data++) edge(x,y,w,E[x]);
E[y]=new((void*)data++) edge(y,x,0,E[y]);
E[x]->back = E[y];
E[y]->back = E[x];
}
S=1;
T=N;
}
template <class T>
const T& min(const T& a,const T& b){
return a<b?a:b;
}
void push(edge* e){
int u=e->x,v=e->y;
int delta=min(P[u],e->w);
e->w -= delta;
e->back->w += delta;
e->f += delta;
e->back->f = -(e->f);
P[u]-=delta;
P[v]+=delta;
}
void relabel(int u){
int minlabel=OO;
for(edge* e=E[u];e;e=e->next){
if(e->w==0)continue;
int v=e->y;
if(H[v]<minlabel){
minlabel=H[v];
cur[u]=e;
}
}
H[u]=minlabel+1;
}
void discharge(int u){
while(P[u]>0){
if(!cur[u]){
relabel(u);
}
int v=cur[u]->y;
if(cur[u]->w >0 && H[u]==H[v]+1)
push(cur[u]);
cur[u]=cur[u]->next;
}
}
void init_preflow(){
memset(H,0,sizeof(H));
memset(P,0,sizeof(P));
H[S]=N;
P[S]=OO;
for(edge* e=E[S];e;e=e->next)
push(e);
}
struct vertex{
int x;
vertex* next;
vertex(int x,vertex* next):x(x),next(next){}
vertex(){}
void* operator new(size_t, void* p){
return p;
}
};
vertex* build_list(){
vertex* data=new vertex[N-2];
vertex* L=0;
for(int i=N;i>=1;--i)
if(i!=S&&i!=T)
L=new((void*)data++) vertex(i,L);
return L;
}
int maxflow(){
memcpy(cur,E,sizeof(E));
init_preflow();
vertex* L=build_list();//list
vertex *v=L,*pre=0;
while(v){
int old_height=H[v->x];
discharge(v->x);
if(H[v->x]>old_height){
if(pre){
pre->next = v->next;
v->next = L;
L = v;
}
else;//本来就在表头
}
pre = v;
v = v->next;
}
return P[T];
}
int main(){
freopen("ditch.in","r",stdin);
freopen("ditch.out","w",stdout);
input();
int res=maxflow();
printf("%d\n",res);
}
Dinic算法
参考:
http://www.nocow.cn/index.php/Dinic
http://hi.baidu.com/_green_hand_/item/a03e9505e8dc0215acdc7076
SAP算法(by dd_engi):求最大流有一种经典的算法,就是每次找增广路时用BFS找,保证找到的增广路是弧数最少的,也就是所谓的Edmonds-Karp算法。可以证明的是在使用最短路增广时增广过程不超过V*E次,每次BFS的时间都是O(E),所以Edmonds-Karp的时间复杂度就是O(V*E^2)。
如果能让每次寻找增广路时的时间复杂度降下来,那么就能提高算法效率了,使用距离标号的最短增广路算法就是这样的。所谓距离标号,就是某个点到汇点的最少的弧的数量(另外一种距离标号是从源点到该点的最少的弧的数量,本质上没什么区别)。设点i的标号为D[i],那么如果将满足D[i]=D[j]+1的弧(i,j)叫做允许弧,且增广时只走允许弧,那么就可以达到“怎么走都是最短路”的效果。每个点的初始标号可以在一开始用一次从汇点沿所有反向边的BFS求出,实践中可以初始设全部点的距离标号为0,问题就是如何在增广过程中维护这个距离标号。
维护距离标号的方法是这样的:当找增广路过程中发现某点出发没有允许弧时,将这个点的距离标号设为由它出发的所有弧的终点的距离标号的最小值加一。这种维护距离标号的方法的正确性我就不证了。由于距离标号的存在,由于“怎么走都是最短路”,所以就可以采用DFS找增广路,用一个栈保存当前路径的弧即可。当某个点的距离标号被改变时,栈中指向它的那条弧肯定已经不是允许弧了,所以就让它出栈,并继续用栈顶的弧的端点增广。为了使每次找增广路的时间变成均摊O(V),还有一个重要的优化是对于每个点保存“当前弧”:初始时当前弧是邻接表的第一条弧;在邻接表中查找时从当前弧开始查找,找到了一条允许弧,就把这条弧设为当前弧;改变距离标号时,把当前弧重新设为邻接表的第一条弧,还有一种在常数上有所优化的写法是改变距离标号时把当前弧设为那条提供了最小标号的弧。当前弧的写法之所以正确就在于任何时候我们都能保证在邻接表中当前弧的前面肯定不存在允许弧。
还有一个常数优化是在每次找到路径并增广完毕之后不要将路径中所有的顶点退栈,而是只将瓶颈边以及之后的边退栈,这是借鉴了Dinic算法的思想。注意任何时候待增广的“当前点”都应该是栈顶的点的终点。这的确只是一个常数优化,由于当前边结构的存在,我们肯定可以在O(n)的时间内复原路径中瓶颈边之前的所有边。
// 去除流量、残量数组,直接在原图上操作
// 参考了第一份代码
/*
ID:t-x.h1
PROB:ditch
LANG:C++
*/
#include<cstdio>
#define min(a,b) ((a)<(b)?(a):(b))
const int MAX=2100000000,MAXn=200+9;
int n,answer,g[MAXn][MAXn],d[MAXn],gap[MAXn],st=1,ed=n;
int sap(int u,int flow)
{
if(u==ed)
return flow;
int res=flow,i,t;
for(i=1;i<=n;++i)
if(g[u][i] && d[u]==d[i]+1)
{
t=sap(i,min(res,g[u][i]));
g[u][i]-=t,g[i][u]+=t;
if(!(res-=t))
return flow;
}
if(!--gap[d[u]])
d[st]=n;
++gap[++d[u]];
return flow-res;
}
int main()
{
int m,i,t1,t2,t3;
freopen("ditch.in","r",stdin);freopen("ditch.out","w",stdout);
scanf("%d%d",&m,&n);
for(i=1;i<=m;++i)
{
scanf("%d%d%d",&t1,&t2,&t3);
g[t1][t2]+=t3;
}
for(gap[0]=n;d[st]<n;)
answer+=sap(st,MAX);
printf("%d\n",answer);
return 0;
}