poj 2391 Ombrophobic Bovines 二分 拆点 floyd 最大流 sap

缘起

【1】中用EK 被T掉的. 很不甘心, 后来了解到大家用的都是dinic或者sap这种主流姿势过的. 遂学习sap。 题目啥的见【1】好了.

分析

不得不说~ SAP真的是牛人想出来的非常高明的主意!!! 该算法敏锐的主意到了EK 为什么低效. 并且给出了极为优雅的解决方式。SAP是EK的优化, 但是本质依旧是EK

首先，之前【2】中展示了EK算法. 而EK的本质就是改进了FK，FK每次在残余网络上搜索那是相当盲目——毕竟dfs. 而EK比FK好点儿——起码EK有一个明确的标准——我每次在残余网络上就找最少的搜索步数（一条弧就是一步）搜到汇.

回到刚刚问的问题——为什么EK会慢?

有读者可能会说: 当然慢啦~ EK每次跑bfs才找到一条增广路径. 如果能一次性把所有增广路径找到的话, 哼哼~

这是一个问题，但是不是本质的. 即EK每次只找出一条残余网络的增广路径并不是导致EK效率低下的本质原因. 因为你一次找多条不就是原先的多次bfs合并成一次而已吗? 这对EK的加速其实没有本质增益作用的.

那EK为什么会慢? 因为EK每次的bfs是盲目的.

怎么讲? 看下图（这里引鉴了陈小玉老师的《趣学算法》的7.3.7节的图, 这本书讲这里讲的很有启发性, 推荐看一下, 秒懂~）

注意到虚线, 理解bfs的童鞋就立马会知道, 上图就是bfs搜索的过程（一条弧在bfs搜索中的视作步长为1的搜索）—— 就像是地震波那样一圈一圈的. 然后最终找到一条增广路 1->2->4->6 但是期间我们将3、5这两个和此条增广路径无关的顶点也搜过了. 这其实就是一种浪费. 而sap对ek的优化的本质——ek之所以会搜到3和5, 就是因为我们在bfs的时候拓展顶点的时候只要是可行弧（即容量>0而且邻接即可）就将其装入队列导致的.

而sap限定bfs搜索的时候只去搜索除了可行条件之外, 还限定只搜索满足

1
level[当前点] = level[拓展点]+1

的拓展点. 其中level[i]表示i到汇的最少搜索步长. 则这样就能通过bfs快速搜索到汇. 这样就提高了每次找到增广路径的速度. 而且sap每次都找到多条增广路径. 而不是调用一次（优化过的）bfs就找到一条. 当然, 这只是锦上添花 而已. sap优化ek的核心本质就在于level的引入.

那么还有一个问题, 如果按照上述原则找不到拓展点怎么办? 这样会导致搜索无法前进的. 可能明明有可行弧, 但就因为不满足上述level之间的关系而导致无法拓展. 很简单，没有的话, 我们就自己造. 当前节点的level不是不满足上述level的关系吗? 那我们就修改当前节点的level值, 修改方式是令当前节点的level为与当前节点邻接的点（即存在可行弧）的level的最小值+1(为什么是最小值? 因为求的是最短增广路径啊~ sap依旧还是bfs啊), 如果没有邻接点,则令当前节点的level值为inf（其实取节点数量就够了，因为任何一个顶点到汇的距离不会超出顶点个数-1, 你现在取值为顶点个数的话, 则该点将不会出现在任何增广路径中.）

通过这种方式重新赋值level[当前点]的话, 则当前点需要退回它上一个顶点（即它怎么来的），因为它上一个顶点的选择就不一样了. 所以要重新选择.

算法什么时候结束呢? 一种情况是源的level值>=节点数. 则源都不会加入增广路径了, 那肯定不会再有增广路径了. 另外一种情况是如果取某个level值的节点数量变成0了. 则意味着断层出现, 则再也不能出现增广路径了, 此时也要返回. 因此需要维护一个数组gap, 维护level取值为某个值的节点的个数.

综合上述认知, 不难搓一个sap的板子. 即所谓的gap优化.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
//#include "stdafx.h"
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <queue>
using namespace std;
//#define LOCAL
typedef long long LL;
const LL inf = ~0ull>>4, maxn = 205;

LL f,p,t, head[maxn<<1], rhead[maxn<<1],rcap[maxn*maxn<<1], cnt,rcnt, gg[maxn][maxn], tot1, tot2, ub, pre[maxn<<1], flow[maxn<<1], gap[maxn<<1], level[maxn<<1]; // gap[i]是level[j]=i的j的个数, level[i]表示顶点i距离汇的最短距离

struct Arc
{
	LL from,to,nxt,cap;
}g[maxn*maxn<<1]; 

void addArc(LL from, LL to, LL cap)
{
	g[cnt].from = from, g[cnt].to = to, g[cnt].nxt = head[from], g[cnt].cap = cap;
	head[from] = cnt++;
	g[cnt].from = to, g[cnt].to = from, g[cnt].nxt = head[to], g[cnt].cap = 0;
	head[to] = cnt++;
}

void floyd()
{
	for (LL k = 1; k<=f; k++)
	{
		for (LL i = 1; i<=f; i++)
		{
			for (LL j = 1; j<=f; j++)
			{
				gg[i][j] = min(gg[i][j],gg[i][k]+gg[k][j]);
			}
		}
	}
}

void build(LL x)
{
	for (LL i = 1; i<=f; i++)
	{
		for (LL j = i+1; j<=f; j++)
		{
			if (gg[i][j]<=x)
			{
				addArc(i,f+j,inf);
				addArc(j, i+f, inf);
			}
		}
	}
}

LL feasible(LL cur)
{
	for (int h = head[cur],to; ~h; h = g[h].nxt)
	{
		to = g[h].to;
		if (g[h].cap>0 && level[cur]==level[to]+1) // 如果是可行弧而且level函数满足条件
		{
			return h; // 返回下一个sap要搜索的顶点(注意, 返回的是弧的编号,因为本题用的是链式向前星构图)
		}
	}
	return -1;
}

void updata(LL s, LL t, LL d)
{
	LL ans = inf, h;
	while(t!=s)
	{
		h = pre[t];
		g[h].cap-=d;
		g[h^1].cap+=d;
		t = g[h].from;
	}
}

LL retreat(LL cur)
{
	LL ans = t+1;
	for (LL h = head[cur],to; ~h; h = g[h].nxt)
	{
		to = g[h].to;
		if (g[h].cap>0 && ans>level[to]+1) // 找最小的level[i]+1(因为sap依旧是求最短增广路啊~ 只是用gap优化了过程而已)
		{
			ans = level[to]+1;
		}
	}
	return ans;
}

LL sap(LL s, LL t)
{
	memset(gap, 0, sizeof(gap));
	memset(level, 0, sizeof(level));
	flow[s] = inf;
	LL ans = 0, cur = s, nxtArc, nxt; // cur是当前点, nxtArc是下一步转移的弧, nxt是cur借助nxtArc到达的顶点
	while(level[s]<=t) // 第一个退出点， h[s]>t的话, 表明不再可能有增广路, 因为s到t一共t+1个顶点, 最多t步长搜索就够了, 你现在达到了t+1
	{
		nxtArc = feasible(cur); // 获取通往下一个可行点的弧的编号
		if (~nxtArc) // 如果找到了
		{
			nxt = g[nxtArc].to;
			pre[nxt] = nxtArc;
			flow[nxt] = min(flow[cur], g[nxtArc].cap);
			cur = nxt; // 移动
			if (cur==t) // 找到了一条增广路
			{
				ans+=flow[t];
				updata(s,t,flow[t]); // 更新残余网络
				cur = s; // 重头再试图找出一条增广路
			}
		}
		else // 如果从cur出发没找到可行弧, 就要自己造了
		{
			if (!--gap[level[cur]]) // 断层出现, 自减的原因在于level[cur]将修改为不同的值(注意, 一定是不同的值,如果还是相同的值的话,则当时用feasible求nxt救不会得到-1了)
			{
				return ans;
			}
			gap[level[cur]=retreat(cur)]++; // 自己造
			if (cur!=s) // 重新造了, 如果当前节点不是源, 则回退一步. 因为pre[cur]的feasible函数需要重新计算
			{
				cur = g[pre[cur]].from; // 退回到上一个顶点
			}
		}
	}
	return ans;
}

bool ck(LL x)
{
	build(x); 
	return sap(0,t)==tot1;
}

void revover()
{
	cnt = rcnt;
	for (LL i = 0; i<cnt; i++)
	{
		g[i].cap = rcap[i];
	}
	memcpy(head, rhead, sizeof(head)); 
}

LL kk()
{
	LL lo = 0, hi = ub, mid, ans = -1;
	while(lo<=hi)
	{
		mid = lo+hi>>1;
		if (ck(mid))
		{
			ans = mid;
			hi = mid-1;
		}
		else
		{
			lo = mid+1;
		}
		revover();
	}
	return ans;
}

int main()
{
#ifdef LOCAL
	freopen("d:\\data.in", "r", stdin);
	//	freopen("d:\\my.out", "w", stdout);
#endif
	scanf("%lld%lld", &f, &p);
	memset(head, -1, sizeof(head));
	t = 2*f+1; 
	for(LL i = 1,a,b; i<=f; i++)
	{
		scanf("%lld%lld", &a, &b);
		tot1+=a, tot2+=b;
		if (a)
		{
			addArc(0, i, a);
		}
		if (b)
		{
			addArc(f+i, t, b);
		}
		addArc(i,f+i, inf); 
	}
	rcnt = cnt; 
	memcpy(rhead, head, sizeof(head)); 
	for (LL i = 0; i<cnt; i++)
	{
		rcap[i] = g[i].cap;
	} 
	if (tot1>tot2)
	{
		puts("-1"); 
		return 0;
	}
	for (LL i = 1;i<=f; i++)
	{
		for (LL j = i; j<=f; j++)
		{
			gg[i][j]=gg[j][i] = i==j?0:inf;
		}
	}
	for (LL i = 1,a,b,c; i<=p; i++)
	{
		scanf("%lld%lld%lld", &a, &b, &c); 
		gg[a][b] = gg[b][a] = min(gg[a][b], c); 
	}  
	floyd(); 
	for (LL i = 1; i<=f; i++)
	{
		for (LL j = 1; j<=f; j++)
		{
			if (gg[i][j]<inf)
			{
				ub = max(ub, gg[i][j]);
			}
		}
	} 
	printf("%lld\n", kk());
	return 0;
}

ac情况

Status Accepted
Time 516ms
Memory 2984kB
Length 3956
Lang C++
Submitted 2019-09-29 12:15:06
Shared
RemoteRunId 20909745

ps: 这里多说一句，一般国外把EK叫做sap, 而本文叫做isap（improved sap）. 而国内一般EK就叫EK, 国外的isap就是本文的sap. 其实国外把EK叫sap是有道理的, 因为sap不就是shortest augment path么? 而EK不就是bfs找最短增广路么? 这里的最短的含义是搜索次数.

参考

【1】https://yfsyfs.github.io/2019/09/28/poj-2391-Ombrophobic-Bovines-%E4%BA%8C%E5%88%86-%E6%8B%86%E7%82%B9-floyd-%E6%9C%80%E5%A4%A7%E6%B5%81/

【2】https://yfsyfs.github.io/2019/09/25/poj-1273-Drainage-Ditches-EK%E6%A8%A1%E6%9D%BF%E9%A2%98-%E6%9C%80%E5%A4%A7%E6%B5%81/

【3】https://yfsyfs.github.io/2019/09/20/%E6%B4%9B%E8%B0%B7-3376-%E3%80%90%E6%A8%A1%E6%9D%BF%E3%80%91%E7%BD%91%E7%BB%9C%E6%9C%80%E5%A4%A7%E6%B5%81-Ford-Fulkerson%E7%AE%97%E6%B3%95/