最大流算法之三:ISAP <转>
(2009-08-14 19:24:27)
10
11 return f
procedure Advance(i)
1 设 (i,j) 为从 i 出发的一条允许弧
2 pi(j) <-- i
3 i <-- j
procedure Retreat(i)
1 d(i) <-- 1 + min{d(j):(i,j)属于残量网络Gf}
2 if i != s
3
procedure Augment
1 pi 中记录为当前找到的增广路 P
2 delta <-- min{rij:(i,j)属于P}
3 沿路径 P 增广 delta 的流量
4 更新残量网络 Gf
算法效率
O(E^2*V)
模板:
#include<cstdio>
#include<memory>
using namespace std;
const int maxnode = 1024;
const int infinity = 2100000000;
struct edge{
}*Net[maxnode];
int dist[maxnode]= {0}, numbs[maxnode] = {0}, src, des, n;
void rev_BFS(){
}
int maxflow(){
}
int main(){
}
ISAP 是图论求最大流的算法之一,它很好的平衡了运行时间和程序复杂度之间的关系,因此非常常用。
约定
我们使用邻接表来表示图,表示方法可以见文章带权最短路 Dijkstra, SPFA, Bellman-Ford, ASP, Floyd-Warshall 算法分析或二分图的最大匹配、完美匹配和匈牙利算法的开头(就不重复贴代码了)。在下文中,图的源点(source)表示为 s ,汇点(sink)表示为 t ,当前节点为 u 。建图时,需要建立双向边(设反向的边容量为0)才能保证算法正确。
引入
求解最大流问题的一个比较容易想到的方法就是,每次在残量网络(residual network)中任意寻找一条从 s 到 t 的路径,然后增广,直到不存在这样的路径为止。这就是一般增广路算法(labeling algorithm)。可以证明这种不加改进的贪婪算法是正确的。假设最大流是 f ,那么它的运行时间为 O( f⋅∣E∣) 。但是,这个运行时间并不好,因为它和最大流 f 有关。
人们发现,如果每次都沿着残量网络中的最短增广路增广,则运行时间可以减为 O(∣E∣2⋅∣V∣) 。这就是最短增广路算法。而 ISAP 算法则是最短增广路算法的一个改进。其实,ISAP 的意思正是「改进的最短增广路」 (Improved Shortest Augmenting Path)。
顺便说一句,上面讨论的所有算法根本上都属于增广路方法(Ford-Fulkerson method)。和它对应的就是大名鼎鼎的预流推进方法(Preflow-push method)。其中最高标号预流推进算法(Highest-label preflow-push algorithm)的复杂度可以达到 O(∣V∣2∣E∣−−−−√) 。虽然在复杂度上比增广路方法进步很多,但是预流推进算法复杂度的上界是比较紧的,因此有时差距并不会很大。
算法解释
概括地说,ISAP 算法就是不停地找最短增广路,找到之后增广;如果遇到死路就 retreat,直到发现 s, t 不连通,算法结束。找最短路本质上就是无权最短路径问题,因此采用 BFS 的思想。具体来说,使用一个数组 d ,记录每个节点到汇点 t 的最短距离。搜索的时候,只沿着满足 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 算法的框架内容就讲完了。对于代码本身,还有几个优化和实现的技巧需要说明。
- 算法执行之前需要用 BFS 初始化 d 数组,方法是从 t 到 s 逆向进行。
- 算法主体需要维护一个「当前节点」 u ,执行这个节点的前进、retreat 等操作。
- 记录路径的方法非常简单,声明一个数组 p ,令 p[i] 等于增广路上到达节点 i 的边的序号(这样就可以找到从哪个顶点到的顶点 i )。需要路径的时候反向追踪一下就可以了。
- 判断残量网络中 s, t 不连通的条件,就是 d[s]≥∣V∣ 。这是因为当 s, t 不连通时,最终残量网络中 s 将没有任何邻接点,对 s 的 retreat 将导致上面条件的成立。
- GAP 优化。GAP 优化可以提前结束程序,很多时候提速非常明显(高达 100 倍以上)。GAP 优化是说,进入 retreat 环节后, u, t 之间的连通性消失,但如果 u 是最后一个和 t 距离 d[u] (更新前)的点,说明此时 s, t 也不连通了。这是因为,虽然 u, t 已经不连通,但毕竟我们走的是最短路,其他点此时到 t 的距离一定大于 d[u] (更新前),因此其他点要到 t ,必然要经过一个和 t 距离为 d[u] (更新前)的点。GAP 优化的实现非常简单,用一个数组记录并在适当的时候判断、跳出循环就可以了。
- 另一个优化,就是用一个数组保存一个点已经尝试过了哪个邻接边。寻找增广的过程实际上类似于一个 BFS 过程,因此之前处理过的邻接边是不需要重新处理的(残量网络中的边只会越来越少)。具体实现方法直接看代码就可以,非常容易理解。需要注意的一点是,下次应该从上次处理到的邻接边继续处理,而非从上次处理到的邻接边的下一条开始。
最后说一下增广过程。增广过程非常简单,寻找增广路成功(当前节点处理到 t )后,沿着你记录的路径走一遍,记录一路上的最小残量,然后从 s 到 t 更新流量即可。
实现
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
|
int
source
;
// 源点
int
sink
;
// 汇点
int
p
[
max_nodes
]
;
// 可增广路上的上一条弧的编号
int
num
[
max_nodes
]
;
// 和 t 的最短距离等于 i 的节点数量
int
cur
[
max_nodes
]
;
// 当前弧下标
int
d
[
max_nodes
]
;
// 残量网络中节点 i 到汇点 t 的最短距离
bool
visited
[
max_nodes
]
;
// 预处理, 反向 BFS 构造 d 数组
bool
bfs
(
)
{
memset
(
visited
,
0
,
sizeof
(
visited
)
)
;
queue
<
int
>
Q
;
Q
.
push
(
sink
)
;
visited
[
sink
]
=
1
;
d
[
sink
]
=
0
;
while
(
!
Q
.
empty
(
)
)
{
int
u
=
Q
.
front
(
)
;
Q
.
pop
(
)
;
for
(
iterator_t
ix
=
G
[
u
]
.
begin
(
)
;
ix
!=
G
[
u
]
.
end
(
)
;
++
ix
)
{
Edge
&e
=
edges
[
(
*
ix
)
^
1
]
;
if
(
!
visited
[
e
.
from
]
&&
e
.
capacity
>
e
.
flow
)
{
visited
[
e
.
from
]
=
true
;
d
[
e
.
from
]
=
d
[
u
]
+
1
;
Q
.
push
(
e
.
from
)
;
}
}
}
return
visited
[
source
]
;
}
// 增广
int
augment
(
)
{
int
u
=
sink
,
df
=
__inf
;
// 从汇点到源点通过 p 追踪增广路径, df 为一路上最小的残量
while
(
u
!=
source
)
{
Edge
&e
=
edges
[
p
[
u
]
]
;
df
=
min
(
df
,
e
.
capacity
-
e
.
flow
)
;
u
=
edges
[
p
[
u
]
]
.
from
;
}
u
=
sink
;
// 从汇点到源点更新流量
while
(
u
!=
source
)
{
edges
[
p
[
u
]
]
.
flow
+
=
df
;
edges
[
p
[
u
]
^
1
]
.
flow
-
=
df
;
u
=
edges
[
p
[
u
]
]
.
from
;
}
return
df
;
}
int
max_flow
(
)
{
int
flow
=
0
;
bfs
(
)
;
memset
(
num
,
0
,
sizeof
(
num
)
)
;
for
(
int
i
=
0
;
i
<
num_nodes
;
i
++
)
num
[
d
[
i
]
]
++
;
int
u
=
source
;
memset
(
cur
,
0
,
sizeof
(
cur
)
)
;
while
(
d
[
source
]
<
num_nodes
)
{
if
(
u
==
sink
)
{
flow
+
=
augment
(
)
;
u
=
source
;
}
bool
advanced
=
false
;
for
(
int
i
=
cur
[
u
]
;
i
<
G
[
u
]
.
size
(
)
;
i
++
)
{
Edge
&
e
=
edges
[
G
[
u
]
[
i
]
]
;
if
(
e
.
capacity
>
e
.
flow
&&
d
[
u
]
==
d
[
e
.
to
]
+
1
)
{
advanced
=
true
;
p
[
e
.
to
]
=
G
[
u
]
[
i
]
;
cur
[
u
]
=
i
;
u
=
e
.
to
;
break
;
}
}
if
(
!
advanced
)
{
// retreat
int
m
=
num_nodes
-
1
;
for
(
iterator_t
ix
=
G
[
u
]
.
begin
(
)
;
ix
!=
G
[
u
]
.
end
(
)
;
++
ix
)
if
(
edges
[
*
ix
]
.
capacity
>
edges
[
*
ix
]
.
flow
)
m
=
min
(
m
,
d
[
edges
[
*
ix
]
.
to
]
)
;
if
(
--
num
[
d
[
u
]
]
==
0
)
break
;
// gap 优化
num
[
d
[
u
]
=
m
+
1
]
++
;
cur
[
u
]
=
0
;
if
(
u
!=
source
)
u
=
edges
[
p
[
u
]
]
.
from
;
}
}
return
flow
;
}
题目大意:n个工作,2个CPU,每个工作都可以工作在2个CPU上,并且对于工作i,在第一个CPU上工作的代价为ai,在第二个CPU上工作的代价是bi,再给m对工作,每对工作有数据交换,如果这2个工作工作在同一个CPU上,数据交换无代价,否则数据交换代价为w。求n个工作完成的最小代价。 题目分析:有2个CPU,每个工作只能工作在一个CPU上,所以我们可以将所有的工作分成2类:工作在第一个CPU上和工作在第二个CPU上的。于是就将所有的点分成了2个集合,每个集合中有若干个工作,并且2个CPU只能在2个集合中。很显然有割的性质,此题求一个最小割,根据最大流最小割定理,求最大流即可。 关于建图:要求最小代价。以2个CPU为源点和汇点。假设A为源点,B为汇点。对于工作i,源点与i建一条边,边权为bi,i与汇点建一条边,边权为ai。因为要求最小割,所以我们要想象割的过程,如果让工作i工作在A CPU上,那么就要割断工作i和B CPU的联系,那么割断这个联系的代价便是工作i工作在A CPU上的代价,所以i到汇点建边的边权要为ai。对于m个工作关系,同样想象割的过程:如果工作i,j工作在同一块CPU上,那么代价为0,如果工作在不同的CPU上,那么这2个工作必须分开,于是工作i和j直接就要有一个割,割的代价就是w,所以对于每对有联系的工作,之间相互建边,边权为w。 这题数据量比较大,测试模版比较不错。 详情请见代码:
|
题目大意:略。
题目分析:网络流模版题。不过数据很弱,只能测很烂的模版。
第一道网络流
详情请见代码:
- #include <iostream>
- #include<cstdio>
- #include<cstring>
- #include<algorithm>
- using namespace std;
- const int N = 205;
- const int M = 410;
- const int inf = 0x3f3f3f3f;
- struct node
- {
- int to,next,pre,c,f;
- }arc[M];
- int num;
- int head[N];
- int que[N];//bfs用
- int sta[N];//保存当前弧
- int rpath[N];//保存反向弧
- int cnt[N];
- int dis[N];
- int m,n;
- void build(int s,int e,int cap)//建图
- {
- arc[num].to = e;
- arc[num].c = cap;
- arc[num].f = 0;
- arc[num].next = head[s];
- head[s] = num ++;
- arc[num - 1].pre = num;//反向弧
- arc[num].pre = num - 1;
- arc[num].to = s;
- arc[num].c = 0;
- arc[num].f = 0;
- arc[num].next = head[e];
- head[e] = num ++;
- }
- void re_Bfs()
- {
- int i,front,rear;
- front = rear = 0;
- for(i = 1;i <= n;i ++)
- {
- dis[i] = inf;
- cnt[i] = 0;
- }
- que[rear ++] = n;
- cnt[0] = 1;
- dis[n] = 0;
- while(front != rear)
- {
- int u = que[front ++];
- for(i = head[u];i != -1;i = arc[i].next)
- {
- if(arc[arc[i].pre].c == 0 || dis[arc[i].to] < inf)
- continue;
- dis[arc[i].to] = dis[u] + 1;
- cnt[dis[arc[i].to]] ++;
- que[rear ++] = arc[i].to;
- }
- }
- }
- int ISAP()
- {
- int i,u,v,ret = 0;
- u = 1;
- for(i = 1;i <= n;i ++)
- sta[i] = head[i];
- while(dis[1] < n)
- {
- if(u == n)
- {
- int curflow = inf;
- for(i = 1;i != n;i = arc[sta[i]].to)
- curflow = min(curflow,arc[sta[i]].c);
- for(i = 1;i != n;i = arc[sta[i]].to)
- {
- arc[sta[i]].c -= curflow;
- arc[arc[sta[i]].pre].c += curflow;
- arc[sta[i]].f += curflow;
- arc[arc[sta[i]].pre].f -= curflow;
- }
- ret += curflow;
- u = 1;
- }
- for(i = sta[u];i != -1;i = arc[i].next)//寻找允许弧
- if(arc[i].c > 0 && dis[arc[i].to] + 1 == dis[u])
- break;
- if(i != -1)
- {
- sta[u] = i;
- rpath[arc[i].to] = arc[i].pre;
- u = arc[i].to;
- }
- else
- {
- if((--cnt[dis[u]]) == 0)//gap优化
- break;
- sta[u] = head[u];
- int Min = N;
- for(i = head[u];i != -1;i = arc[i].next)
- if(arc[i].c > 0)
- Min = min(Min,dis[arc[i].to]);
- dis[u] = Min + 1;
- cnt[dis[u]] ++;
- if(u != 1)
- u = arc[rpath[u]].to;
- }
- }
- return ret;
- }
- int main()
- {
- int i,u,v,c;
- while(scanf("%d%d",&m,&n) != EOF)
- {
- memset(head,-1,sizeof(head));
- num = 0;
- while(m --)
- {
- scanf("%d%d%d",&u,&v,&c);
- build(u,v,c);
- }
- re_Bfs();
- printf("%d\n",ISAP());
- }
- return 0;
- }