1

算法学习笔记(8.1): 网络最大流算法 EK, Dinic, ISAP - jeefies

 1 year ago
source link: https://www.cnblogs.com/jeefy/p/17052660.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

网络最大流

前置知识以及更多芝士参考下述链接
网络流合集链接:网络流


最大流,值得是在不超过管道(边)容量的情况下从源点到汇点最多能到达的流量

抽象一点:使 ∑(S,v)∈Ef(S,v) 最大的流函数被称为网络的最大流,此时的流量被称为网络的最大流量


有了最大流量,就可以通过奇奇怪怪的建模解决很多令人摸不着头脑的题

例如二分图

对于一张二分图,经过建模之后我们可以这样画

其中左部点集 A={1,2,3,4},右部点集 B={5,6,7,8}, 其中源点为 0,汇点为 9

建模过程:

新增一个源点 S 和一个汇点 T, 从 S 到每一个左部点连有向边,从每一个右部点到 T 连有向边,把原二分图的每条边看作从左部点到右部点的有向边,形成了一张 n+2 个点 n+m 条边的网络。其中每一条边的容量都为 1。

不难发现,二分图的最大匹配数就等于网络的最大流量。求出最大流后,所有有有”流“经过的点,边就是匹配点,匹配边。

进一步的:如果要求二分图多重匹配,依据题目信息改变连接汇点和源点的边的容量即可

计算最大流的算法很多,这里主要讲解 EK(Edmonds−Karp) ,Dinic 和 ISAP 算法。


EK 增广路算法

这里的增广路与二分图里面的增广路有不一样了 Q^Q

增广路:若一条源点 S 到 T 的路径上各边的剩余容量都严格大于 0,则称这条路径为增广路。显然,可以让一个流沿着增广路流过,使得网络的流量增大。

而EK的思路就是不断进行广度优先搜索寻找增广路,直到不存在增广路为止。

而在搜索的时候,我们只考虑图中所有 f(x,y)<c(x,y) 的边(或者说是有剩余容量的边)。在寻找路径的同时,还要记录其前驱结点以及路径上的最小容量minf, 在找到一条增广路后,则网络的流量可以增加 minf。

但是,考虑到斜对称的性质,由于我们需要把增广路上的所有边的剩余容量 e(u,v) 减去minf,所以要在对其反边容量 e(v,u) 加上 minf。

初始化的时候,如果是有向边 (u,v),则 e(u,v)=c(u,v),e(v,u)=0,如果是无向边,则 e(u,v)=e(v,u)=c(u,v)。

?: 为什么会出现无向边,网络流不是有向图吗?

考虑双向道路马路,既可以顺流,又可以逆着。

例如:[ICPC-Beijing 2006] 狼抓兔子 - 洛谷

这道题需要用到最小割,顺便说一下,最小割 = 最大流

?: 为什么使用BFS,而不是DFS?

因为DFS可能会绕圈圈……在讲述DInic的时候我会再提及

复杂度:复杂度上界为 O(nm2),然而实际上远远达不到这个上界,效率还行,可以处理 103∼104 规模的网络

--《算法竞赛进阶指南》

我不会证明,下面的两个算法也不会 Q^Q

这里给出一种参考代码

【模板】网络最大流 - 洛谷

提交记录:记录详情

#include <iostream>
#include <cstring>
#include <algorithm>
#include <deque>

using std::deque;

const int N = 2e3 + 7, M = 5e5 + 7, INF = 0x7F7F7F7F;

int n, m, s, t;
int to[M], nex[M], wi[M] = {INF};
int head[N], tot = 1;

void add(int u, int v, int w) {
    to[++tot] = v, nex[tot] = head[u], wi[tot] = w, head[u] = tot;
}

void read() {
    scanf("%d %d %d %d", &n, &m, &s, &t);

    int u, v, w;
    for (int i = 0; i < m; ++i) {
        scanf("%d %d %d", &u, &v, &w);
        add(u, v, w);
        add(v, u, 0);
    }
}

#define min(x, y) ((x)<(y)?(x):(y))
#define pop() que.pop_front()
#define top() que.front()
#define push(x) que.push_back(x);
#define empty() que.empty()

static int inq[N], it = 0;
// px记录上一个点,pe记录遍历过来的边
static int px[N], pe[N];
inline int bfs() {
    deque<int> que;

    push(s); inq[s] = ++it;

    int x, y;
    while (!empty()) {
        x = top(); pop();
        for (int i = head[x]; i; i = nex[i]) {
            if ((inq[(y = to[i])] ^ it) && wi[i]) {
                px[y] = x, pe[y] = i;
                if (y == t) return 1; // 找到增广路了,
                inq[y] = it; push(y);
            }
        }
    }
    return 0; // 到不到了,没有增广路了
}

void work(long long & res) {
    while (bfs()) {
        int val = INF;
        for (int x = t; x ^ s; x = px[x]) {
            val = min(val, wi[pe[x]]);
        }

        for (int x = t; x ^ s; x = px[x]) {
            wi[pe[x]] -= val;
            // 处理反边的时候利用了成对变换的方法!
            wi[pe[x] ^ 1] += val;
        }

        res += val;
    }
}

int main() {
    read();

    long long res = 0;
    work(res);
    printf("%lld\n", res);
    return 0;
}

Dinic

考虑到 EK 算法每一次在残量网络上只找出来的一条增广路,太慢了,所以有了更优化的东西 Dinic?歌姬吧

先引入一点点概念:

深度:在搜索树上的深度(BFS搜索时的层数)

残量网络:网络中所有节点以及剩余容量大于 0 的边构成的子图

分层图:依据深度分层的一段段图……或者说在残量网络上,所有满足 dep[u]+1=dep[v] 的边 (u,v) 构成的子图。

分层图显然是一张有向无环图

Dinic 算法不断重复下述过程,直到在残量网络中,S 不能到达 T

  • 利用BFS求出分层图

  • 在分层图上DFS寻找增广路,在回溯的时候实时更新剩余容量。另外,每个点可以同时流出到多个结点,每个点也可以接收多个点的流。

?: 这里为什么可以使用DFS

由于我们分了层,意味着DFS只会向更深的地方搜索,而不会在同一层乱跳,甚至搜索到前面。这也是为什么EK用BFS更优秀

复杂度:一般来说,时间复杂度为 O(n2m),可以说是不仅简单,而且容易实现的高翔算法之一,一般能够处理 104∼105 规模的网络。特别的,用此算法求二分图的最大匹配时只需要 O(mn), 实际上表现会更好。

没有当前弧优化:提交详情

有当前弧优化:记录详情

// 重复内容已省略

int dis[N], vis[N], vt = 0;
int now[N]; // 用于当前弧优化
// return true if exists non-0 road to t
bool bfs() {
    memset(dis, 0, sizeof(dis)); dis[s] = 1;

    deque<int> que;
    que.push_back(s);
    while (que.size()) {
        int x = que.front(); que.pop_front();
        now[x] = head[x]; // 更新当前弧
        for (int y, i = head[x]; i; i = nex[i]) {
            if (!dis[y = to[i]] && wi[i]) {
                dis[y] = dis[x] + 1;
                que.push_back(y);
                if (y == t) return true;
            }
        }
    }
    return false;
}

#define min(x, y) ((x) < (y) ? (x) : (y))

long long dinic(int x, long long maxflow) {
    if (x == t) return maxflow;
    long long rest = maxflow, k;
    for (int y, i = now[x]; i && rest; i = nex[i]) {
        now[x] = i; // 更新当前弧
        // 要在更深的一层,以及需要有剩余流量
        if (dis[y = to[i]] == dis[x] + 1 && wi[i]) {
            k = dinic(y, min(rest, wi[i]));
            if (!k) dis[y] = 0;
            wi[i] -= k, wi[i ^ 1] += k;
            rest -= k;
        }
    }
    return maxflow - rest;
}

int main() {
    read();
    long long maxflow = 0, flow;
    while (bfs()) {
        while (flow = dinic(s, INF)) maxflow += flow;
    } 
    printf("%lld\n", maxflow);
}

?: 当前弧优化是个啥玩意

注意到如果我们每一次遍历后,对于当前边 (u,v),不可能再有流量流过这条边,所以我们可以暂时的删除这条边……注意,只是暂时,每一分层的时候是需要考虑这条边的,因为这条边的剩余流量不一定为 0


某位大佬的博客上说这是究极最大流算法之一。还有一个HLPP(最高标记预留推进),思路完全与这几个方法不同,不依赖于增广路,我会把它放在另外的文章中单独讲。我可不会告诉你们是我不会优化,太笨了,看不懂大佬的优化

题目链接:【模板】最大流 加强版 / 预流推进 - 洛谷

这是我的:记录详情 4.77s

这是大佬的:记录详情 185ms

由于Dinic需要多次BFS……所以有些不满足的数学家决定优化常数……于是有了ISAP,只需要一次BFS的东西……

可恶,竟然没有找到不用gap优化的写法 T^T

ISAP算法从某种程度上是SAP算法和Dinic的融合

SAP算法就是所谓的EK算法……ISAP也就是Improved SAP……但是主体怎么跟DInic几乎一模一样!

算法流程如下:

  1. 从 T 开始进行BFS,直到 S ,标记深度,同时记录当前深度有多少个

  2. 利用DFS不断寻找增广路,思路与Dinic类似

  3. 每次回溯结束后,将所在节点深度加一(远离 T 一点),同时更新深度记录。如果出现了断层(有一个深度没有点了)那么结束寻找。

?: 为什么需要深度加一

由于我们在便利过一次过后,这个点不可能再向更靠近 T 的点送出流量,所以只能退而求其次,给自己同层的结点送流量。

怎么跟Dinic一摸一样啊,关键是也可以用当前弧优化,只是我用写的是vetor存图……用不了

参考代码……

提交题目还是【模板】网络最大流 - 洛谷

记录详情

!! 竟然在最优解第二页 O-O

对于下面代码做出一些解释

?: 为什么终止条件是 dep[s] > n

考虑如果是dep[s] <= n 的情况,由于只有n个点,意味着只要最大深度 <n 那么要么是有连续的层,要么断层了(此时我们在DFS中会将dep[s]设为n+1

如果最大深度 >n 所以一定会有一个深度是没有点的,意味着一定出现了断层,也就是流量无法到达了

所以,更新答案之后就可以结束循环了

// 写这个的时候,借鉴了写HLPP最优解的大佬写快读的方法……
template<typename T>
inline void read(T &x) {
    char c, f(0); x = 0;
    do if ((c = getchar()) == '-') f = true; while (isspace(c));
    do x = (x<<3) + (x<<1) + (c ^ 48), c = getchar(); while (isdigit(c));
    if (f) x = -x;
}
template <typename T, typename ...Args> inline void read(T &t, Args&... args) { read(t), read(args...); }

typedef long long Data;
using namespace std;

const int N = 207, M = 5007;

struct Edge {
    int to;
    size_t rev; // 反边的位置,用int也没问题
    Data flow;
    Edge(int to, size_t rev, Data f) : to(to), rev(rev), flow(f) {}
};

class ISAP {
public:
    int n, m, s, t;
    vector<int> dep;
    int q[N * 2], gap[N * 2];

    // vector< vector<Edge> > v;
    vector<Edge> v[N * 2];

    ISAP(int n, int m, int s, int t) : n(n), m(m), s(s), t(t) {
        input();
    } 

    inline void input() {
        // v.resize(n + 1);
        for (int x, y, f, i(0); i ^ m; ++i) {
            read(x, y, f);
            v[x].push_back(Edge(y, v[y].size(), f));
            v[y].push_back(Edge(x, v[x].size() - 1, 0));
        }
    }

    inline void init() {
        dep.assign(n + 1, -1);
        dep[t] = 0, gap[0] = 1;

        // 如果要用手写队列,要开大一点……避免玄学RE,虽然理论上N就够了
        register int qt(0), qf(0);
        q[qt++] = t;
        int x, y;
        while (qf ^ qt) {
            x = q[qf++];
            for (auto &e : v[x]) {
                if (dep[(y = e.to)] == -1) // if dep[y] != -1
                    ++gap[(dep[y] = dep[x] + 1)], q[qt++] = y;
            }
        } // bfs end
    }

    inline Data sap(int x, Data flow) {
        if (x == t) return flow;

        Data rest = flow;
        int y, f;
        for (auto &e : v[x]) {
            if (dep[(y = e.to)] + 1 == dep[x] && e.flow) {
                f = sap(y, min(e.flow, rest));
                if (f) {
                    e.flow -= f, v[e.to][e.rev].flow += f;
                    rest -= f;
                }
                if (!rest) return flow; // flow all used
            }
        }

        // change dep
        if (--gap[dep[x]] == 0) dep[s] = n + 1; // can not reach to t
        ++gap[++dep[x]]; // ++depth
        return flow - rest;
    }

    inline Data calc() {
        Data maxflow(0);
        static const Data INF(numeric_limits<Data>::max());
        // dep[s]最大为n,为一条链的时候
        while (dep[s] <= n) {
            // 如果要当前弧优化,在这里需要重置当前弧的now!
            maxflow += sap(s, INF);
        }
        return maxflow;
    }
};

int main() {
    int n, m, s, t;
    read(n, m, s, t);

    static ISAP isap(n, m, s, t);
    isap.init();
    printf("%lld\n", isap.calc());

    return 0;
}

?: 如果我想用vector存图实现当前弧优化怎么整

在sap函数的主体部分

for (int & i = now[x]; i < G[x].size(); ++i) {
    Edge & e = G[x][i];
    if (dep[(y = e.to)] + 1 == dep[x] && e.flow) {
        f = sap(y, min(rest, e.flow));
        if (f) {
            rest -= f, e.flow -= f;
            G[e.to][e.rev].flow += f;
        }
        if (!rest) return flow;
    }
}

在calc不部分

while (dep[s] <= n) {
    now.assign(n, 0);
    maxflow += sap(s, INF);
}

然后……就搞定了QwQ

作者有话说

一般来说,如果图非常稠密(边数远远大于点数),当前弧优化的力度就非常大了

如:Zoj3229 Shoot the Bullet|东方文花帖|【模板】有源汇上下界最大流 - 洛谷

这个专题我会放在网络流的其他部分详解,敬请期待……

写了当前弧优化的Dinic能轻松过……没写全TLE

虽然没写当前弧优化的ISAP能更快的过前三个点,但最后一个点过不了……QwQ

我没有试过当前弧优化的ISAP

更新:有当前弧优化的ISAP可以过

但是如果边数不多,当前弧优化可能就成了负优化了……所以需要根据题目数据合理使用


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK