网络流算法。Edmond-Karp，Dinic，ISAP的实现与性能分析

网络流很有用啊， 对于有些规则很复杂的图，很有可能就可以用网络流来求解。
首先，网络流是求在一个有向图中，边存在有流量限制也即容量（有些还存在费用问题）。给定一个起点，一个终点，求从起点到终点通过各种途径能传送的最大流量（在最小费用下的最大流量）。
初步，我们不考虑费用问题。 
先按照正常的人类思维来考虑这个问题：
<ol>
<li>先把图画出来，然后从起点往终点连线；</li>
<li>连到不能连后再看有没有哪条路不走会得到更好的结果，对之前的结果加以调整；</li>
<li>计算出每条路径的流量值</li>
<li>得到出最大流。</li>
</ol>
可以看出，有一个很重要的过程就是，如果我们选择了一个路径，我们还需要一个后悔的途径。 
于是，引入了反向边的概念。 
如果我们从a流向b 100单位的流量，同时记录下从b流向a -100单位的流量。如果这条边的流量不应该用在这个路径上时，真正应该使用这条边流量的路径可以通过反向边来调整流量结果。
如果我们找到了从起点到终点的一条路径，我们应该如何计算流量？ 
计算路径上所有边容量的最小值，然后这些边的流量同时加上这个最小值。 
这种每条路径慢慢往上增加流量的做法叫做增广路。 
可以得知，如果还存在一条路径，其容量减去流量的delta值不为0，那么总的流量还能再加上delta。
<hr>
增广路和反向边就是求解网络流问题的关键所在。
要实现反向边很容易，只需要在向图中插入<code>(from,to,cap,flow)</code>的边的同时插入<code>(to,from,0,0)</code>的边即可。 
每次在给一条边增加流量的同时，给其对应边减少相应的流量
而实现增广路，可以通过找到起点到终点的最短路，然后使用前置节点数组从终点返回起点计算delta，最后再次从终点到起点增加流量。 
而找最短路的方法则是使用BFS。 
不断重复BFS，直到起点和终点不连通（流量等于容量时，这条路就不能走了），获得的流量值的和即是最大流。
将这个思路实现出来，就是Edmond-Karp
<ol>
<li>使用BFS找到增广路</li>
<li>根据前置节点获取Delta</li>
<li>更新这条路径的流量（包括反向边流量）</li>
<li>重复上面的过程，直到起点和终点不连通。</li>
</ol>
<hr>
如果有一条起点到终点的路径，中间有很多分分合合的支路，那么使用Edmond-Karp将会把大量的时间用在BFS找最短路上。 
很可能上一次已经找到了一条增广路，下一条增广路其实只需要退一步换个方向即可，但是在Edmond-Karp里仍然要跑一次BFS。 
而这种“退一步，换个方向走”则很容易联想到DFS。
Dinic则是使用BFS分层（计算每个节点到起点的位置），然后使用DFS一层一层往下走，找增广路，找到后返回上一层继续往下找，直到该节点往后没有增广路，进行下一次BFS。 
这里有一个当前弧优化，在返回上一层继续往下找时，其实可以直接从上一条路径的下一条路径开始找就行，没有必要从第一条开始找。 
在图很复杂的时候，这里可以省下大量的时间。
将Dinic的算法实现就是：
<ol>
<li>根据BFS分层</li>
<li>根据分层结果进行DFS</li>
<li>保证按照层次往后进行DFS，来寻找增广路</li>
<li>找到增广路后，返回上一层，并且从当前路径的下一条路径继续找增广路。</li>
<li>当DFS找不到新的增广路后，再次进行新的BFS分层</li>
<li>当起点和终点不联通时结束程序</li>
</ol>
<h2 id="poj-1273postpoj1273">相关代码可以见<a href="/post/POJ/1273/">poj 1273</a></h2>
那么又可以发现，我已经DFS了，是不是DFS的同时也可以做一部分BFS的工作，甚至完全替代掉BFS呢? 
每次找到出度为0的节点的时候，就地更新分层来替代BFS的工作，这样就成为了ISAP算法。
如果某一层的节点数目为0，说明起点到终点不连通，可以直接结束，这是GAP优化
因此ISAP是一种对Dinic的进一步优化，两者时间复杂度是一样的，但是结合上各种剪枝优化，以及省去BFS的常数优化，其时间性能远好于Dinic。
思路如下：
<ol>
<li>先进行一次BFS预分层（到终点到的距离）</li>
<li>统计各层的节点个数（GAP优化用）</li>
<li>DFS遍历找增广路。</li>
<li>如果找到增广路，则按照之前的算法计算。</li>
<li>如果找不到能走的路径，则更新当前点的层数为周围能到达节点的最小值+1（这里不能走是因为不符合分层的走法，但是其实还是存在通路的。如果完全不存在通路，直接把距离设置成最大，这个点已经没有再考虑的必要了），然后返回上一层，继续往后找增广路（当前弧优化）</li>
<li>如果发现有一层的点数为0，则结束程序（GAP优化）</li>
<li>直到从起点到终点找不到任何路径，结束程序。</li>
</ol>
其中，ISAP的DFS改为循环写法， 可以进一步优化常数。
下面是一个ISAP模板。
<pre style="background-color:#fff">class maxFlow {
 struct Edge {
 int cap, flow;
 Edge(int _cap = 0, int _flow = 0) : cap(_cap), flow(_flow) {}
 };
 static const int INF = 0x7FFFFFFF;
 static const bool DEBUG = false;
 int s, v, n;

 Edge **edges;
 int *dis, *num, *pre, *cur;

 public:
 void addEdge(int from, int to, int cap) {
 edges[from][to].cap += cap;
 }

 int ISAP() {
 bfs();
 layerCalc();
 return dfs();
 }
 maxFlow(int s, int v, int n) {
 this-&gt;s = s;
 this-&gt;v = v;
 this-&gt;n = n;

 edges = new Edge *[n];
 for (int i = 0; i &lt; n; ++i) {
 edges[i] = new Edge[n];
 memset(edges[i], 0, sizeof(Edge) * n);
 }

 dis = new int[n];
 num = new int[n + 1];
 pre = new int[n];
 cur = new int[n];
 }
 ~maxFlow() {
 for (int i = 0; i &lt; n; ++i)
 delete[] edges[i];
 delete[] edges;

 delete[] dis;
 delete[] num;
 delete[] pre;
 delete[] cur;
 }

 private:
 queue&lt;int&gt; Q;

 void bfs() {
 while (!Q.empty())
 Q.pop();
 memset(dis, 0, sizeof(int) * n);
 Q.push(v);
 dis[v] = 1;
 while (!Q.empty()) {
 int t = Q.front();
 Q.pop();
 for (int i = 0; i &lt; n; ++i) {
 Edge &amp;e = edges[i][t];
 if (e.cap &gt; e.flow &amp;&amp; !dis[i]) {
 dis[i] = dis[t] + 1;
 Q.push(i);
 }
 }
 }
 }

 void layerCalc() {
 memset(num, 0, sizeof(int) * (n + 1));
 for (int i = 0; i &lt; n; ++i)
 ++num[dis[i]];
 }

 int Augumemt() {
 int t = v, delta = INF;
 while (t != s) {
 int &amp;lastNode = pre[t];
 Edge &amp;e = edges[lastNode][t];
 delta = min(delta, e.cap - e.flow);
 t = lastNode;
 }
 t = v;
 while (t != s) {
 int &amp;lastNode = pre[t];
 Edge &amp;e = edges[lastNode][t];
 Edge &amp;e2 = edges[t][lastNode];
 e.flow += delta;
 e2.flow -= delta;
 t = lastNode;
 }
 return delta;
 }

 int dfs() {
 memset(pre, 0, sizeof(int) * n);
 memset(cur, 0, sizeof(int) * n);

 int flow = 0;
 int t = s;

 while (dis[s] &lt;= n) {
 if (t == v) {
 flow += Augumemt();
 t = s;
 }

 int finish = true;
 for (int i = cur[t]; i &lt; n; ++i) {
 Edge &amp;e = edges[t][i];
 if (e.cap &gt; e.flow &amp;&amp; dis[t] == dis[i] + 1) {
 finish = false;
 pre[i] = t;
 cur[t] = i;
 t = i;
 break;
 }
 }

 if (finish) {
 int m = n;
 for (int i = 0; i &lt; n; i++) {
 Edge &amp;e = edges[t][i];
 if (e.cap &gt; e.flow)
 m = min(m, dis[i]);
 }
 if (--num[dis[t]] == 0)
 break;
 ++num[dis[t] = m + 1];
 cur[t] = 0;
 if (t != s)
 t = pre[t];
 }
 }
 return flow;
 }
};
</pre>

EK、Dinic和ISAP