Network Flow(III) - Mincost Maxflow

网络流笔记:最小费用最大流。

Theory

经过每条边,单位流量需支付费用。

最大流量的前提下,费用最小。

满足两个最优化限制的话,经典的思路其实就只有两个了:

  1. 不断选取费用最小的路增广,直到不能增广为止。
  2. 找出一个可行解后,慢慢调整。

步步为营,只维护可行性和最优性其中之一。(后者其实我也不会)

当然说是经典,也就是说还有更现代的。像是松弛算法网络单纯形什么的我既不会,也不是本文探讨重点。不过由于放宽了最优化限制,想必一定很快,但是也很难理解,代码也长。

Classic ZKW Algorithm

考虑不断选取费用最小的路增广,直到不能增广为止。

那其实和普通最大流相比,区别就是一个标号的过程。之前是只有可以走就贪心标号(其实那时叫分层),现在必须走通往终点的最短路。

可以用 KM 修改顶标的方式在每一次增广后修改顶标。

不能直接用于负权图,且在流量不大,费用不小,增广路较长的图表现得不够优秀。因为在零权网络中每次只能增加一条边,重复标号反复尝试增广而次次不能增广。

而且 KM 难写。

Primal-Dual ZKW Algorithm

我们不用 KM 了,不就是用最短路标号嘛,像 Dinic 一样每次用 SPFA 求标号不就好了吗。

牺牲复杂度换取实现简单也是很划算的。

注意一条边如果被走过,那一定比没走过其他边更优。所以只要没断,边权可以视为 $ 0 $ .

LOJ 102 (4126, 772)

#include <cstdio>
#include <queue>
#include <cstring>
#include <algorithm>
using namespace std;
template<class T> void read(T& x)
{
    char c = getchar(); T p = 1, n = 0;
    while(c < '0' || c > '9') {if(c == '-') p = -1; c = getchar();}
    while(c >= '0' && c <= '9'){n = n * 10 + c - '0'; c = getchar();}
    x = p * n;
}
template<class T,class U>void read(T&x,U&y){read(x),read(y);}
template<class T,class U,class V>void read(T&x,U&y,V&z){read(x),read(y),read(z);}
const int maxn = 400 + 5, maxm = 15000 + 5, INF = 0x3f3f3f3f;
struct E{int v, n, f, c, w; } e[maxm << 1];
int h[maxn], ec = 2;
inline void add_edge(int u,int v,int c,int w){e[ec]=(E){v,h[u],0,c,w};h[u]=ec++;e[ec]=(E){u,h[v],0,0,-w};h[v]=ec++;}
deque<int> Q;
bool inq[maxn]; int d[maxn];
int n, m, c, cost;
bool bfs()
{
    memset(inq, 0, sizeof inq); memset(d, INF, sizeof d);
    d[n] = 0; Q.push_back(n); inq[n] = true;
    while(!Q.empty())
    {
        int u = Q.front(); Q.pop_front();
        for(int i = h[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(e[i ^ 1].f < e[i ^ 1].c && d[v] > d[u] - e[i].w)
        {
            d[v] = d[u] - e[i].w;
            if(!inq[v])
            {
                if(!Q.empty() && d[v] < d[Q.front()]) Q.push_front(v); else Q.push_back(v);
                inq[v] = true;
            }
        }
        inq[u] = false;
    }
    for(int i = 2; i < ec; ++i) e[i].w += d[e[i].v] - d[e[i ^ 1].v];
    c += d[1];
    return d[1] != INF;
}
bool vis[maxn]; int cur[maxn];
int dfs(int u = 1, int a = INF)
{
    if(u == n || !a) return cost += a * c, a;
    int flow = 0, f; vis[u] = true;
    for(int&i = cur[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(e[i].f < e[i].c && !e[i].w && !vis[v] && (f = dfs(v, min(a, e[i].c - e[i].f))))
    {
        e[i].f += f, e[i ^ 1].f -= f, flow += f, a -= f;
        if(!a) return flow;
    }
    return flow;
}
int flow()
{
    int ret = 0, f;
    while(bfs()) 
    {
        memcpy(cur, h, sizeof h); 
        do memset(vis, 0, sizeof vis), ret += (f = dfs()); while(f);
    }
    return ret;
}
int main()
{
    read(n, m);
    int s, t, c, w; while(m--) read(s, t), read(c, w), add_edge(s, t, c, w);
    int f = flow();
    printf("%d %d\n", f, cost);
    return 0;
}

Network Flow(V) - Min-Cut

网络流笔记:最小割。

Min-Cut Max-Flow Theorem

这个定理给出求一个图最小割的大小的方法。

Judgement

问题:如何判断一条边是恒在最小割集中、可能在最小割集中、不在最小割集中呢。

对残量网络找 SCC. 设各个点所在的 SCC 编号为 $ id_u $ , 有以下结论成立。

  1. $ id_s \neq id_t $.
  2. $(u,v)$ 能出现在最小割集中 $ \Leftrightarrow id_u \neq id_v $.
  3. $(u,v)$ 必定出现在最小割集中 $ \Leftrightarrow id_u = id_s \wedge id_v = id_t $.

由反证法显然。只需注意到将 SCC 视为一个点后各点之间均满流。

Pairwise Min-Cut

问题:无向图中,求任意两点间最小割。

A Naïve approach: $ \Theta(|V|^2) $

Advanced

显然我们要挖掘最小割的性质。

Theorem 对任意不同的三点 $ a, b, c $ ,有 $ mincut(a,b) \geq \min\{ mincut(a,c),mincut(b,c) \}$. 事实上这三个值总是两个相同的较小值和一个最大值.

Proof. 考虑点 $ a,b $, 且 $ mincut(a,b) $ 为三者最小. 不失一般性假设 $ c $ 在 $ a,b $ 割靠近 $ b $ 一边,则易有 $ mincut(a,b) \ge mincut(a,c) $, 从而 $ mincut(a,b) = mincut(a,c) $ , 对称情况同理。 QED

推论 对任意不同的点,有$ mincut(u,v) \geq \min\{ mincut(u,w_1), mincut(w_1,w_2), \cdots,mincut(w_n,v) \}$.

回忆起最小生成树的贪心正确性的证法,这个形式一模一样。

令 $ (s, t) $ 为 $ (u,v) $ 树上路径的最小边,则 $ mincut(u,v)=mincut(s,t) $. 即两点间最小割为最小割树上路径最小值。

Divide-and-Conquer

选定两个点求最小割,然后将原来的点集分成了两个集合。

每个集合接着求,直到集合只剩小于等于一个点。

BZOJ 4519

无向图中,所有点对最小割共有多少种不同取值。

Sol. 在最小割树上记录边权即可。

Code (1128, 2876)

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <queue>
#include <map>
using namespace std;
template<class T> void read(T& x)
{
    char c = getchar(); T p = 1, n = 0;
    while(c < '0' || c > '9') {if(c == '-') p = -1; c = getchar();}
    while(c >= '0' && c <= '9'){n = n * 10 + c - '0'; c = getchar();}
    x = p * n;
}
template<class T,class U>void read(T&x,U&y){read(x),read(y);}
template<class T,class U,class V>void read(T&x,U&y,V&z){read(x),read(y),read(z);}
const int maxn = 1000 + 10, maxm = 10000 + 10, INF = 0x3f3f3f3f;
struct E{int v, n, f; }e[maxm << 1];
int h[maxn], ec = 2;
inline void add_edge(int u,int v,int c){e[ec]=(E){v,h[u],c};h[u]=ec++;e[ec]=(E){u,h[v],c};h[v]=ec++;}
int d[maxn], s, t;
queue<int> Q;
bool bfs()
{
    memset(d, INF, sizeof d);
    d[s] = 0; Q.push(s);
    while(!Q.empty())
    {
        int u = Q.front(); Q.pop();
        for(int i = h[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(d[v] == INF && e[i].f)
            d[v] = d[u] + 1, Q.push(v);
    }
    return d[t] != INF;
}
int cur[maxn];
int dfs(int u = s, int a = INF)
{
    if(u == t || !a) return a;
    int flow = 0, f;
    for(int &i = cur[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(d[v] == d[u] + 1 && (f = dfs(v, min(a, e[i].f))))
    {
        e[i].f -= f, e[i ^ 1].f += f, a -= f, flow += f;
        if(!a) return flow;
    }
    return flow;
}
bool flag[maxn];
void go(int u)
{
    flag[u] = true;
    for(int i = h[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(!flag[v] && e[i].f)
        go(v);
}
int num[maxn], tmp[maxn], tmp2[maxn];
map<int, bool> m;
int ans = 0;
void solve(int l, int r)
{
    if(l == r) return ;
    for(int i = 2; i < ec; i += 2) e[i].f = e[i ^ 1].f = (e[i].f + e[i ^ 1].f) >> 1;
    int ret = 0;
    s = num[l], t = num[r];
    while(bfs()) 
    memcpy(cur, h, sizeof cur), ret += dfs();
    if(!m[ret]) ++ans, m[ret] = true;
    memset(flag, false, sizeof flag);
    go(s);
    int tc = 0, td = 0;
    for(int i = l; i <= r; ++i) if(flag[num[i]]) tmp[tc++] = num[i]; else tmp2[td++] = num[i];
    for(int i = l; i <= r; ++i) if(i - l < tc) num[i] = tmp[i - l]; else num[i] = tmp2[i - l - tc];
    solve(l, l + tc - 1);
    solve(l + tc, r);
}
int main()
{
    int n, m; read(n, m);
    int u, v, w; while(m--) read(u, v, w), add_edge(u, v, w);
    for(int i = 1; i <= n; ++i) num[i] = i;
    solve(1, n);
    printf("%d\n", ans);
    return 0;
}

BZOJ 2229

无向图中,所有点对最小割不大于 $ x $ 的有多少对,交换不计。

Sol. 求出最小割树后统计两两间答案即可。

也可以像边分那样按边权从小到大枚举搞一发。

Code (1028, 1816)

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <queue>
#include <map>
using namespace std;
template<class T> void read(T& x)
{
    char c = getchar(); T p = 1, n = 0;
    while(c < '0' || c > '9') {if(c == '-') p = -1; c = getchar();}
    while(c >= '0' && c <= '9'){n = n * 10 + c - '0'; c = getchar();}
    x = p * n;
}
template<class T,class U>void read(T&x,U&y){read(x),read(y);}
template<class T,class U,class V>void read(T&x,U&y,V&z){read(x),read(y),read(z);}
const int maxn = 150 + 10, maxm = 3000 + 10, INF = 0x3f3f3f3f;
struct E{int v, n, f;} e[maxm << 1];
int h[maxn], ec = 2;
inline void add_edge(int u,int v,int c){e[ec]=(E){v,h[u],c},h[u]=ec++,e[ec]=(E){u,h[v],c},h[v]=ec++;}
int s, t, n;
queue<int> Q; int d[maxn];
int cur[maxn];
bool bfs()
{
    memset(d, INF, sizeof d);
    d[s] = 0; Q.push(s);
    while(!Q.empty())
    {
        int u = Q.front(); Q.pop();
        for(int i = h[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(e[i].f && d[v] == INF)
            d[v] = d[u] + 1, Q.push(v);
    }
    return d[t] != INF;
}
int dfs(int u = s, int a = INF)
{
    if(u == t || !a) return a;
    int flow = 0, f;
    for(int &i = cur[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(d[v] == d[u] + 1 && (f = dfs(v, min(a, e[i].f))))
    {
        e[i].f -= f, e[i ^ 1].f += f, a -= f, flow += f;
        if(!a) break;
    }
    return flow;
}
bool flag[maxn];
void go(int u)
{
    flag[u] = true;
    for(int i = h[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(!flag[v] && e[i].f)
        go(v);
}
int a[maxn], b[maxn], ans[maxn][maxn];
void dc(int l = 1, int r = n)
{
    if(l >= r) return;
    for(int i = 2; i < ec; i += 2) e[i].f = e[i ^ 1].f = (e[i].f + e[i ^ 1].f) >> 1;
    s = a[l], t = a[r];
    int ret = 0;
    while(bfs()) memcpy(cur, h, sizeof h), ret += dfs();
    memset(flag, false, sizeof flag); go(s);
    for(int i = 1; i <= n; ++i) if(flag[i]) for(int j = 1; j <= n; ++j) if(!flag[j])
        ans[i][j] = ans[j][i] = min(ans[i][j], ret);
    int _l = l, _r = r;
    for(int i = l; i <= r; ++i) if(flag[a[i]]) b[_l++] = a[i]; else b[_r--] = a[i];
    for(int i = l; i <= r; ++i) a[i] = b[i];
    dc(l, _l - 1); dc(_r + 1, r);
}
int main()
{
    int t; read(t);
    while(t--)
    {
        memset(h, 0, sizeof h); ec = 2; memset(ans, INF, sizeof ans);
        int m; read(n, m);
        int u, v, c; while(m--) read(u, v, c), add_edge(u, v, c);
        for(int i = 1; i <= n; ++i) a[i] = i;
        dc();
        int q; read(q);
        while(q--)
        {
            int x, cnt = 0; read(x);
            for(int i = 1; i <= n; ++i) for(int j = 1; j < i; ++j) if(ans[i][j] <= x)
                ++cnt;
            printf("%d\n", cnt);
        }
        puts("");
    }
    return 0;
}

CF 343E

无向图中,给出 $ n $ 个点的一个排列,使得相邻两点作为源汇点的最大流的总和最大。

Sol. 求出最小割树。首先答案肯定不会大于树边和,而注意到树边和总是可以取到。我们考虑选出目前联通块最小的边,它两边分成了两棵子树,我们可以保证让左右子树的路径只经过该边一次,往下递归考虑同理。

Code (62, 2100)

#include <cstdio>
#include <cstring>
#include <queue>
#include <algorithm>
using namespace std;
template<class T> void read(T& x)
{
    char c = getchar(); T p = 1, n = 0;
    while(c < '0' || c > '9') {if(c == '-') p = -1; c = getchar();}
    while(c >= '0' && c <= '9'){n = n * 10 + c - '0'; c = getchar();}
    x = p * n;
}
template<class T,class U>void read(T&x,U&y){read(x),read(y);}
template<class T,class U,class V>void read(T&x,U&y,V&z){read(x),read(y),read(z);}
const int maxn = 200 + 5, maxm = 1000 + 5, INF = 0x3f3f3f3f;
namespace FLOW
{
struct E{int v, n, f;} e[maxm << 1];
int h[maxn], ec = 2;
inline void add_edge(int u,int v,int c){e[ec]=(E){v,h[u],c},h[u]=ec++,e[ec]=(E){u,h[v],c},h[v]=ec++;}
queue<int> Q;
int d[maxn], s, t;
bool bfs()
{
    memset(d, INF, sizeof d);
    d[s] = 0; Q.push(s);
    while(!Q.empty())
    {
        int u = Q.front(); Q.pop();
        for(int i = h[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(e[i].f && d[v] == INF)
            d[v] = d[u] + 1, Q.push(v);
    }
    return d[t] != INF;
}
int cur[maxn];
int dfs(int u = s, int a = INF)
{
    if(u == t || !a) return a;
    int flow = 0, f;
    for(int &i = cur[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(d[v] == d[u] + 1 && (f = dfs(v, min(a, e[i].f))))
    {
        e[i].f -= f, e[i ^ 1].f += f, a -= f, flow += f;
        if(!a) break;
    }
    return flow;
}
}
namespace TREE
{
struct E{int v, n, w, f;}e[maxn << 1];
int h[maxn], ec = 2;
inline void add_edge(int u,int v,int w){e[ec]=(E){v,h[u],w,0},h[u]=ec++,e[ec]=(E){u,h[v],w,0},h[v]=ec++;}
int ti, w;
void dfs(int u, int f = 0)
{
    for(int i = h[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(!e[i].f && v != f)
    {
        dfs(v, u);
        if(e[i].w < w) w = e[i].w, ti = i;
    }
}
void dc(int u = 1)
{
    ti = INF, w = INF;
    dfs(u);
    int i = ti;
    if(i == INF) {printf("%d ", u);return;}
    e[i].f = e[i ^ 1].f = true;
    dc(e[i].v); dc(e[i ^ 1].v);
}
}
bool f[maxn];
void go(int u)
{
    f[u] = true;
    using namespace FLOW;
    for(int i = h[u], v = e[i].v; i; i = e[i].n, v = e[i].v) if(!f[v] && e[i].f)
        go(v);
}
int a[maxn], b[maxn], ret = 0;
void dc(int l, int r)
{
    if(l == r) return;
    using namespace FLOW;
    for(int i = 2; i < ec; i += 2) e[i].f = e[i ^ 1].f = (e[i].f + e[i ^ 1].f) >> 1;
    s = a[l], t = a[r];
    int w = 0; while(bfs()) memcpy(cur, h, sizeof h), w += dfs();
    TREE::add_edge(s, t, w); ret += w;
    memset(f, false, sizeof f); go(s);
    int _l = l, _r = r;
    for(int i = l; i <= r; ++i) if(f[a[i]]) b[_l++] = a[i]; else b[_r--] = a[i];
    for(int i = l; i <= r; ++i) a[i] = b[i];
    dc(l, _l - 1), dc(_r + 1, r);
}
int main()
{
    int n, m; read(n, m);
    int u, v, c; while(m--) read(u, v, c), FLOW::add_edge(u, v, c);
    for(int i = 1; i <= n; ++i) a[i] = i; dc(1, n);
    printf("%d\n", ret);
    TREE::dc();
    return 0;
}

Max Flow(i) - Stoer-Wagner Algorithm

考虑如下问题:求全局最小割。

形式化,求 $ \min_{s,t\in V} mincut(s,t) $ .

First Glimpse

使用最大流算法一般是 $ O(|V|^2|E|) $ 的,而枚举点对则是 $ O(|V|^2) $.

用上最小割树的分析,也需要 $ O(|V|) $.

基于最大流算法的下界似乎就是四次方的了。为了得到更优的解法,需要我们更加关注多一些性质。

Stoer-Wagner

令 $ w(A,t) := \sum_{u\in t} w(u,t) $ .

初始有 $ A = \emptyset $ . 每次选出 $ t : \max w(A,t) $ ,$ A \leftarrow A \cup \{t\} $ ,直到 $ A = V $. 接着合并倒数两个加入的点。

在 $ |V|>1 $ 时,继续循环。

可以证明全局最小割的值一定在最后加进集合的点与集合的函数值中。

这个算法是 $ O(|V|^3) $ .

值得注意的是这个算法与 Prim 求最大生成树的区别。

POJ 2914

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
template<class T> void read(T& x)
{
    char c = getchar(); T p = 1, n = 0;
    while(c < '0' || c > '9') {if(c == '-') p = -1; c = getchar();}
    while(c >= '0' && c <= '9'){n = n * 10 + c - '0'; c = getchar();}
    x = p * n;
}
template<class T,class U>void read(T&x,U&y){read(x),read(y);}
template<class T,class U,class V>void read(T&x,U&y,V&z){read(x),read(y),read(z);}
const int maxn = 500 + 5, INF = 0x3f3f3f3f;
int G[maxn][maxn], cur[maxn], v[maxn];
bool vis[maxn];
int main()
{
    int n, m; 
    while(~scanf("%d%d", &n, &m))
    {
        memset(G, 0, sizeof G);
        for(int u, v, w; m--;) read(u, v, w), G[u][v] += w, G[v][u] += w;
        int ans = INF;
        for(int i = 0; i < n; ++i) v[i] = i;
        while(n > 1)
        {
            int pre = 0;
            memset(vis, 0, sizeof vis);
            memset(cur, 0, sizeof cur);
            for(int i = 1; i < n; ++i)
            {
                int k = -1;
                for(int j = 1; j < n; ++j)
                {
                    if(!vis[v[j]])
                    {
                        cur[v[j]] += G[v[pre]][v[j]];
                        if(k == -1 || cur[v[k]] < cur[v[j]]) k = j;
                    }
                }
                vis[v[k]] = true;
                if(i == n - 1)
                {
                    int s = v[pre], t = v[k];
                    ans = min(ans, cur[t]);
                    for(int j = 0; j < n; ++j) G[s][v[j]] += G[v[j]][t], G[v[j]][s] += G[v[j]][t];
                    v[k] = v[--n];
                }
                pre = k;
            }
        }
        printf("%d\n", ans);
    }
    return 0;
}

Matrix Tree Theorem

矩阵数定理。

【此处应有知识及证明】

线性代数。

其实我根本就不会证233。

SPOJ High

BZOJ 3534

LA 7138