ISAP算法学习笔记

Posted by Panda2134's Blog on March 7, 2018

ISAP,即改进的SAP算法。

算法基于:如果每次在残量网络中从汇点开始反向BFS, s-t 距离是单调不减的。

流程如下:

首先用BFS从汇点出发计算每个点到汇点的距离标号。

Advance:

每次走满足 $d_v = d_u - 1$ 的路径,即走最短路。维护每个点的当前弧。

(实践证明,只有ISAP需要当前弧优化……原始对偶加了之后反而更慢……Dinic加了之后有时候快很多,有时候慢一倍……)

Retreat:

当沿着允许弧无法增广的时候,说明需要更新距离标号,产生新的允许弧。

于是令 $d_u = \min \{ d_v \} + 1$ . 再从 s 开始走允许弧。特别地,当残量网络中 $u(u \neq t)$ 没有后继时,说明任何允许弧最终都不应该到达 $u$ 点,因为从 $u$ 点继续走无法到 $t$ ,更无法增广。于是此时令 $d_u = n$ ,即可排除 $u$ 点。( $n$ 个点的无权图中最大距离标号小于等于 $n-1$ )

GAP优化:设每个距离标号 $d$ 对应的点有 $num_d$ 个。如果某次 Retreat 操作后,$num_d = 0$ ,说明残量网络中 s-t 不连通,算法立刻结束。原因很简单:在这次修改距离标号之后,$d_s$ 只会增大。在从 $s$ 到 $t$ 的过程中,距离标号逐渐由 $d_s$ 一步一步地减小到0. 如果在这个过程中距离标号出现了断层(“gap”),s-t 最短路长度一定是无穷大,即残量网络中 s-t 不连通。于是当前流即为最大流。

代码

rk1,偷税

为啥第二名打表的都比我慢2333333

/*
 * ISAP.cpp - LOJ最大流模板
 *
 * Copyright (C) 2018 - panda_2134
 *
 * This is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with This. If not, see <http://www.gnu.org/licenses/>.
 */

#include <bits/stdc++.h>
using namespace std;

typedef long long int64;
const int MAXN = 100, MAXM = 5000;
const int64 LL_INF = 0x3f3f3f3f3f3f3f3fLL;

namespace Maxflow {
	struct Edge { 
		int u, v; 
		int64 flow, cap; 
		int next; 
	};
	int n, s, t, e_ptr = 1, head[MAXN+10]; Edge E[(MAXM+10)<<1];
	int p[MAXN+10], cur[MAXN+10], num[MAXN+10], d[MAXN+10], vis[MAXN+10];
	
	void AddEdge(int u, int v, int64 cap) {
		E[++e_ptr] = (Edge) { u, v, 0, cap, head[u] }; head[u] = e_ptr;
		E[++e_ptr] = (Edge) { v, u, 0,   0, head[v] }; head[v] = e_ptr;
	}
	
	int Hd, Tl, Q[MAXN+10];
	void BFS() {
	    Hd = 1, Tl = 0;
		memset(d, 0, sizeof(d));
		memset(vis, 0, sizeof(vis));
		Q[++Tl] = t; vis[t] = true; d[t] = 0;
		while(Hd <= Tl) {
			int u = Q[Hd++]; 
			for(int j=head[u]; j; j=E[j].next) {
				int v = E[j].v; int64 f = E[j^1].flow, c = E[j^1].cap;
				if(f < c && !vis[v]) {
					vis[v] = true; Q[++Tl] = v;
				}
			}
		}
	}
	
	int64 Augment() {
		int64 aug = LL_INF;
		for(int u = t; u != s; u = E[p[u]].u) 
			aug = min(aug, E[p[u]].cap - E[p[u]].flow);
		for(int u = t; u != s; u = E[p[u]].u) {
			E[p[u]].flow += aug;
			E[p[u]^1].flow -= aug;
		}
		return aug;
	}
	
	int64 ISAP() {
		int u = s; int64 flow = 0;
		BFS(); memcpy(cur, head, sizeof(head));
		for(int i = 1; i <= n; i++) ++num[d[i]];
		while(d[s] < n) {
			bool ok = false;
			if(u == t) {
				flow += Augment();
				u = s;
				continue;
			}
			for(int& j=cur[u]; j; j=E[j].next) { // Advance
				int v = E[j].v; int64 f = E[j].flow, c = E[j].cap;
				if(f < c && d[v] == d[u] - 1) {
					p[v] = j; u = v;
					ok = true;
					break;
				}
			}
			if(!ok) { // Retreat
				int slk = n-1;
				for(int j=head[u]; j; j=E[j].next) {
					int v = E[j].v; int64 f = E[j].flow, c = E[j].cap;
					if(f < c) slk = min(slk, d[v]);
				}
				if(--num[d[u]] == 0) break; // gap
				++num[d[u] = slk + 1];
				cur[u] = head[u]; // 重置当前弧
				if(u != s) u = E[p[u]].u;
			}
		}
		return flow;
	}
}

template<typename T>
inline void readint(T& x) {
	T f=1, r=0; char c=getchar();
	while(!isdigit(c)) { if(c=='-')f=-1; c=getchar(); }
	while(isdigit(c)) { r=r*10+c-'0'; c=getchar(); }
	x = f*r;
}

int main() {
	using namespace Maxflow;
	int m, u, v; int64 c;
	readint(n); readint(m); readint(s); readint(t);
	for(int i = 1; i <= m; i++) {
		readint(u); readint(v); readint(c);
		AddEdge(u, v, c);
	}
	cout << ISAP();
}