网络流--最大流(自用,勿看)
注意:这是一篇个人学习笔记,如果有人因为某些原因点了进来并且要看一下,请一定谨慎地阅读,因为可能存在各种奇怪的错误,如果有人发现错误请指出谢谢!
EK
每次用BFS寻找最短(边数最少)增广路,增广,直到没有增广路
可以证明,最多只需要O(nm)次增广
证明:
1 #include<cstdio>
2 #include<queue>
3 #include<algorithm>
4 #include<cstring>
5 using namespace std;
6 int n,m,s,t;
7 int flow[10100];
8 int pre[10100];//记录找增广路时每个结点的前驱结点到该结点的边
9 void read(int& x)
10 {
11 int c=getchar();
12 x=0;
13 while(c<'0'||c>'9') c=getchar();
14 while(c>='0'&&c<='9')
15 {
16 x=(x<<3)+(x<<1)+(c^'0');
17 c=getchar();
18 }
19 }
20 struct Edge
21 {
22 int from,to,cap,next;//cap表示容量
23 //Edge(int u,int v,int c):from(u),to(v),cap(c){}
24 }edge[201000];
25 int num_edge,first[10100];
26 int make_edge(int u,int v,int w)
27 {
28 edge[num_edge].from=u;
29 edge[num_edge].to=v;
30 edge[num_edge].next=first[u];
31 first[u]=num_edge;
32 edge[num_edge++].cap=w;
33 edge[num_edge].from=v;
34 edge[num_edge].to=u;
35 edge[num_edge].next=first[v];
36 first[v]=num_edge;
37 edge[num_edge++].cap=0;//由于要使用^1,所以必须Num_edge从0开始
38 }
39 int EK(int s,int t)
40 {
41 //整个算法是不断bfs找出一条增广路,记找出的增广路最大流量为flow[t],
42 //然后把增广路上前向边容量加上flow[t],反向边流量减去flow[t],直到没有增广路
43 int index,k,sumflow=0;
44 for(;;)
45 {
46 memset(pre,-1,sizeof(pre));
47 queue<int> q;
48 q.push(s);
49 flow[s]=0x3f3f3f3f;
50 //pre[s]=-1;
51 while(!q.empty())
52 {
53 index=q.front();
54 q.pop();
55 if(index==t) break;//bfs到汇点,就是找到了增广路
56 k=first[index];
57 while(k!=0)
58 {
59 if(pre[edge[k].to]==-1&&edge[k].cap)//&&edge[k].to!=s//如果pre[边指向的点]不为0则表示边指向的点已经搜索到过
60 {
61 pre[edge[k].to]=k;
62 flow[edge[k].to]=min(edge[k].cap,flow[index]);//这条增广路的最大流量是前面路径最大流量与这一条边最大流量的最小值
63 q.push(edge[k].to);
64 }
65 k=edge[k].next;
66 }
67 }
68 //这个bfs是为了找出一条增广路
69 if(pre[t]==-1) break;//如果找不到增广路了就退出
70 index=t;
71 while(index!=s)
72 {
73 edge[pre[index]].cap-=flow[t];//flow[t]就是这整一条增广路的最大流量
74 edge[pre[index]^1].cap+=flow[t];
75 index=edge[pre[index]].from;
76 }
77 sumflow+=flow[t];
78 }
79 return sumflow;
80 }
81 int main()
82 {
83 int i,u,v,w;
84 read(n);read(m);read(s);read(t);
85 for(i=1;i<=m;i++)
86 {
87 read(u);read(v);read(w);
88 if(u==v) continue;//||u==t||v==s
89 make_edge(u,v,w);
90 }
91 printf("%d",EK(s,t));
92
93 return 0;
94 }
ISAP
转自看一下
ISAP简化的描述是:
程序开始时用一个反向 BFS 初始化所有顶点的距离标号,之后从源点开始,进行如下三种操作:
(1)当前顶点 i 为终点时增广
(2) 当前顶点有满足 dist[i] = dist[j] + 1 的出弧时前进
(3) 当前顶点无满足条件的出弧时重标号并回退一步。
算法解释
概括地说,ISAP 算法就是不停地找最短增广路,找到之后增广;如果遇到死路就 retreat,直到发现s, t不连通,算法结束。找最短路本质上就是无权最短路径问题,因此采用 BFS 的思想。具体来说,使用一个数组d,记录每个节点到汇点tt的最短距离。搜索的时候,只沿着满足d[u]=d[v]+1的边u→v(这样的边称为允许弧)走。显然,这样走出来的一定是最短路。
原图存在两种子图,一个是残量网络,一个是允许弧组成的图。残量网络保证可增广,允许弧保证最短路(时间界较优)。所以,在寻找增广路的过程中,一直是在残量网络中沿着允许弧寻找。因此,允许弧应该是属于残量网络的,而非原图的。换句话说,我们沿着允许弧,走的是残量网络(而非原图)中的最短路径。当我们找到沿着残量网络找到一条增广路,增广后,残量网络肯定会变化(至少少了一条边),因此决定允许弧的d数组要进行相应的更新(顺便提一句,Dinic 的做法就是每次增广都重新计算d数组)。然而,ISAP 「改进」的地方之一就是,其实没有必要马上更新d数组。这是因为,去掉一条边只可能令路径变得更长,而如果增广之前的残量网络存在另一条最短路,并且在增广后的残量网络中仍存在,那么这条路径毫无疑问是最短的。所以,ISAP 的做法是继续增广,直到遇到死路,才执行 retreat 操作。
说到这里,大家应该都猜到了,retreat 操作的主要任务就是更新d数组。那么怎么更新呢?非常简单:假设是从节点u找遍了邻接边也没找到允许弧的;再设一变量m,令m等于残量网络中u的所有邻接点的d数组的最小值,然后令d[u]等于m+1即可。这是因为,进入 retreat 环节说明残量网络中u和 t已经不能通过(已过时)的允许弧相连,那么u和t实际上在残量网络中的最短路的长是多少呢?(这正是d的定义!)显然是残量网络中u的所有邻接点和t的距离加1的最小情况。特殊情况是,残量网络中u根本没有邻接点。如果是这样,只需要把d[u]设为一个比较大的数即可,这会导致任何点到u的边被排除到残量网络以外。(严格来说只要大于等于∣V∣即可。由于最短路一定是无环的,因此任意路径长最大是∣V∣−1)。修改之后,只需要把正在研究的节点u沿着刚才走的路退一步,然后继续搜索即可。
讲到这里,ISAP 算法的框架内容就讲完了。对于代码本身,还有几个优化和实现的技巧需要说明。
1.算法执行之前需要用 BFS 初始化d数组,方法是从t到s逆向进行。
2.算法主体需要维护一个「当前节点」u,执行这个节点的前进、retreat 等操作。
3.记录路径的方法非常简单,声明一个数组p,令p[i]等于增广路上到达节点ii的边的序号(这样就可以找到从哪个顶点到的顶点ii)。需要路径的时候反向追踪一下就可以了。
4.判断残量网络中s,t不连通的条件,就是d[s]≥∣V∣ 。这是因为当s,t不连通时,最终残量网络中s将没有任何邻接点,对s的 retreat 将导致上面条件的成立。
5.GAP 优化。GAP 优化可以提前结束程序,很多时候提速非常明显(高达 100 倍以上)。GAP 优化是说,进入 retreat 环节后,u,t之间的连通性消失,但如果u是最后一个和t距离d[u](更新前)的点,说明此时s,t也不连通了。这是因为,虽然u,t已经不连通,但毕竟我们走的是最短路,其他点此时到t的距离一定大于d[u](更新前),因此其他点要到t,必然要经过一个和t距离为d[u](更新前)的点。GAP 优化的实现非常简单,用一个数组记录并在适当的时候判断、跳出循环就可以了。
6.另一个优化,就是用一个数组保存一个点已经尝试过了哪个邻接边。寻找增广的过程实际上类似于一个 BFS 过程,因此之前处理过的邻接边是不需要重新处理的(残量网络中的边只会越来越少)。具体实现方法直接看代码就可以,非常容易理解。需要注意的一点是,下次应该从上次处理到的邻接边继续处理,而非从上次处理到的邻接边的下一条开始。
最后说一下增广过程。增广过程非常简单,寻找增广路成功(当前节点处理到tt)后,沿着你记录的路径走一遍,记录一路上的最小残量,然后从ss到tt更新流量即可。
实例
1 int source; // 源点 2 int sink; // 汇点 3 int p[max_nodes]; // 可增广路上的上一条弧的编号 4 int num[max_nodes]; // 和 t 的最短距离等于 i 的节点数量 5 int cur[max_nodes]; // 当前弧下标 6 int d[max_nodes]; // 残量网络中节点 i 到汇点 t 的最短距离 7 bool visited[max_nodes]; 8 9 // 预处理, 反向 BFS 构造 d 数组 10 bool bfs() 11 { 12 memset(visited, 0, sizeof(visited)); 13 queue<int> Q; 14 Q.push(sink); 15 visited[sink] = 1; 16 d[sink] = 0; 17 while (!Q.empty()) { 18 int u = Q.front(); 19 Q.pop(); 20 for (iterator_t ix = G[u].begin(); ix != G[u].end(); ++ix) { 21 Edge &e = edges[(*ix)^1]; 22 if (!visited[e.from] && e.capacity > e.flow) { 23 visited[e.from] = true; 24 d[e.from] = d[u] + 1; 25 Q.push(e.from); 26 } 27 } 28 } 29 return visited[source]; 30 } 31 32 // 增广 33 int augment() 34 { 35 int u = sink, df = __inf; 36 // 从汇点到源点通过 p 追踪增广路径, df 为一路上最小的残量 37 while (u != source) { 38 Edge &e = edges[p[u]]; 39 df = min(df, e.capacity - e.flow); 40 u = edges[p[u]].from; 41 } 42 u = sink; 43 // 从汇点到源点更新流量 44 while (u != source) { 45 edges[p[u]].flow += df; 46 edges[p[u]^1].flow -= df; 47 u = edges[p[u]].from; 48 } 49 return df; 50 } 51 52 int max_flow() 53 { 54 int flow = 0; 55 bfs(); 56 memset(num, 0, sizeof(num)); 57 for (int i = 0; i < num_nodes; i++) num[d[i]]++; 58 int u = source; 59 memset(cur, 0, sizeof(cur)); 60 while (d[source] < num_nodes) { 61 if (u == sink) { 62 flow += augment(); 63 u = source; 64 } 65 bool advanced = false; 66 for (int i = cur[u]; i < G[u].size(); i++) { 67 Edge& e = edges[G[u][i]]; 68 if (e.capacity > e.flow && d[u] == d[e.to] + 1) { 69 advanced = true; 70 p[e.to] = G[u][i]; 71 cur[u] = i; 72 u = e.to; 73 break; 74 } 75 } 76 if (!advanced) { // retreat 77 int m = num_nodes - 1; 78 for (iterator_t ix = G[u].begin(); ix != G[u].end(); ++ix) 79 if (edges[*ix].capacity > edges[*ix].flow) 80 m = min(m, d[edges[*ix].to]); 81 if (--num[d[u]] == 0) break; // gap 优化 82 num[d[u] = m+1]++; 83 cur[u] = 0; 84 if (u != source) 85 u = edges[p[u]].from; 86 } 87 } 88 return flow; 89 }
转自看一下3
这几天由于种种原因经常接触到网络流的题目,这一类型的题给人的感觉,就是要非常使劲的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%左右的下降,如果乐意修正的话可以进行改良,当然不修正不影响算法正确性)。
1 typedef struct {int v,next,val;} edge; 2 const int MAXN=20010; 3 const int MAXM=500010; 4 5 edge e[MAXM]; 6 int p[MAXN],eid; 7 8 inline void init(){memset(p,-1,sizeof(p));eid=0;} 9 10 //有向 11 inline void insert1(int from,int to,int val) 12 { 13 e[eid].v=to; 14 e[eid].val=val; 15 e[eid].next=p[from]; 16 p[from]=eid++; 17 18 swap(from,to); 19 20 e[eid].v=to; 21 e[eid].val=0; 22 e[eid].next=p[from]; 23 p[from]=eid++; 24 } 25 26 //无向 27 inline void insert2(int from,int to,int val) 28 { 29 e[eid].v=to; 30 e[eid].val=val; 31 e[eid].next=p[from]; 32 p[from]=eid++; 33 34 swap(from,to); 35 36 e[eid].v=to; 37 e[eid].val=val; 38 e[eid].next=p[from]; 39 p[from]=eid++; 40 } 41 42 int n,m;//n为点数 m为边数 43 int h[MAXN]; 44 int gap[MAXN]; 45 46 int source,sink; 47 inline int dfs(int pos,int cost) 48 { 49 if (pos==sink) 50 { 51 return cost; 52 } 53 54 int j,minh=n-1,lv=cost,d; 55 56 for (j=p[pos];j!=-1;j=e[j].next) 57 { 58 int v=e[j].v,val=e[j].val; 59 if(val>0) 60 { 61 if (h[v]+1==h[pos]) 62 { 63 if (lv<e[j].val) d=lv; 64 else d=e[j].val; 65 66 d=dfs(v,d); 67 e[j].val-=d; 68 e[j^1].val+=d; 69 lv-=d; 70 if (h[source]>=n) return cost-lv; 71 if (lv==0) break; 72 } 73 74 if (h[v]<minh) minh=h[v]; 75 } 76 } 77 78 if (lv==cost) 79 { 80 --gap[h[pos]]; 81 if (gap[h[pos]]==0) h[source]=n; 82 h[pos]=minh+1; 83 ++gap[h[pos]]; 84 } 85 86 return cost-lv; 87 88 } 89 90 int sap(int st,int ed) 91 { 92 93 source=st; 94 sink=ed; 95 int ret=0; 96 memset(gap,0,sizeof(gap)); 97 memset(h,0,sizeof(h)); 98 99 gap[st]=n; 100 101 while (h[st]<n) 102 { 103 ret+=dfs(st,INT_MAX); 104 } 105 106 return ret; 107 }