重修 MCMF 最小费用最大流
Primal-Dual 原始对偶算法
Upd 2023年10月24日
我们尝试将 \(f\) 次 spfa 改成 dijkstra 来优化复杂度。
但是边权不一定是非负的，考虑搬用全源最短路（Johnson）的做法（这里能找到），给每个点设一个势能 \(h_i\)（\(h_S=0\)）使得每条边加上势能差后都非负。
但是网络流的图是会变的，会有反向边出现，我们只需要每次增广之后将所有点势能 \(h_i\) 加上 \(dis_i\)，其中 \(dis_i\) 为新边权跑出来的最短路。
证明略，可看 OI-wiki
由于初始求势能要一次 spfa，复杂度 \(O(nm+mf\log)\)。
稠密图可以用 \(O(n^2)\) dijkstra，复杂度 \(O(nm+n^2f)\)，具体应用有 二分图最大权匹配。
用途
求在满足最大流的情况下最小的费用。
算法思路
每次找到一条费用之和最小的增广路径增广（反向边的费用是正向边的相反数）。
易证明贪心正确，或参考 OI Wiki 的证明。
具体实现
每次在残量网络中跑 Bellman-Ford 或 SPFA（按费用跑最短路），取路径上流量最小值增广，修改答案并改动网络上的这条路径。
注意我们不像 Dinic 一样按树形增广，而是像 EK 一样单条路径。树形也行，不过时间一样，何必呢。
时间复杂度
\(O(nmf)\) 点数、边数、最大流的积。
最坏可以到达 \(O(n^32^{n/2})\)，其实不刻意卡的话是伪多项式复杂度的。
由于 SPFA 已经玄学了，MCMF 时间复杂度更是玄学。
例题 & Code
P3381 【模板】最小费用最大流
#include<bits/stdc++.h>
using namespace std;
#define IOS ios::sync_with_stdio(0);cin.tie(0);cout.tie(0)
#define For(i,j,k) for(int i=j;i<=k;i++)
#define Rof(i,j,k) for(int i=j;i>=k;i--)
#define N 5010
#define M 50010
#define inf 0x7fffffff
struct edge{int to,flow,cost,nxt;}e[2*M];//2 times of edges!!!
int head[N],tot=1;
int n,m,S,T;
inline void adde(int x,int y,int flow,int cost){
e[++tot]=(edge){y,flow, cost,head[x]}; head[x]=tot;
e[++tot]=(edge){x,0 ,-cost,head[y]}; head[y]=tot;
}
int dis[N];//min toll to i (which SPFA works)
int flow[N];//the flow of the path with the min toll from S to i (just ONE path)
int pre[N],lst[N];//the last node and edge of the path with the min toll
queue<int> q;//for SPFA
bool ins[N];//in the queue or not to reduce the useless 'q.push()'
bool spfa(){
while(!q.empty()) q.pop();
For(i,1,n) ins[i]=pre[i]=0;
For(i,1,n) dis[i]=inf;
dis[S]=0;
ins[S]=1;
flow[S]=inf;
q.push(S);
int x;
while(!q.empty()){
x=q.front();
q.pop();
ins[x]=0;
for(int i=head[x];i;i=e[i].nxt){
int to=e[i].to;
assert(e[i].flow>=0);
if(e[i].flow==0 || dis[to]<=dis[x]+e[i].cost) continue;
dis[to]=dis[x]+e[i].cost;
pre[to]=x;
lst[to]=i;
flow[to]=min(flow[x],e[i].flow);
if(!ins[to]){
ins[to]=1;
q.push(to);
}
}
}
return dis[T]<inf;
}
void mcmf(){
int mf=0;
int mc=0;
while(spfa()){//find ONE augmentation path at a time to augment
mf+=flow[T];
mc+=flow[T]*dis[T];
int x=T;
while(x!=S){
e[lst[x] ].flow-=flow[T];
e[lst[x]^1].flow+=flow[T];
x=pre[x];
}
}
cout<<mf<<" "<<mc<<endl;
}
signed main(){IOS;
cin>>n>>m>>S>>T;
int x,y,w,c;
while(m--){
cin>>x>>y>>w>>c;
adde(x,y,w,c);
}
mcmf();
return 0;}
去注释版：
#include<bits/stdc++.h>
using namespace std;
#define IOS ios::sync_with_stdio(0);cin.tie(0);cout.tie(0)
#define For(i,j,k) for(int i=j;i<=k;i++)
#define Rof(i,j,k) for(int i=j;i>=k;i--)
#define N 5010
#define M 50010
#define inf 0x7fffffff
struct edge{int to,flow,cost,nxt;}e[2*M];
int head[N],tot=1;
int n,m,S,T;
inline void adde(int x,int y,int flow,int cost){
e[++tot]=(edge){y,flow, cost,head[x]}; head[x]=tot;
e[++tot]=(edge){x,0 ,-cost,head[y]}; head[y]=tot;
}
int dis[N];
int flow[N];
int pre[N],lst[N];
queue<int> q;
bool ins[N];
bool spfa(){
while(!q.empty()) q.pop();
For(i,1,n) ins[i]=pre[i]=0;
For(i,1,n) dis[i]=inf;
dis[S]=0;
ins[S]=1;
flow[S]=inf;
q.push(S);
int x;
while(!q.empty()){
x=q.front();
q.pop();
ins[x]=0;
for(int i=head[x];i;i=e[i].nxt){
int to=e[i].to;
assert(e[i].flow>=0);
if(e[i].flow==0 || dis[to]<=dis[x]+e[i].cost) continue;
dis[to]=dis[x]+e[i].cost;
pre[to]=x;
lst[to]=i;
flow[to]=min(flow[x],e[i].flow);
if(!ins[to]){
ins[to]=1;
q.push(to);
}
}
}
return dis[T]<inf;
}
void mcmf(){
int mf=0;
int mc=0;
while(spfa()){
mf+=flow[T];
mc+=flow[T]*dis[T];
int x=T;
while(x!=S){
e[lst[x] ].flow-=flow[T];
e[lst[x]^1].flow+=flow[T];
x=pre[x];
}
}
cout<<mf<<" "<<mc<<endl;
}
signed main(){IOS;
cin>>n>>m>>S>>T;
int x,y,w,c;
while(m--){
cin>>x>>y>>w>>c;
adde(x,y,w,c);
}
mcmf();
return 0;}
AcWing382. K取方格数
算是 传纸条 的加强版吧。
#include<bits/stdc++.h>
using namespace std;
#define IOS ios::sync_with_stdio(0);cin.tie(0);cout.tie(0)
#define fir first
#define sec second
#define mkp make_pair
#define pb emplace_back
#define mem(x,y) memset(x,y,sizeof(x))
#define For(i,j,k) for(int i=j;i<=k;i++)
#define Rof(i,j,k) for(int i=j;i>=k;i--)
#define ckmx(a,b) a=max(a,b)
#define ckmn(a,b) a=min(a,b)
#define fin(s) freopen(s,"r",stdin)
#define fout(s) freopen(s,"w",stdout)
#define file(s) fin(s".in");fout(s".out")
#define cerr cerr<<'_'
#define debug cerr<<"Passed line #"<<__LINE__<<endl
template<typename T>T ov(T x){cerr<<"Value: "<<x<<endl;return x;}
#define N 52
#define inf 0x3f3f3f3f
//Since the algorithm template is the minimum cost flow,
//and the topic requires the maximum value,
//the edge weights take the opposite number.
struct edge{int to,flow,cost,nxt;}e[8*N*N];
int head[2*N*N],tot=1;
int n,k,S,T,lim;
inline int num(int x,int y,bool z){return (((x-1)*n+y)<<1)+z;}
inline void adde(int x,int y,int flow,int cost){
e[++tot]=(edge){y,flow, cost,head[x]}; head[x]=tot;
e[++tot]=(edge){x,0 ,-cost,head[y]}; head[y]=tot;
}
int dis[2*N*N];//min toll to i (which SPFA works)
int flow[2*N*N];//the flow of the path with the min toll from S to i (just ONE path)
int pre[2*N*N],lst[2*N*N];//the last node and edge of the path with the min toll
queue<int> q;//for SPFA
bool ins[2*N*N];//in the queue or not to reduce the useless 'q.push()'
bool spfa(){
while(!q.empty()) q.pop();
For(i,1,lim) ins[i]=pre[i]=0;
For(i,1,lim) dis[i]=inf;
dis[S]=0;
ins[S]=1;
flow[S]=inf;
q.push(S);
int x;
while(!q.empty()){
x=q.front();
q.pop();
ins[x]=0;
for(int i=head[x];i;i=e[i].nxt){
int to=e[i].to;
assert(e[i].flow>=0);
if(e[i].flow==0 || dis[to]<=dis[x]+e[i].cost) continue;
dis[to]=dis[x]+e[i].cost;
pre[to]=x;
lst[to]=i;
flow[to]=min(flow[x],e[i].flow);
if(!ins[to]){
ins[to]=1;
q.push(to);
}
}
}
return dis[T]<inf;
}
int mcmf(){
int mf=0;//useless in this problem
int mc=0;
while(spfa()){//find ONE augmentation path at a time to augment
mf+=flow[T];
mc+=flow[T]*dis[T];
int x=T;
while(x!=S){
e[lst[x] ].flow-=flow[T];
e[lst[x]^1].flow+=flow[T];
x=pre[x];
}
}
return mc;
}
signed main(){IOS;
cin>>n>>k;
int x;
For(i,1,n) For(j,1,n){
cin>>x;
adde(num(i,j,0),num(i,j,1),1,-x);
adde(num(i,j,0),num(i,j,1),k,0);//'k' means inf
}
For(i,1,n) For(j,1,n-1) adde(num(i,j,1),num(i,j+1,0),k,0);
For(i,1,n-1) For(j,1,n) adde(num(i,j,1),num(i+1,j,0),k,0);
S=1,lim=T=num(n,n,1)+1;
adde(S,num(1,1,0),k,0);//but only here not,JUST k
adde(num(n,n,1),T,k,0);
cout<<(-mcmf())<<endl;//DONT forget to *(-1)
return 0;}
作者：ShaoJia，欢迎分享本文，转载时敬请注明原文来源链接。
浙公网安备 33010602011771号