算法学习笔记(16):网络流最大流

网络

是一种特殊的带权有向图，有以下几个特殊概念：

容量：边的权值，流在边上的最大量。
源点：没有入度的点
汇点：没有出度的点
流：每一条边都有对应的流量，流量不超过边的容量。除去汇点和源点，每个点的净流量等于0

本文考虑网络最大流问题。

网络流最大流

在最大流问题中，先假设只有一个源点和一个汇点。

对于网络G，给每条边指定对应的流量，使得源点和汇点的流量最大。

此时称这个分配 ”f( 边 ) = 流量“ 为最大流。

自然思路

建立一个图，每条边的值表示这个边的空闲流量，初始化和原图一样，因为一开始没有我们设定的流量。

进行以下操作的循环：

找到一个从源点到汇点的简单路径(无环)（任意算法）。
对于这条简单路径，找到这条路中的最小容量，在空闲流量中减去这个最小流量。
对于剩余容量归 0 的边将其删除或视作没有边。

直到找不到简单路径为止。

缺点

然而这个思路并不正确，并不能找到真正的最大流。

反例说明

考虑一个简单的网络：

源点 S 到中间节点 A1 有一条容量为 5 的边。
源点 S 到中间节点 A2 有一条容量为 2 的边。
A1到A2有一条容量为5的边。
A1到 T 有一条容量为5的边。
A1到 T 有一条容量为2的边。

方案1：

如果我们先找到 S - A2 - T ，两个都为2，所以都减去2.

再找到S- A1 - T 两个边都是5，所以都减去5.

总流量为7 = 2 + 5

方案2:

先找到了 S - A1 - A2 - T, 5 5 2,最小为2 ，都减去2，A2-T的边没有了

在找到 S - A1 - T 3 3 5，最小为3，都减去3

没有路径了

总流量为 5 = 2 + 3

显然方案2并不是最大流答案，这也被称为堵塞流，不能新增新的流量，当然，最大流也是堵塞流。

Ford–Fulkerson 增广

Ford–Fulkerson 增广是计算最大流的一类算法的总称。该方法运用贪心的思想，通过寻找增广路来更新并求解最大流。

这个算法可以说基于上面给的例子的自然思路。让我们找到自然思路中的问题，也就是当一条路把另一条路的容量抢走，导致不能拿到最多的量，一旦我们能撤销原来的量，就可以找到真正的最大量。

如何把这个撤销操作具现化呢？

Ford–Fulkerson 算法采用反向边的思路，例如：

对于一条边 A - B，假设我们已经找到了一条S1 - … - A1 - B1 - … - T1，把A - B之间的流量都占用了。此时，如果我们找到一条 S2 - … - B2 - A2 - … - T2的路径。(数字序号便于区分流量)

就可以有S2 - … - (B2 = B1) - … T1 , 有 S1 - … - (A1 = A2) - … - T2.

将原本不能走的路径，通过分配变得可行。

因此Ford–Fulkerson 算法步骤为：

找到一个从源点到汇点的简单路径(无环)（任意算法）。
对于这条简单路径，找到这条路中的最小容量，在空闲流量中减去这个最小流量。
对于剩余容量归 0 的边将其删除或视作没有边。
添加原路径的反向边，其权值为减去的流量值。边如果存在，则合并权值（新增！）

直到找不到边为止。

这个算法成功解决了最大流问题，得到了真正的最大流答案

缺点

遇到最坏情况会花费很多的时间，例如例子：

源点 S 到中间节点 A1 有一条容量为 100 的边。
源点 S 到中间节点 A2 有一条容量为 100 的边。
A1 到 A2有一条容量为 1 的边。
A1 到 T 有一条容量为 100 的边。
A1 到 T 有一条容量为 100 的边。

如果每次找路径都通过了这个A1 - A2的边，将每次只对最大网络流增加1，就会导致这个简单的网络流需要循环整整100次，如果最大网络流为F，最坏情况的复杂度有O(FV)

Edmonds–Karp 算法

Edmonds–Karp 算法属于Ford–Fulkerson 算法的一种特例，两者的思路一样。

不同点在于，在寻找最短路的第一步，EK算法将图看做无权图，寻找无权图中的最短路径。算法的好处在于不再依赖于最大流的大小，M为边的数量，N为点的数量，复杂度为O(M2N)。

每次找最短路使用BFS需要O(M)，被证明最多只要M*N轮（证明略）

代码：

struct Edge
{
    int from, to, cap, flow;

    Edge(int u, int v, int c, int f) : from(u), to(v), cap(c), flow(f) {}
};

struct EK
{   
    static const int maxn = 250;

    int n, m;             // n：点数，m：边数
    vector<Edge> edges;   // edges：所有边的集合
    vector<int> G[maxn];  // G：点 x -> x 的所有边的编号
    int a[maxn];          // a：点 x -> BFS 过程中最近接近点 x 的边给它的最大流
    int pathFrom[maxn];   // p：点 x -> BFS 过程中最近接近点 x 的边
                          

    void init(int n){
        for (int i = 0; i <= n; i++)
            G[i].clear();
        edges.clear();
    }

    void AddEdge(int from, int to, int cap){
        edges.push_back(Edge(from, to, cap, 0));
        edges.push_back(Edge(to, from, 0, 0));  // 提前建立反向边 权值是0
        m = edges.size();
        G[from].push_back(m - 2); // 插入了边的编号
        G[to].push_back(m - 1); 
    }

    int Maxflow(int s, int t){
        int flow = 0;
        while(true){
            memset(a, 0, sizeof(a));
            queue<int> BfsQ;
            BfsQ.push(s);
            a[s] = INF;
            while (not BfsQ.empty()){
                int x = BfsQ.front();
                BfsQ.pop();
                for (int edgeNo : G[x]){ // 遍历以 x 作为起点的边
                    Edge &e = edges[edgeNo];
                    if (!a[e.to] && e.cap > e.flow){
                        pathFrom[e.to] = edgeNo; // edgeNo 是最近接近点 e.to 的边
                        a[e.to] = min(a[x], e.cap - e.flow); // 最近接近点 e.to 的边赋给它的流
                        BfsQ.push(e.to);
                    }
                }
                if (a[t])
                    break; // 如果汇点接受到了流，就退出 BFS
            }
            if (!a[t])
                break; // 如果汇点没有接受到流，说明源点和汇点不在同一个连通分量上
            for (int u = t; u != s; u = edges[pathFrom[u]].from){  // 通过 u 追寻 BFS 过程中 s -> t 的路径               
                edges[pathFrom[u]].flow += a[t];     // 增加路径上边的 flow 值
                edges[pathFrom[u] ^ 1].flow -= a[t]; // 减小反向路径的 flow 值
            }
            flow += a[t];
        }
        return flow;
    }
};

增广成为了算法的重要概念，找增广路也就是不断地从反悔边和可走边中找到一条简单路径，而只要不断去寻找增广路，就可以得到正确的答案。

Dinic 算法

dinic算法修改了获得增广路的方式以优化复杂度，使得复杂度更多基于点。

引入了新概念分层图和堵塞流。

堵塞流前面已经提过，也就是一个流使得网络不能再新增流量，最大流一定是堵塞流，但堵塞流不一定是最大流。

分层图将图分成多个层，每个层表示为从起点出发，走几步可以到达的点。

类似于BFS分层，分层图不保留一层之间的点的边，只保留不同层之间的边，并且还是从小层到大层的边。

算法步骤：

建立剩余流量图的分层图
找到分层图的一个堵塞流
在剩余流量图中进行对应的减少以及建立反向边

直到分层图找不到堵塞流为止。

为什么使用分层图

我们已知算法的本质是基于不断寻找增广路而得到，根据EK算法，为了提高效率，我们通常选择最短路作为我们的增广路，如果使用分层图，可以使得我们获得的堵塞流中建立的增广路都是最短路径，通过一次预处理的分层而大大降低DFS寻找堵塞流的复杂度。

当前弧优化

在找到分层图的堵塞流中，我们需要提高效率以保证复杂度。

我们通常需要反复寻找增广路径。当前弧优化通过记录每个节点当前正在尝试的出边（称为“当前弧”），以避免重复检查已经确定不能作为增广路径的一些边，从而加速算法的执行。

对于原本的dfs，可能会重复尝试已经无效了的边，使用当前弧优化，如下面代码的current数组，始终指向有效的边，避免dfs再从已经无效的边进行检查。

这个被维护的current数组也就是下一个有效边的指针，我们又称作有效弧指针，故作有效弧优化。

代码：

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

const int MAX_NODES = 250;
const int INF = 0x3f3f3f3f;
const int MAX_EDGES = 10000;

struct MaxFlow {
    // 边结构体
    struct Edge {
        int to, next, capacity, flow;
    } edges[MAX_EDGES];

    int head[MAX_NODES], edge_count = 0;
    int node_count, source, sink;
    long long max_flow = 0;
    int depth[MAX_NODES], current[MAX_NODES];

    // 初始化
    void init() {
        memset(head, -1, sizeof head);
        edge_count = 0;
    }

    // 添加边 u->v，容量为capacity
    void addEdge(int u, int v, int capacity) {
        edges[edge_count] = {v, head[u], capacity, 0};
        head[u] = edge_count++;
        edges[edge_count] = {u, head[v], 0, 0};
        head[v] = edge_count++;
    }

    // BFS 构建分层图
    bool bfs() {
        queue<int> q;
        memset(depth, 0, sizeof(int) * (node_count + 1));

        depth[source] = 1;
        q.push(source);
        while (!q.empty()) {
            int u = q.front();
            q.pop();
            for (int i = head[u]; i != -1; i = edges[i].next) {
                int v = edges[i].to;
                if (!depth[v] && edges[i].capacity > edges[i].flow) {
                    depth[v] = depth[u] + 1;
                    q.push(v);
                }
            }
        }
        return depth[sink];
    }

    // DFS 寻找增广路
    int dfs(int u, int flow) {
        if (u == sink || flow == 0) return flow;

        int total_flow = 0;
        for (int& i = current[u]; i != -1; i = edges[i].next) {
            int v = edges[i].to;
            if (depth[v] == depth[u] + 1) {
                int d = dfs(v, min(flow - total_flow, edges[i].capacity - edges[i].flow));
                if (d > 0) {
                    edges[i].flow += d;
                    edges[i ^ 1].flow -= d;
                    total_flow += d;
                    if (total_flow == flow) return total_flow;
                }
            }
        }
        return total_flow;
    }

    // Dinic 算法求最大流
    void dinic() {
        while (bfs()) {
            memcpy(current, head, sizeof(int) * (node_count + 1));
            max_flow += dfs(source, INF);
        }
    }
} maxFlow;

复杂度

最多循环n - 1轮，每次复杂度为O(MN)，在最坏情况下是O(N2M)

MPM 算法

MPM算法和Danic类似，只是在每一轮寻找增广路不同。

所以这里注重于如何找到找到分层图的一个堵塞流，MPM算法将这个步骤优化：

MPM算法注重于点，引入新关键词点容量，点剩余流量，参考点。

点容量：这个点的所有入边容量之和，这个点的所有出边容量之和，这两者的较小者就是点容量。
点剩余容量：都改成剩余容量。
参考点：每次所有点中剩余容量最小的点。

MPM算法认为，参考点每次都可以因流的增加而使参考点的剩余流量为0。

因为参考点是容量最小的，所以必定有路径能从S到参考点以及参考点到T，把路径删去点的容量，让点的剩余流量为0，再删去点。

总结一遍步骤：

建立剩余流量图的分层图
建立所有点的点容量，找到参考点。
S到参考点以及参考点到T的路径流量删去对应量，使参考点的剩余容量归零0.
删去参考点，找到堵塞流为止回到第二步
在剩余流量图中进行对应的减少以及建立反向边

复杂度

可以发现每次找堵塞流是基于顶点的，所以每次寻找堵塞流的步骤复杂度为O(V2)，应为寻找参考点使用了V 次，而每次删去流也用了V次。这是一个O(V3)的最大流算法

from: https://cp-algorithms.com/graph/mpm.html

struct MPM{
    struct FlowEdge{
        int v, u;
        long long cap, flow;
        FlowEdge(){}
        FlowEdge(int _v, int _u, long long _cap, long long _flow)
            : v(_v), u(_u), cap(_cap), flow(_flow){}
        FlowEdge(int _v, int _u, long long _cap)
            : v(_v), u(_u), cap(_cap), flow(0ll){}
    };
    const long long flow_inf = 1e18;
    vector<FlowEdge> edges;
    vector<char> alive;
    vector<long long> pin, pout;
    vector<list<int> > in, out;
    vector<vector<int> > adj;
    vector<long long> ex;
    int n, m = 0;
    int s, t;
    vector<int> level;
    vector<int> q;
    int qh, qt;
    void resize(int _n){
        n = _n;
        ex.resize(n);
        q.resize(n);
        pin.resize(n);
        pout.resize(n);
        adj.resize(n);
        level.resize(n);
        in.resize(n);
        out.resize(n);
    }
    MPM(){}
    MPM(int _n, int _s, int _t){resize(_n); s = _s; t = _t;}
    void add_edge(int v, int u, long long cap){
        edges.push_back(FlowEdge(v, u, cap));
        edges.push_back(FlowEdge(u, v, 0));
        adj[v].push_back(m);    // 同样添加边的编号
        adj[u].push_back(m + 1);
        m += 2;
    }
    // 创建分层图
    bool bfs(){
        while(qh < qt){
            int v = q[qh++];
            for(int id : adj[v]){
                if(edges[id].cap - edges[id].flow < 1)continue;
                if(level[edges[id].u] != -1)continue;
                level[edges[id].u] = level[v] + 1;
                q[qt++] = edges[id].u;
            }
        }
        return level[t] != -1;
    }
    // 点容量
    long long pot(int v){
        return min(pin[v], pout[v]);
    }
    // 如果已经让一个参考点剩余容量为0，那就再也不用考虑这个点和它的边了
    void remove_node(int v){
        for(int i : in[v]){
            int u = edges[i].v;
            auto it = find(out[u].begin(), out[u].end(), i);
            out[u].erase(it);
            pout[u] -= edges[i].cap - edges[i].flow;
        }
        for(int i : out[v]){
            int u = edges[i].u;
            auto it = find(in[u].begin(), in[u].end(), i);
            in[u].erase(it);
            pin[u] -= edges[i].cap - edges[i].flow;
        }
    }
    // 上述说明的将参考点到起点或者终点删去对应流大小
    void push(int from, int to, long long f, bool forw){
        qh = qt = 0;
        ex.assign(n, 0);
        ex[from] = f;
        q[qt++] = from;
        while(qh < qt){
            int v = q[qh++];
            if(v == to)
                break;
            long long must = ex[v];
            auto it = forw ? out[v].begin() : in[v].begin();
            while(true){
                int u = forw ? edges[*it].u : edges[*it].v;
                long long pushed = min(must, edges[*it].cap - edges[*it].flow);
                if(pushed == 0)break;
                if(forw){
                    pout[v] -= pushed;
                    pin[u] -= pushed;
                }
                else{
                    pin[v] -= pushed;
                    pout[u] -= pushed;
                }
                if(ex[u] == 0)
                    q[qt++] = u;
                ex[u] += pushed;
                edges[*it].flow += pushed;
                edges[(*it)^1].flow -= pushed;
                must -= pushed;
                if(edges[*it].cap - edges[*it].flow == 0){
                    auto jt = it;
                    ++jt;
                    if(forw){
                        in[u].erase(find(in[u].begin(), in[u].end(), *it));
                        out[v].erase(it);
                    }
                    else{
                        out[u].erase(find(out[u].begin(), out[u].end(), *it));
                        in[v].erase(it);
                    }
                    it = jt;
                }
                else break;
                if(!must)break;
            }
        }
    }
    
    // 计算最大流量
    long long flow(){
        long long ans = 0;
        while(true){
            pin.assign(n, 0);
            pout.assign(n, 0);
            level.assign(n, -1);
            alive.assign(n, true);
            level[s] = 0;
            qh = 0; qt = 1;
            q[0] = s;
            if(!bfs())
                break;
            for(int i = 0; i < n; i++){
                out[i].clear();
                in[i].clear();
            }
            for(int i = 0; i < m; i++){
                if(edges[i].cap - edges[i].flow == 0)
                    continue;
                int v = edges[i].v, u = edges[i].u;
                if(level[v] + 1 == level[u] && (level[u] < level[t] || u == t)){
                    in[u].push_back(i);
                    out[v].push_back(i);
                    pin[u] += edges[i].cap - edges[i].flow;
                    pout[v] += edges[i].cap - edges[i].flow;
                }
            }
            pin[s] = pout[t] = flow_inf;
            while(true){
                int v = -1;
                for(int i = 0; i < n; i++){
                    if(!alive[i])continue;
                    if(v == -1 || pot(i) < pot(v))
                        v = i;
                }
                if(v == -1)
                    break;
                if(pot(v) == 0){
                    alive[v] = false;
                    remove_node(v);
                    continue;
                }
                long long f = pot(v);
                ans += f;
                push(v, s, f, false); // 删去参考点到S，true和false表示方向
                push(v, t, f, true);  // 删去参考点到T
                alive[v] = false;
                remove_node(v);       // 删去节点以及对应的边
            }
        }
        return ans;
    }
};

ISAP 算法

通过改进层次图的构建和流量的增广来加速算法的运行。ISAP的关键优化在于使用了一种叫做分层推进（Relabel-to-Front）的技术，以及通过维护每个节点的高度信息来减少不必要的搜索和迭代。

对于原来的ISAP算法，每次删去堵塞流，都需要重新创建分层图，这一步是效率极低的，因为整个遍历了一遍图，实际层数变化的只有其中的几个点。

而变化的原因也是出于找到的流废掉了几条边，使得层数变化，那能否在DFS的过程中就解决这个问题。

显然是在dfs回溯的过程中才能修改层数，所以是从终点往起点进行回溯。那么层数的更新就成为了一个问题，如果从后往前，则最前面的层数只能越来越减少，甚至变成负数，因此我们改变保存层数的方式，保存点到终点的距离，也就是反向的层数。

因此第一步应当是建立初始的层数，也就是对于反图从终点进行BFS标号。

有以下步骤：

对于反图从终点进行BFS标号
当前顶点 V为终点时增广
当前顶点有满足 dist[V] = dist[j] + 1 的出弧时前进
当前顶点无满足条件的出弧时重标号并回退一步。对V点的重标号操作为 dist[i] = 1 + min。min为对于V相邻的dis最小（离汇点最近）的点。
整个循环当源点 s 的距离标号 dist[s] >= n 时结束。

GAP优化

记录每个层的点的数量，gap[i]数组表示距离标号为i的点有多少个。如果重标号使得gap数组中原标号数目变为0，即修改则无法到达，则算法结束。

const int maxn = 1000; // 最大节点数，根据需要调整
const int INF = 1e9; // 无穷大

struct Edge {
    int from, to, cap, flow;

    Edge(int u, int v, int c, int f) : from(u), to(v), cap(c), flow(f) {}
};

bool operator<(const Edge& a, const Edge& b) {
    return a.from < b.from || (a.from == b.from && a.to < b.to);
}

struct ISAP {
    int n, m, s, t;
    vector<Edge> edges;
    vector<int> G[maxn];
    bool vis[maxn];
    int d[maxn];
    int cur[maxn];
    int p[maxn];
    int num[maxn];

    // 添加边
    void AddEdge(int from, int to, int cap) {
        edges.emplace_back(from, to, cap, 0);
        edges.emplace_back(to, from, 0, 0);
        m = edges.size();
        G[from].push_back(m - 2);
        G[to].push_back(m - 1);
    }

    // 广度优先搜索建立层次图
    bool BFS() {
        memset(vis, 0, sizeof(vis));
        queue<int> Q;
        Q.push(t);
        vis[t] = 1;
        d[t] = 0;
        while (!Q.empty()) {
            int x = Q.front();
            Q.pop();
            for (int i : G[x]) {
                Edge& e = edges[i ^ 1];
                if (!vis[e.from] && e.cap > e.flow) {
                    vis[e.from] = 1;
                    d[e.from] = d[x] + 1;
                    Q.push(e.from);
                }
            }
        }
        return vis[s];
    }

    // 初始化
    void init(int n) {
        this->n = n;
        for (int i = 0; i < n; i++) G[i].clear();
        edges.clear();
    }

    // 增广路径
    int Augment() {
        int x = t, a = INF;
        while (x != s) {
            Edge& e = edges[p[x]];
            a = min(a, e.cap - e.flow);
            x = edges[p[x]].from;
        }
        x = t;
        while (x != s) {
            edges[p[x]].flow += a;
            edges[p[x] ^ 1].flow -= a;
            x = edges[p[x]].from;
        }
        return a;
    }

    // 求最大流
    int Maxflow(int s, int t) {
        this->s = s;
        this->t = t;
        int flow = 0;
        BFS();
        memset(num, 0, sizeof(num));
        for (int i = 0; i < n; i++) num[d[i]]++;
        int x = s;
        memset(cur, 0, sizeof(cur));
        while (d[s] < n) {
            if (x == t) {
                flow += Augment();
                x = s;
            }
            bool ok = false;
            for (int& i = cur[x]; i < G[x].size(); i++) {
                Edge& e = edges[G[x][i]];
                if (e.cap > e.flow && d[x] == d[e.to] + 1) {
                    ok = true;
                    p[e.to] = G[x][i];
                    cur[x] = i;
                    x = e.to;
                    break;
                }
            }
            if (!ok) {
                int m = n - 1;
                for (int i : G[x]) {
                    Edge& e = edges[i];
                    if (e.cap > e.flow) m = min(m, d[e.to]);
                }
                if (--num[d[x]] == 0) break;
                num[d[x] = m + 1]++;
                cur[x] = 0;
                if (x != s) x = edges[p[x]].from;
            }
        }
        return flow;
    }
};

Mai Icy