网络流6：费用流基础
网络流6 费用流基础（含消圈和原始对偶）
网络流3、4、5分别是最大流应用、最小割及其简单应用和最大流与最小割转化，但是我懒了还没写……
定义
（前置知识：网络流本体）
在费用流中，每条边除了容量外，还有一个“单位流量费用” $w(u,v) $。当 \((u,v)\) 的流量为 \(f(u,v)\) 时，其价值为 $f(u,v)\times w(u,v) $。
同前面对流量的扩展定义，价值也有这种性质，即 \(w(u,v)=-w(v,u)\)。这样一条边不论以谁为起点，都可以正确的计算价值。
做法
存在一个贪心的算法．它的思路是每次寻找单位费用最小的增广路进行增广，直到图上不存在增广路为止。可以用 spfa 实现。
但是，如果存在负圈，我们的 spfa 就会暴毙。此时，需要先使用消圈算法消去图上的负圈。
老样子证明懒了，详见 OI wiki。
EK算法求费用流
dinic 算法在权值值域较大且均匀分布的情况下，是比 EK 要慢的，所以这里主要讲 EK。
本质就是将原本 EK 算法的增广路 dfs 函数改成了 spfa。
下面是最小费用最大流，对边权取负就可以当最大费用最大流
namespace Flow{
#define Tpl int
#define Tpw int
struct Edge{
int v,nx;Tpl c;Tpw w;
}e[M*2];
int hd[N],ecnt;
void add(int u,int v,Tpl c,Tpw w){
e[++ecnt]={v,hd[u],c,w};hd[u]=ecnt;
e[++ecnt]={u,hd[v],0,-w};hd[v]=ecnt;
}
void init(){ecnt=-1;memset(hd,-1,sizeof(hd));}
int S,T;
Tpw dis[N];
int pre[N];
Tpl incf[N];
queue <int> q;
bool vis[N];
bool spfa(){
while(!q.empty()) q.pop();
memset(dis,0x3f,sizeof(dis));
memset(incf,0,sizeof(incf));
memset(vis,0,sizeof(vis));
incf[S]=INF,dis[S]=0,q.push(S);
while(!q.empty()){
int u=q.front();
q.pop();
vis[u]=0;
for(int i=hd[u];~i;i=e[i].nx){
int v=e[i].v;Tpw &w=e[i].w;Tpl &c=e[i].c;
if(!c||dis[v]<=dis[u]+w) continue;
dis[v]=dis[u]+w;
incf[v]=min(incf[u],c);
pre[v]=i;
if(!vis[v]) q.push(v),vis[v]=1;
}
}
return incf[T]>0;
}
pair<Tpl,Tpw> EK(int s,int t){
S=s,T=t;
Tpl flow=0;Tpw cost=0;
while(spfa()){
Tpl t=incf[T];
flow+=t,cost+=(Tpw)t*dis[T];
for(int i=T;i!=S;i=e[pre[i]^1].v){
e[pre[i]].c-=t;
e[pre[i]^1].c+=t;
}
}
return {flow,cost};
}
int jiliguala=(init(),1);
#undef Tpl
#undef Tpw
}using Flow::init;using Flow::add;using Flow::EK;
最短路优化扩展&消圈算法
有些时候，我们会遇到负圈，由于增广路给予最短路，我们就会挂掉。
为了解决这种情况，我们可以先试着将负费用边流满，转化为有源汇上下界费用流。由于贪心，开局流满负圈边的情况下，我们不会重新造出来一个负圈，因此可以解决负圈的问题。
也顺便在此给出上下界费用流&消圈算法的模板（含 EK 本体和建图）
namespace Flow{
#define Tpl int
#define Tpw int
struct Edge{
int v,nx;Tpl c;Tpw w;
}e[M*2];
int hd[N],ecnt;
void add(int u,int v,Tpl c,Tpw w){
e[++ecnt]={v,hd[u],c,w};hd[u]=ecnt;
e[++ecnt]={u,hd[v],0,-w};hd[v]=ecnt;
}
void init(){ecnt=-1;memset(hd,-1,sizeof(hd));}
Tpl del(){Tpl c=e[ecnt].c;e[ecnt].c=e[ecnt^1].c=0;return c;}
int S,T;
Tpw dis[N];
int pre[N];
Tpl incf[N];
queue <int> q;
bool vis[N];
bool spfa(){
while(!q.empty()) q.pop();
memset(dis,0x3f,sizeof(dis));
memset(incf,0,sizeof(incf));
memset(vis,0,sizeof(vis));
incf[S]=INF,dis[S]=0,q.push(S);
while(!q.empty()){
int u=q.front();
q.pop();
vis[u]=0;
for(int i=hd[u];~i;i=e[i].nx){
int v=e[i].v;Tpw &w=e[i].w;Tpl &c=e[i].c;
if(!c||dis[v]<=dis[u]+w) continue;
dis[v]=dis[u]+w;
incf[v]=min(incf[u],c);
pre[v]=i;
if(!vis[v]) q.push(v),vis[v]=1;
}
}
return incf[T]>0;
}
pair<Tpl,Tpw> EK(int s,int t){
S=s,T=t;
Tpl flow=0;Tpw cost=0;
while(spfa()){
Tpl t=incf[T];
flow+=t,cost+=(Tpw)t*dis[T];
for(int i=T;i!=S;i=e[pre[i]^1].v){
e[pre[i]].c-=t;
e[pre[i]^1].c+=t;
}
}
return {flow,cost};
}
int jiliguala=(init(),1);
#undef Tpl
#undef Tpw
}using Flow::init;using Flow::add;using Flow::EK;using Flow::del;
int W[N];
void mian()
{
int s,t;
read(n);read(m);read(s),read(t);
int S=0,T=n+1;
int resf=0,resw=0;
while(m--){
int u,v,w,c;
read(u),read(v),read(c),read(w);
if(w>0) add(u,v,c,w);
else{
add(v,u,c,-w);
W[u]+=c,W[v]-=c;
resw+=w*c;
}
}
For(i,1,n){
if(W[i]<0) add(S,i,-W[i],0);
else add(i,T,W[i],0);
}
add(t,s,INF,0);
auto p=EK(S,T);
resf=del(),resw+=p.second;
p=EK(s,t);
resf+=p.first,resw+=p.second;
printf("%d %d\n",resf,resw);
}
Primal-Dual 原始对偶算法
咕咕咕。 开始干活！
由于无良出题人会闲没事干发癫，可能会去卡 spfa。
关于 spfa，它死了。
为了解决这种情况，我们想要更换最短路算法。我们的第一反应肯定是换成 Dijkstra（dij） 算法。但网络上存在单位费用为负的边，因此无法直接使用 dij 算法。由此，我们引入了，万众瞩目的——Johnson 全源最短路（的思路）！
先简单介绍一下 Johnson 全源最短路的思路。我们随便找一个点作为源点，跑一个 Bellman-Ford（等能处理负圈边的）最短路，得到每个点的最短路 \(h_i\)，然后将每个点的势能设置为 \(h_i\)，并将边 \((u,v)\) 的边权设置为 $ w(u,v)'\leftarrow w(u,v)+h_u-h_v $。根据 $h_v\le h_u+d(u,v) $ 的性质，得知这样操作后边权为正。
假如说有一条最短路 $p_1\rightarrow p_2\rightarrow p_3\rightarrow p_4 \rightarrow \dots \rightarrow p_n $，那么其修改边权后的最短路长度为 $ d'(p_1, p_n) = (w(p_1, p_2) + h_{p_1} - {p_2}) + (w(p_2, p_3) + h_{p_2} - h_{p_3}) + (w(p_3, p_4) + h_{p_3} - h_{p_4}) + \dots + (w(p_{n-1}, p_n) + h_{p_{n-1}} - h_{p_n}) = d(p_1,p_n) + h_{p_1} - h_{p_n} $，也就是只与两点的势能有关。复杂度 $ O(nm \log m) $
我们通过为每个点设计一个势能，将网络上所有边的费用（下面简称为边权）全部变为非负值，从而用 dij 计算最短路。问题在于，每次增广后图的形态会发生变化，因此势能需要动态更新。
直接给出解决方案：每次跑完 dij 后，使用新的最短路 \(d_i\) 更新势能：$ h_i \leftarrow h_i+d_i $，这显然也能代表最短路，问题在于如何证明其边权仍为真。
证明类似：
- 对于原有的边，增广前，有 $ d_i+(w(i,j)+h_i-h_j)-d_j \ge 0 $，因此 $w(i,j) + (h_i+d_i) - (h_j+d_j) \ge 0 $
- 对于反向边，如果原本就存在则和上面相同，而新的边肯定来自于增广，有 $ d_i+(w(i,j)+h_i-h_j)-d_j = 0 $（否则不会在增广路上），和上文类似简单变形后能得出 $ d_i+(w(i,j)+h_i-h_j)-d_j = 0 $
综上，边权非负。
至此，我们学会了 Primal-Dual 原始对偶算法的所有内容，复杂度 $O(nm + fm \log m) $ 。如果图有特殊性质，第一次的最短路可以换成 bfs 之类的。
namespace Flow{
#define Tpl int
#define Tpw int
struct Edge{
int v,nx;Tpl c;Tpw w;
}e[M*2];
int hd[N],ecnt;
void add(int u,int v,Tpl c,Tpw w){
e[++ecnt]={v,hd[u],c,w};hd[u]=ecnt;
e[++ecnt]={u,hd[v],0,-w};hd[v]=ecnt;
}
void init(){ecnt=-1;memset(hd,-1,sizeof(hd));}
// Tpl del(){Tpl c=e[ecnt].c;e[ecnt].c=e[ecnt^1].c=0;return c;}//版子题就没有删边操作了
int S,T;
Tpw dis[N],h[N];
int pre[N];
Tpl incf[N];
bool vis[N];
queue <int> q;
void spfa(){
while(!q.empty()) q.pop();
memset(dis,0x3f,sizeof(dis));
memset(incf,0,sizeof(incf));
memset(vis,0,sizeof(vis));
incf[S]=INF,dis[S]=0,q.push(S);
while(!q.empty()){
int u=q.front();
q.pop();
vis[u]=0;
for(int i=hd[u];~i;i=e[i].nx){
int v=e[i].v;Tpw &w=e[i].w;Tpl &c=e[i].c;
if(!c||dis[v]<=dis[u]+w) continue;
dis[v]=dis[u]+w;
incf[v]=min(incf[u],c);
pre[v]=i;
if(!vis[v]) q.push(v),vis[v]=1;
}
}
}
priority_queue <pair <Tpw,int> > pq;
bool dij(){
while(!pq.empty()) pq.pop();
memset(dis,0x3f,sizeof(dis));
memset(incf,0,sizeof(incf));
memset(vis,0,sizeof(vis));
incf[S]=INF,dis[S]=0,pq.push({0,S});
while(!pq.empty()){
int u=pq.top().second;pq.pop();
if(vis[u]) continue;vis[u]=1;
for(int i=hd[u];~i;i=e[i].nx){
int v=e[i].v;Tpw w=e[i].w+h[u]-h[v];Tpl &c=e[i].c;
if(!c||dis[v]<=dis[u]+w) continue;
dis[v]=dis[u]+w;
incf[v]=min(incf[u],c);
pre[v]=i;
pq.push({-dis[v],v});
}
}
return incf[T]>0;
}
pair<Tpl,Tpw> Primal_Dual(int l,int r,int s,int t){
S=s,T=t;
spfa();
For(i,l,r) h[i]=dis[i];
Tpl flow=0;Tpw cost=0;
while(dij()){
Tpl t=incf[T];
flow+=t,cost+=(Tpw)t*(dis[T]-h[S]+h[T]);
for(int i=T;i!=S;i=e[pre[i]^1].v){
e[pre[i]].c-=t;
e[pre[i]^1].c+=t;
}
For(i,l,r) h[i]+=dis[i];
}
return {flow,cost};
}
int jiliguala=(init(),1);
#undef Tpl
#undef Tpw
}using Flow::init;using Flow::add;using Flow::Primal_Dual;
参考资料
我说 OI Wiki 神了有没有懂的 孩子们网络流写了 3 篇了，每次参考资料都有 OI Wiki。那你这不就是把 OI Wiki 的资料搬过来了吗
update on 7/14
更新了 Primal-Dual 原始对偶算法的内容。
浙公网安备 33010602011771号