网络与流( Network and Flow \text{Network and Flow} Network and Flow)
请基本了解一下 网络流 的各种定义。
- 流网络(简称网络)是一张有向图 G = ( V , E ) G=(V,E) G=(V,E) 。其中有两个特殊点,源 s s s 和 汇 t t t 。
- ∀ ⟨ u , v ⟩ ∈ E \forall \langle u,v \rangle \in E ∀⟨u,v⟩∈E,存在一个 容量 c ( u , v ) ⩾ 0 c(u,v)\geqslant 0 c(u,v)⩾0;如果 ⟨ u , v ⟩ ∉ E \langle u,v \rangle \notin E ⟨u,v⟩∈/E ,定义 c ( u , v ) = 0 c(u,v)=0 c(u,v)=0 。
- 流量 是一个实值函数 f : V × V ↦ R f:V\times V\mapsto\reals f:V×V↦R,满足 f ( u , v ) ⩽ c ( u , v ) , f ( u , v ) = − f ( v , u ) , ∑ v ∈ V f ( u , v ) = 0 ( u ∉ { s , t } ) f(u,v)\leqslant c(u,v),\;f(u,v)=-f(v,u),\;\sum_{v\in V}f(u,v)=0\;(u\notin\{s,t\}) f(u,v)⩽c(u,v),f(u,v)=−f(v,u),∑v∈Vf(u,v)=0(u∈/{s,t}) 。
- 一个流 f f f 的 流值 ∣ f ∣ |f| ∣f∣ 定义为 ∑ u ∈ V f ( s , u ) \sum_{u\in V}f(s,u) ∑u∈Vf(s,u) 。
- 最大流 即是这个网络中流值最大的一个流。
- 残余流量 r ( u , v ) = c ( u , v ) − f ( u , v ) ⩾ 0 r(u,v)=c(u,v)-f(u,v)\geqslant 0 r(u,v)=c(u,v)−f(u,v)⩾0,注意这里 ⟨ u , v ⟩ ∈ E \langle u,v\rangle\in E ⟨u,v⟩∈E 不一定成立。
- 流 f f f 的 残留网络 G f = { ⟨ u , v ⟩ ∈ V × V ∣ r ( u , v ) > 0 } G_f=\{\langle u,v\rangle \in V\times V\;|\;r(u,v)>0\} Gf={⟨u,v⟩∈V×V∣r(u,v)>0} 。
- 流 f f f 的 增广路径 P f P_f Pf 是 G f G_f Gf 中的一条简单路径,起点为 s s s、终点为 t t t 。
那么网络流有什么性质呢?
- 可加性:将流 f f f 加上 G f G_f Gf 中的一个流,仍然是 G G G 的一个流。也就是说 f f f 可以加上一个增广路径 P f P_f Pf 。
- 最大流最小割定理:先证明流值总是不超过割容量。 ∣ f ∣ = ∑ i f ( s , i ) = ∑ i ∈ S ∑ j f ( i , j ) = ∑ i ∈ S ∑ j ∈ T f ( i , j ) ⩽ ∑ i ∈ S ∑ j ∈ T c ( i , j ) |f|=\sum_{i}f(s,i)=\sum_{i\in S}\sum_{j}f(i,j)=\sum_{i\in S}\sum_{j\in T}f(i,j)\leqslant\sum_{i\in S}\sum_{j\in T}c(i,j) ∣f∣=∑if(s,i)=∑i∈S∑jf(i,j)=∑i∈S∑j∈Tf(i,j)⩽∑i∈S∑j∈Tc(i,j) 。再说明可以取等。残余网络 G f G_f Gf 中 s s s 不可达 t t t,取 S S S 为 s s s 可达点的集合,则 cut ( S , T ) \operatorname{cut}(S,T) cut(S,T) 中正向边满流、反向边无流量(否则有正向残余流量),故 ∣ f ∣ = ∣ cut ( S , T ) ∣ |f|=|\operatorname{cut}(S,T)| ∣f∣=∣cut(S,T)∣,取等。
Edmonds–Karp \text{Edmonds–Karp} Edmonds–Karp 算法
这是一个非常普通的算法——找一个增广路径,然后把它增广进去。为了保证时间复杂度,我们找长度最短的一条增广路径。
长度最短,可以通过 b f s \tt{bfs} bfs 求得;增广进去时,应该增广一个大小为 δ ( P f ) = min ⟨ u , v ⟩ ∈ P f r ( u , v ) \delta (P_f)=\min_{\langle u,v\rangle \in P_f}r(u,v) δ(Pf)=min⟨u,v⟩∈Pfr(u,v) 的流。
int delta[MaxN]; /* s到i的增广路径上残余流量最小值 */
int pre[MaxN]; /* s到i的增广路径中的上一个点 */
int r[MaxN][MaxN]; /* 残余流量 */
int s, t, n; /* 源点、汇点、点的总数 */
int bfs(){
memset(delta,0,sizeof delta); /* 未访问 */
queue<int> q; q.push(s);
while(not q.empty()){
int x = q.front(); q.pop();
for(int i=1; i<=n; ++i)
if(delta[i] == 0 and r[x][i] > 0){
delta[i] = min(delta[x],r[x][i]);
pre[i] = x;
if(i == t) return delta[t];
/* 已经找到汇点,提前退出 */
q.push(i);
}
}
return 0; /* 无法增广 */
}
int EK(){
int maxFlow = 0;
while(true){
int d = bfs();
if(d == 0) return maxFlow;
for(int i=t; i!=s; i=pre[i]){
r[pre[i]][i] -= d;
r[i][pre[i]] += d;
/* “反对称性” */
}
maxFlow += d;
}
return maxFlow;
}
时间复杂度证明显然会跟路径长度挂钩(不然为何用 b f s \tt bfs bfs 呢)。当找到的增广路长度始终为某个值时,走过的边不会退流(因为这样的路径长度超过目前长度)。一次增广是 O ( V + E ) \mathcal O(V+E) O(V+E) 的,每次增广必然导致残余网络中一条边消失,故最多增广 E E E 次。所以 O ( V E + E 2 ) = O ( E 2 ) \mathcal O(VE+E^2)=\mathcal O(E^2) O(VE+E2)=O(E2) 时间内,就会让路径长度改变。
同样的道理,本次增广所带来的残留网络中边的改变,无非就是退流嘛。退流不会造出长度更短的路径,所以路径长度单增。最坏要增大
V
V
V 次,从
1
1
1 开始一步一步增大。所以复杂度
O
(
V
E
2
)
\mathcal O(VE^2)
O(VE2)
isap \text{isap} isap 算法
这是对 E K \tt{EK} EK 的改进。我们是这样想的——要是我们维持一个标号,让我们避开了许多坑点,并且维护起来并不难,不就很好吗?
然后, isap \text{isap} isap 算法诞生。我们引入 距离标号 d ( i ) d(i) d(i),满足这些性质:
- d ( t ) = 0 d(t)=0 d(t)=0
- ∀ ⟨ u , v ⟩ ∈ G f , d ( u ) ⩽ d ( v ) + 1 \forall \langle u,v \rangle \in G_f,\;d(u)\leqslant d(v)+1 ∀⟨u,v⟩∈Gf,d(u)⩽d(v)+1
显然 d ( i ) d(i) d(i) 不超过 i i i 到 t t t 的最短路长度(但是并非充要条件)。
所以,已经用不到的点(增广终了),可以将 d d d 值设为 n n n 。
这个 d d d 有什么用?答案:将边集重置为 E ′ = { ⟨ u , v ⟩ ∣ d ( u ) = d ( v ) + 1 } E'=\{\langle u,v \rangle \;|\; d(u)=d(v)+1\} E′={⟨u,v⟩∣d(u)=d(v)+1},在这个边集上求增广路径。
G A P \tt{GAP} GAP 优化:显然 d d d 值是连续的。所以,如果有 “空洞”(断层),那么一定是达到了最大流(即 G f G_f Gf 不存在 s s s 到 t t t 的路径)。
int c[MAXN][MAXN]; // 残留网络
int d[MAXN]; // d[]:距离标号
int vd[MAXN]; // vd[]:标号为i的结点个数
int n; /* 顶点数 */
int dfs(int i,int inFlow,const int &T){
// i:顶点; inFlow:最大有多大的流进入i; T:汇
int j, sum = 0, mind = n-1, delta;
if(i == T) // 到达汇点
return inFlow; /* 返回值为有多大的流进入T */
for(j = 1;j <= n; j++) // 枚举i的邻接点
if(c[i][j] > 0) { // 如果有边到j
if(d[i] == d[j]+1){// (i,j) in E'
delta = dfs(j,min(inFlow-sum,c[i][j]),T);
/* inFlow-sum:在i点剩下的流量; c[i,j]:这条边的容量 */
// 递归增广,返回沿(i,j)的实际增广量
c[i][j] -= delta; // 更新残留网络
c[j][i] += delta; /* 反对称性 */
sum += delta; // sum记录已经增广的流量
if(d[T] == -1)
// 结束,向上一层返回经过i的实际增广量
return sum;
if(sum == inFlow) break;
// 已经到达可增广上界,提前跳出
}
if (d[j] < mind) mind = d[j];
// 更新最小的邻接点标号
}
if(sum == 0) { // 如果从i点无法增广
vd[d[i]] --; // 标号为d[i]的结点数-1
if(vd[d[i]] == 0) // GAP优化
d[T] = -1; /* break标记 */
d[i] = mind + 1; // 更新标号
vd[d[i]] ++; // 新标号的结点数+1
}
return sum; // 向上一层返回经过i的实际增广量
}
int isap(){
int maxFlow = 0;
memset(d,0,sizeof d);
/* 显然,d全部为0是合法的 */
memset(vd,0,sizeof vd);
vd[0] = n; // all vertexes
while(d[T] != -1)
maxFlow += dfs(S,INF);
return maxFlow;
}
模板(不常用)
int c[MaxN][MaxN], d[MaxN], vd[MaxN], n;
int dfs(int x,int inFlow,const int &T){
int sum = 0, minD = n-1, delta;
if(x == T) return inFlow;
for(int i=1; i<=n; ++i)
if(c[x][i] > 0){
if(d[x] == d[i]+1){
delta = dfs(i,min(inFlow-sum,c[x][i]),T);
c[x][i] -= delta, c[i][x] += delta;
if((sum += delta) == inFlow) break;
if(d[T] == -1) return sum;
}
if (d[i] < minD) minD = d[i];
}
if(sum == 0){
if((-- vd[d[x]]) == 0) d[T] = -1;
++ vd[d[x] = minD+1];
}
return sum;
}
int isap(int s,int t){
int maxFlow = 0;
for(int i=1; i<=n; ++i) d[i] = vd[i] = 0;
for(vd[0]=n; d[t]!=-1; maxFlow+=dfs(s,infty,t));
return maxFlow;
}
时间复杂度
O
(
V
2
E
)
\mathcal O(V^2E)
O(V2E),然而我不会证明。
dinic \text{dinic} dinic 算法
是
isap
\text{isap}
isap 的兄弟版本。虽然复杂度一样。
换一种距离标号——让 d ( i ) d(i) d(i) 为 s s s 到 i i i 最短路的距离。这可以 b f s \tt{bfs} bfs 求出。讲完了。
int d[MaxN], q[MaxN];
bool bfs(int s,int t){
for(int i=1; i<=n; ++i) d[i] = -1;
d[s] = 0; int *head = q, *tail = q;
*(tail ++) = s;
while(head != tail){
int x = *(head ++);
for(int i=1; i<=n; ++i)
if(d[i] == -1 and c[x][i] > 0){
d[i] = d[x]+1;
*(tail ++) = i;
}
}
return d[t] != -1; /* 存在增广路 */
}
int dfs_T;
int dfs(int x,int inFlow){
if(x == dfs_T) return inFlow;
int sum = 0;
for(int i=1,delta; i<=n; ++i)
if(d[i] == d[x]+1 and c[x][i] > 0){
delta = dfs(i,min(inFlow-sum,c[x][i]));
c[x][i] -= delta, c[i][x] += delta;
if((sum += delta) == inFlow) break;
}
return sum;
}
int dinic(int s,int t){
int maxFlow = 0; dfs_T = t;
while(bfs(s,t)) maxFlow += dfs(s,infty);
return maxFlow;
}
重要优化
也许都不应该叫做 “优化”,直接叫做 dinic \text{dinic} dinic 的一部分,尤其是 “当前弧优化”。
- 吔屎优化:如果从一个点出发,其流量没有达到上界,则再给它更多的流量,也无法增广了。将其从残留网络中删去(具体来说,将 d d d 赋为 − 1 -1 −1 即可)。
if(sum != inFlow) d[x] = -1;
- 当前弧优化:记录一下,从当前点出发,已经访问到了哪条弧。因为访问过的边 ⟨ x , y ⟩ \lang x,y\rang ⟨x,y⟩ 要么是满流,要么 y y y 无法流出更多流量。
for(int &i=cur[x]; ~i; i=edge[i].nxt)
完整的代码如下(最初版):
int q[MaxN*MaxN<<2], d[MaxN*MaxN];
bool bfs(int s,int t){
int L = 0, R = 0;
for(int i=1; i<=n; ++i)
d[i] = -1;
d[s] = 0; q[++ R] = s;
while(L != R){
int x = q[++ L];
for(int i=head[x]; ~i; i=edge[i].nxt)
if(d[edge[i].to] == -1 and edge[i].val > 0){
d[edge[i].to] = d[x]+1;
q[++ R] = edge[i].to;
}
}
return d[t] != -1;
}
inline long long dfs(int x,long long inFlow,const int &T){
if(x == T) return inFlow;
long long sum = 0, delta;
for(int &i=cur[x]; ~i; i=edge[i].nxt) /* 当前弧 */
if(d[edge[i].to] == d[x]+1 and edge[i].val > 0){
delta = dfs(edge[i].to,min(inFlow-sum,edge[i].val),T);
edge[i].val -= delta, edge[i^1].val += delta;
if((sum += delta) == inFlow) break;
}
if(sum != inFlow) d[x] = -1; /* 吔屎优化 */
return sum;
}
long long dinic(int s,int t){
long long maxFlow = 0;
while(bfs(s,t)){
for(int i=1; i<=n; ++i)
cur[i] = head[i];
maxFlow += dfs(s,infty,t);
}
return maxFlow;
}
实践发展永无止境,解放思想永无止尽,码风是在不断变化的。所以我再贴一份,方便我复制。
struct Edge{
int to, nxt, capa;
Edge() = default;
Edge(int _to,int _nxt,int _capa)
:to(_to),nxt(_nxt),capa(_capa){ }
};
Edge e[MAXM<<1];
int head[MAXN], cntEdge;
void addEdge(int a,int b,int c,int d=0){
e[cntEdge] = Edge(b,head[a],c);
head[a] = cntEdge ++;
e[cntEdge] = Edge(a,head[b],d);
head[b] = cntEdge ++;
}
int dis[MAXN], que[MAXN];
void bfs(int x,const int &n){
memset(dis+1,-1,n<<2);
int *fro = que, *bac = que+1;
for(dis[*fro=x]=0; fro!=bac; ++fro)
for(int i=head[x=*fro]; ~i; i=e[i].nxt)
if(e[i].capa && dis[e[i].to] == -1){
dis[e[i].to] = dis[x]+1;
*(bac ++) = e[i].to;
}
}
int cur[MAXN];
int dfs(int x,int inFlow,const int &sink){
int sum = 0; if(x == sink) return inFlow;
for(int &i=cur[x]; ~i; i=e[i].nxt)
if(e[i].capa && dis[e[i].to] == dis[x]+1){
int d = dfs(e[i].to,min(inFlow-sum,e[i].capa),sink);
e[i].capa -= d, e[i^1].capa += d;
if((sum += d) == inFlow) break;
}
if(sum != inFlow) dis[x] = -1;
return sum;
}
const int INFTY = 0x7fffffff;
int dinic(int source,int sink,int n){
int res = 0;
while(bfs(source,n), dis[sink] != -1){
memcpy(cur+1,head+1,n<<2);
res += dfs(source,INFTY,sink);
}
return res;
}
时间复杂度
复杂度证明看这里:还是和
E
K
\tt EK
EK 类似,考虑增广路径长度不变时,每次增广会砍掉一条边——由于
d
i
n
i
c
\rm dinic
dinic 的标号,退流不会被考虑进来。那么最多砍
E
E
E 次,但是每次砍只需要
O
(
V
)
\mathcal O(V)
O(V) 的复杂度了!因为路径长度是
O
(
V
)
\mathcal O(V)
O(V) 的,
d
f
s
\tt dfs
dfs 回溯时就直接修改了!
O
(
E
)
\mathcal O(E)
O(E) 因为 当前弧优化 而不会乘上去。所以复杂度是
O
[
V
⋅
(
V
E
+
E
)
]
=
O
(
V
2
E
)
\mathcal O[V\cdot(VE+E)]=\mathcal O(V^2E)
O[V⋅(VE+E)]=O(V2E)
这也是为什么我说 当前弧
可能不是优化,而是固有部分。
跑二分图匹配的话,
d
i
n
i
c
\rm dinic
dinic 可以做到
O
(
m
n
)
\mathcal O(m\sqrt n)
O(mn) 的优秀复杂度。简要证明就是,考虑到边的容量为
1
1
1 ,一条增广路径会直接让一个点进入 吔屎优化
,所以一个分层图的完全增广变成了
O
(
V
+
E
)
\mathcal O(V+E)
O(V+E) 。而增广路径长度小于
n
\sqrt{n}
n 时,显然只有
O
(
n
)
\mathcal O(\sqrt{n})
O(n) 次;当路径长度不小于
n
\sqrt{n}
n 时,最大流与当前流的差是残余网络中很多个 点不交的路径,故最多
n
n
=
n
{n\over\sqrt{n}}=\sqrt{n}
nn=n 条路径(也就是剩余流量),所以最多增广这么多次。于是复杂度为
O
[
V
(
V
+
E
)
]
=
O
(
E
V
)
\mathcal O[\sqrt V(V+E)]=\mathcal O(E\sqrt V)
O[V(V+E)]=O(EV)