网络单纯形 学习笔记

· · 题解

可能更好的食用体验。

upd：修了英文与标点间多空格的锅。

upd：修了中文与符号间少空格的锅，补上一处句号并自查无锅。这次再格式不过我就不打扰管理了 orz。

网络单纯形是一种神奇的算法。它可以求解带负圈的费用流，可以过 HLPP 板子，但它的（最坏）复杂度好像是指数级，尽管我并不会证。

感性理解：它和线规算法 simplex 有许多相似之处，而 simplex（最坏）是指数级的。

虽然但是，据 CF（见本文末尾）上所讲，它的平均时间复杂度是 O(VE)，且常数较小（无 LCT 情况下）。SPFA被卡前：这我熟。

网络单纯形算法可用来求解无源汇有负圈上下界最小费用可行流。好像常见的网络流都可以转化为这个，但加入上下界后常数会略大，不知道是不是我写炸了，欢迎大佬修改。

将有源汇最小费用最大流转换为无源汇最小费用可行流的方法：从汇向源连一条费用为 -\infty，流量为 +\infty 的边。实际可取 inf 为其余边费用/流量之和加一。

将上下界最小费用可行流转换为一般的最小费用可行流的方法：先所有边强制推下界的流，记录“溢出”的流，建虚拟点 K，K 向所有“上溢”（流量有多余）的点连边，权为 -\infty；所有“下溢”的点向 K 连边，权为 0。最后检查这些边是否跑满即可。

感性理解：若存在一条路径，从某个“上溢”的点 A 到另一个“下溢”的点 B，则它与 B \rightarrow K \rightarrow A 的路径构成负环，可推流（见下文），这相当于之前我们强制推的流从这条路径回去。

网络单纯形算法的大致思路：求解最小费用可行流时，给定图 G=(V,E)，在算法过程中，我们维护它的生成森林（原图可能不联通）边集 T，每次找一条不在 T 内的边 t \in E \setminus T，若 T \cup \{t\} 包含负环，则在 T 中沿此负环推流，选出一条满流的边删去，将 t 加入 T。重复此过程直到 T 中不含负环，我们就得到了 G 的最小费用可行流。这样的思想使得网络单纯型算法天生就支持带负环的费用流。

线规角度的解释：生成森林相当于初始可行解，向负环推流相当于转动变量。

判断一条边加入生成树后是否有负环

考虑函数 h:V \rightarrow \mathbf{R}，满足对 \forall e \in T, h(e.\text{to}) - h(e.\text{from}) = e.\text{cost}，那么对于两个点 a,b，h(a)-h(b) 就是 T 上 a 到 b 一条路的权值之和。对边 e \in E，若 h(e.\text{to}) - h(e.\text{from}) + e.\text{cost} < 0，则它所在的环路（沿着它的方向）为负环。这样可做到 \Theta(1) 判断是否这条边加入后存在负环。

找到负环后推流

先从 e.\text{from} 和 e.\text{to} 向上跳到 LCA，再从两边分别找能推的最小流值，最后推流。

这里如果用 LCT 维护，则可从 O(n) 优化到 O(\log n)，但常数大。

推流后，别忘了删边并反转被删的链的父子关系，保证生成树还是树。

Code

跑得飞快.jpg

#include <cstdio>
#include <algorithm>
#include <vector>
#define sd std::
#define UP(i,s,e) for(auto i=s; i!=e; ++i)
#define IT(i,x) for(int i=flow::head[x]; i!=flow::EDGE_NIL; i=flow::es[i].nxt)
namespace flow{ // }{{{
typedef long long ll;
constexpr int V = 5e3+100, E = V+5e4+100;
constexpr int EDGE_NIL = 0;
struct Edge{
        int to;
        ll lf, cost;
        int nxt;
} es[E*2+4];
ll sumcost = 0, sumflow = 0;
int is, it, iv;
ll minc, maxf;
int head[V], cnt = (EDGE_NIL|1)+1;
ll pi[V]; // h 函数
int fe[V], mark[V], time = 2; // fe: father edge
int fa[V];
void init(int v, int s, int t){
        is = s, it = t, iv = v+1;
        sd fill(head, head+iv, EDGE_NIL);
        sd fill(mark, mark+iv, 0);
        sd fill(fe, fe+iv, EDGE_NIL);
        time = 2;
        cnt = (EDGE_NIL|1)+1;
        sumcost = sumflow = 0;
        minc = maxf = 0;
}
void addflow(int s, int t, ll f, ll c){
        es[cnt] = (Edge){t, f, c, head[s]}, head[s] = cnt++;
        es[cnt] = (Edge){s, 0, -c, head[t]}, head[t] = cnt++;
        sumflow += f, sumcost += sd abs(c);
}
void mktree(int x, int from_e){
        fe[x] = from_e;
        fa[x] = es[from_e^1].to;
        mark[x] = 1;
        for(int i=head[x]; i!=EDGE_NIL; i=es[i].nxt){
                if(mark[es[i].to] == 1 || es[i].lf == 0) continue;
                mktree(es[i].to, i);
        }
}
ll getpi(int x){ // 获取某个点的 h 值
        if(mark[x] == time) return pi[x];
        mark[x] = time;
        pi[x] = getpi(fa[x]) - es[fe[x]].cost;
        return pi[x];
}
ll pushflow(int e){ // 返回减少的费用
        int rt = es[e].to, lca = es[e^1].to;
        time++;
        while(rt){ // rt 用来标记点
                mark[rt] = time;
                rt = fa[rt];
        }
        while(mark[lca] != time){ // lca 用来找 LCA
                mark[lca] = time;
                lca = fa[lca];
        }
        ll df = es[e].lf; // df 为流量的改变量
        int todel = e, dir = -1; // dir: direction, dir 为 0 表示 es[e].to 方向
        for(int i=es[e^1].to; i!=lca; i=fa[i]){ // 两边向上找能推的最小流
                if(es[fe[i]].lf < df){
                        df = es[fe[i]].lf;
                        todel = fe[i];
                        dir = 1;
                }
        }
        for(int i=es[e].to; i!=lca; i=fa[i]){
                if(es[fe[i]^1].lf < df){
                        df = es[fe[i]^1].lf;
                        todel = fe[i];
                        dir = 0;
                }
        }
        ll dcst = 0; // delta cost
        if(df) { // 推流
                for(int i=es[e].to; i!=lca; i=fa[i]){
                        es[fe[i]].lf += df;
                        es[fe[i]^1].lf -= df;
                        dcst += es[fe[i]^1].cost * df;
                }
                for(int i=es[e^1].to; i!=lca; i=fa[i]){
                        es[fe[i]].lf -= df;
                        es[fe[i]^1].lf += df;
                        dcst += es[fe[i]].cost * df;
                }
                es[e].lf -= df;
                es[e^1].lf += df;
                dcst += es[e].cost * df;
        }
        if(todel == e) return dcst;
        int last = e^dir, lastu = es[e^dir^1].to;
        for(int i=es[e^dir].to; i!=es[todel^1].to; ){
                mark[i]=time-1;
                int i_ = fa[i];
                fa[i] = lastu;
                lastu = i;
                sd swap(fe[i], last);
                last ^= 1;
                i=i_;
        }
        return dcst;
}
void mcmf(){
        ll sfl_ = sumflow, scs_ = sumcost;
        addflow(iv-1, is, sfl_, 0);
        addflow(it, iv-1, sfl_, -scs_-1);
        sumflow = sfl_, sumcost = scs_;
        mktree(iv-1, EDGE_NIL);
        mark[iv-1] = ++time;
        fa[iv-1] = 0;
        bool run = true;
        while(run){
                run = false;
                UP(i, (EDGE_NIL|1)+1, cnt){
                        int s = es[i^1].to, t = es[i].to;
                        if(es[i].lf && es[i].cost + getpi(t) - getpi(s) < 0){
                                run = true;
                                minc += pushflow(i);
                        }
                } // 一圈下来没流可推时算法结束
        }
        maxf = es[cnt-1].lf;
        minc += maxf * (scs_+1);
}
} // {}}}
namespace m { // }{{{
int in, im, is, it;
void work() {
        scanf("%d%d%d%d", &in, &im, &is, &it);
        is--, it--;
        flow::init(in, is, it);
        UP(i, 0, im) {
                int u, v, w, c;
                scanf("%d%d%d%d", &u, &v, &w, &c);
                u--, v--;
                flow::addflow(u, v, w, c);
        }
        flow::mcmf();
        printf("%lld %lld\n", flow::maxf, flow::minc);
}
} // {}}}
int main() { m::work(); return 0; }

参考资料