费用流(自用，勿看)
Posted hehe54321
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了费用流(自用，勿看)相关的知识，希望对你有一定的参考价值。
最小费用最大流
https://www.luogu.org/problemnew/show/P3381
spfa费用流
就是每一条边多一个属性：单位流量的费用
方法是在从0流开始找最大流过程中，保证每一次增广都沿着可行且费用最小的路增广，直到不能增广；最大流只有一个，因此最终一定能得到
就是把EK的bfs换成按费用为边权的spfa。。。
复杂度的话...EK的增广次数上限的分析好像还是适用的（应该吧?），所以是n^2*m^2
注意：
spfa不能在搜到T（u==T）时break
某一条正向边有cost的费用，对应反向边要加上-cost的费用
如果要卡常什么的可以加spfa的两个优化（不过好像复杂度会假掉？然而一般的图上都是快的）
如果改成longlong，注意33行改成sizeof(ll)
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<vector> 5 #include<queue> 6 using namespace std; 7 #define fi first 8 #define se second 9 #define mp make_pair 10 #define pb push_back 11 typedef long long ll; 12 typedef unsigned long long ull; 13 typedef pair<int,int> pii; 14 namespace F 15 { 16 17 struct E 18 { 19 int to,nxt,d,from,cap,flow; 20 }e[101000]; 21 int f1[5010],ne=1; 22 int S,T,n; 23 bool inq[5010]; 24 int d[5010],pre[5010],minn[5010]; 25 int flow,cost; 26 const int INF=0x3f3f3f3f; 27 void solve() 28 { 29 int k,u; 30 flow=cost=0; 31 while(1) 32 { 33 memset(d,0x3f,sizeof(int)*(n+1)); 34 memset(inq,0,sizeof(bool)*(n+1)); 35 queue<int> q; 36 q.push(S);d[S]=0;pre[S]=0;inq[S]=1;minn[S]=INF; 37 while(!q.empty()) 38 { 39 u=q.front();q.pop(); 40 inq[u]=0; 41 //if(u==T) break; 42 for(k=f1[u];k;k=e[k].nxt) 43 if(e[k].cap>e[k].flow&&d[u]+e[k].d<d[e[k].to]) 44 { 45 d[e[k].to]=d[u]+e[k].d; 46 pre[e[k].to]=k; 47 minn[e[k].to]=min(minn[u],e[k].cap-e[k].flow); 48 if(!inq[e[k].to]) 49 { 50 inq[e[k].to]=1; 51 q.push(e[k].to); 52 } 53 } 54 } 55 if(d[T]==INF) break; 56 flow+=minn[T];cost+=d[T]*minn[T]; 57 for(k=pre[T];k;k=pre[e[k].from]) 58 { 59 e[k].flow+=minn[T];e[k^1].flow-=minn[T]; 60 } 61 } 62 } 63 void me(int a,int b,int c,int d) 64 { 65 e[++ne].to=b;e[ne].nxt=f1[a];f1[a]=ne;e[ne].from=a; 66 e[ne].cap=c;e[ne].d=d; 67 e[++ne].to=a;e[ne].nxt=f1[b];f1[b]=ne;e[ne].from=b; 68 e[ne].cap=0;e[ne].d=-d; 69 } 70 71 } 72 73 int n,m,S,T; 74 int main() 75 { 76 int i,a,b,c,d; 77 scanf("%d%d%d%d",&n,&m,&S,&T);F::S=S;F::T=T;F::n=n; 78 for(i=1;i<=m;i++) 79 { 80 scanf("%d%d%d%d",&a,&b,&c,&d); 81 F::me(a,b,c,d); 82 } 83 F::solve(); 84 printf("%d %d",F::flow,F::cost); 85 return 0; 86 }
Primal-Dual 原始对偶算法（费用流）
先看这个johnson算法介绍：https://blog.csdn.net/howarli/article/details/73824179
说白了，就是可以证明：给每个点一个权值h[i]，然后使得图中每一条边(a,b)的长度增加h[b]-h[a]，最终在新图中随便跑两点(u,v)间的最短路d‘，设原图中两点间最短路为d，那么d‘=d+h[u]-h[v]
只要构造一种方法，使得每一条边加上这个权值后非负，就能一次spfa后都只跑dijkstra了
算法的流程是：
spfa-->更新h[i]-->增广-->dijkstra-->更新h[i]-->增广-->dijkstra...直到某一次spfa/dijkstra发现没有增广路时跳出
具体方法是：
每一次最短路，从任一点sx开始向每个点跑（只考虑残量网络中边）（sx为S或T，具体见最后）
更新h[i]时，将每一个h[i]加上前一次最短路时的d[i]（表示起始点sx到i的最短路）
跑最短路时，查询(u,v,w)的边权时就改成w+h[u]-h[v]
复杂度是：n*m^2*log
原因：
其实是不难的。。却因为理不清思路还有搞复杂了，用了很多时间
目标就是：证明任意一次跑最短路（除了第一次）时，任意边(u,v,w)满足w+h[u]-h[v]>=0；h的取值只要满足这一个条件即可，并没有其他限制（当初就是因为总是去想h的其他意义，导致证不出来）
不太好考虑，从最开始考虑
第一遍spfa+更新：根据最短路的性质，更新完后显然h[u]+w>=h[v]，即w+h[u]-h[v]>=0
第一遍增广：设d[i]为第一遍spfa时sx到i的不带h影响的最短路；在增广中有一些边会变成其反边，设这条反边是(v,u,-w)，由于增广只会走最短路上的边，可知d[u]+w=d[v]，那么增广后-w+d[v]-d[u]=-(w+d[u]-d[v])=0，因此增广不会产生新负权边
第二遍最短路：由于此时如果给边权加上h的影响的话，图中并没有负权边，因此可以用dijkstra算法跑出最短路；设d[i]为这一遍做出来的sx到i的带h的影响的最短路；设d‘[i]为这一遍的sx到i的不带h的影响的最短路（当然并没有算出来，只是用于证明）；根据johnson算法的那个证明，d[i]=d‘[i]+h[sx]-h[i]=d‘[i]-h[i]
第二遍更新：每个h[i]加上一个d[i]，即加上d‘[i]-h[i]，因此每个h[i]变成d‘[i]，即当前图不带h的影响的sx到i的最短路
可以发现每一次增广之前，h[i]都能成为当前图（不带h的影响的）sx到i的最短路，又可以用来辅助增广路
可以发现这样的证明对于任意一轮的操作都是生效的（只要上一轮符合），可以当做一个归纳法的证明
证明参考：https://www.luogu.org/problemnew/solution/P3381
实现：
可以用多路增广，当然仍然只能增广在(S,T)最短路上的点；注意到如果不加更多的限制条件，且残量网络中有零费用环，那么会陷入无限递归（而且即使判掉这个，这个“多路增广”我感觉跟爆搜也没什么区别(?猜的，并不清楚)），因此可以干脆写成一次增广每个点最多经过一次（我没有更好的方法）；这样子一次可能增广不完全部路径，因此每次要增广的时候不断进行多路增广直到无法增广；（这样的多路增广感觉很不好啊，但是网上要么是单路增广，要么就是这样(?有别的也说不定)，也许真的只能这样？）
同一轮中，每一条增广路的总费用都相等（毕竟都是最短路），等于这一轮时的h[S]/h[T]（毕竟这等于当前图实际上S到T的最短路）；这样就可以计算每一次增广增加的费用了
（sx可以设为S；sx也可以设为T，当然这样的话要对反图跑最短路，听说会快一点，另外前面多路增广时的判断最短路要改一下）
实现参考：https://panda-2134.blog.luogu.org/solution-p3381
多路增广版本：
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<vector> 5 #include<queue> 6 using namespace std; 7 #define fi first 8 #define se second 9 #define mp make_pair 10 #define pb push_back 11 typedef long long ll; 12 typedef unsigned long long ull; 13 typedef pair<int,int> pii; 14 namespace F 15 { 16 17 struct E 18 { 19 int to,nxt,d,from,cap,flow; 20 }e[101000]; 21 int f1[5010],ne=1; 22 int S,T,n; 23 bool inq[5010],*vis=inq; 24 int d[5010]; 25 int h[5010]; 26 int flow,cost; 27 bool spfa() 28 { 29 int u,k,k1; 30 memset(d,0x3f,sizeof(int)*(n+1)); 31 memset(inq,0,sizeof(bool)*(n+1)); 32 queue<int> q; 33 q.push(T);d[T]=0;inq[T]=1; 34 while(!q.empty()) 35 { 36 u=q.front();q.pop(); 37 inq[u]=0; 38 for(k1=f1[u];k1;k1=e[k1].nxt) 39 { 40 k=k1^1; 41 if(e[k].cap>e[k].flow&&d[u]+e[k].d<d[e[k].from]) 42 { 43 d[e[k].from]=d[u]+e[k].d; 44 if(!inq[e[k].from]) 45 { 46 inq[e[k].from]=1; 47 q.push(e[k].from); 48 } 49 } 50 } 51 } 52 return d[S]!=0x3f3f3f3f; 53 } 54 bool dij() 55 { 56 int u,k,k1;pii t; 57 memset(d,0x3f,sizeof(int)*(n+1)); 58 memset(vis,0,sizeof(bool)*(n+1)); 59 priority_queue<pii,vector<pii>,greater<pii> > q; 60 q.push(mp(0,T));d[T]=0; 61 while(!q.empty()) 62 { 63 t=q.top();q.pop(); 64 u=t.se; 65 if(vis[u]) continue; 66 vis[u]=1; 67 for(k1=f1[u];k1;k1=e[k1].nxt) 68 { 69 k=k1^1; 70 if(e[k].cap>e[k].flow&&d[u]+e[k].d+h[u]-h[e[k].from]<d[e[k].from]) 71 { 72 d[e[k].from]=d[u]+e[k].d+h[u]-h[e[k].from]; 73 q.push(mp(d[e[k].from],e[k].from)); 74 } 75 } 76 } 77 return d[S]!=0x3f3f3f3f; 78 } 79 void update() 80 { 81 for(int i=1;i<=n;i++) h[i]+=d[i]; 82 } 83 int dfs(int u,int x) 84 { 85 if(u==T||x==0) return x; 86 int flow=0,f; 87 vis[u]=1; 88 for(int k=f1[u];k;k=e[k].nxt) 89 if(!vis[e[k].to]&&e[k].cap>e[k].flow&&h[e[k].to]==h[u]-e[k].d) 90 { 91 f=dfs(e[k].to,min(x-flow,e[k].cap-e[k].flow)); 92 e[k].flow+=f;e[k^1].flow-=f;flow+=f; 93 if(flow==x) return flow; 94 } 95 return flow; 96 } 97 void augment() 98 { 99 int f; 100 while(1) 101 { 102 memset(vis,0,sizeof(int)*(n+1)); 103 f=dfs(S,0x3f3f3f3f); 104 if(!f) break; 105 flow+=f;cost+=f*h[S]; 106 } 107 } 108 void solve() 109 { 110 flow=cost=0; 111 memset(h,0,sizeof(int)*(n+1)); 112 if(!spfa()) return; 113 do 114 { 115 update(); 116 augment(); 117 }while(dij()); 118 } 119 void me(int a,int b,int c,int d) 120 { 121 e[++ne].to=b;e[ne].nxt=f1[a];f1[a]=ne;e[ne].from=a; 122 e[ne].cap=c;e[ne].d=d; 123 e[++ne].to=a;e[ne].nxt=f1[b];f1[b]=ne;e[ne].from=b; 124 e[ne].cap=0;e[ne].d=-d; 125 } 126 127 } 128 129 int n,m,S,T; 130 int main() 131 { 132 int i,a,b,c,d; 133 scanf("%d%d%d%d",&n,&m,&S,&T);F::S=S;F::T=T;F::n=n; 134 for(i=1;i<=m;i++) 135 { 136 scanf("%d%d%d%d",&a,&b,&c,&d); 137 F::me(a,b,c,d); 138 } 139 F::solve(); 140 printf("%d %d",F::flow,F::cost); 141 return 0; 142 }
以上版本改了pbds的堆，跑的很快：
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<vector> 5 #include<queue> 6 #include<ext/pb_ds/assoc_container.hpp> 7 #include<ext/pb_ds/priority_queue.hpp> 8 using namespace std; 9 #define fi first 10 #define se second 11 #define mp make_pair 12 #define pb push_back 13 typedef long long ll; 14 typedef unsigned long long ull; 15 typedef pair<int,int> pii; 16 namespace F 17 { 18 19 struct E 20 { 21 int to,nxt,d,from,cap,flow; 22 }e[101000]; 23 int f1[5010],ne=1; 24 int S,T,n; 25 bool inq[5010],*vis=inq; 26 int d[5010]; 27 int h[5010]; 28 int flow,cost; 29 typedef __gnu_pbds::priority_queue<pii,greater<pii> > pq; 30 pq::point_iterator it[5010]; 31 //const int INF=0x3f3f3f3f; 32 bool spfa() 33 { 34 int u,k,k1; 35 memset(d,0x3f,sizeof(int)*(n+1)); 36 memset(inq,0,sizeof(bool)*(n+1)); 37 queue<int> q; 38 q.push(T);d[T]=0;inq[T]=1; 39 while(!q.empty()) 40 { 41 u=q.front();q.pop(); 42 inq[u]=0; 43 for(k1=f1[u];k1;k1=e[k1].nxt) 44 { 45 k=k1^1; 46 if(e[k].cap>e[k].flow&&d[u]+e[k].d<d[e[k].from]) 47 { 48 d[e[k].from]=d[u]+e[k].d; 49 if(!inq[e[k].from]) 50 { 51 inq[e[k].from]=1; 52 q.push(e[k].from); 53 } 54 } 55 } 56 } 57 return d[S]!=0x3f3f3f3f; 58 } 59 bool dij() 60 { 61 int i,u,k,k1;pii t; 62 memset(d,0x3f,sizeof(int)*(n+1)); 63 memset(vis,0,sizeof(bool)*(n+1)); 64 pq q; 65 for(i=1;i<=n;i++) it[i]=q.end(); 66 it[T]=q.push(mp(0,T));d[T]=0; 67 while(!q.empty()) 68 { 69 t=q.top();q.pop(); 70 u=t.se; 71 if(vis[u]) continue; 72 vis[u]=1; 73 for(k1=f1[u];k1;k1=e[k1].nxt) 74 { 75 k=k1^1; 76 if(e[k].cap>e[k].flow&&d[u]+e[k].d+h[u]-h[e[k].from]<d[e[k].from]) 77 { 78 d[e[k].from]=d[u]+e[k].d+h[u]-h[e[k].from]; 79 if(it[e[k].from]!=q.end()) q.modify(it[e[k].from],mp(d[e[k].from],e[k].from)); 80 else it[e[k].from]=q.push(mp(d[e[k].from],e[k].from)); 81 } 82 } 83 } 84 return d[S]!=0x3f3f3f3f; 85 } 86 void update() 87 { 88 for(int i=1;i<=n;i++) h[i]+=d[i]; 89 } 90 int dfs(int u,int x) 91 { 92 if(u==T||x==0) return x; 93 int flow=0,f; 94 vis[u]=1; 95 for(int k=f1[u];k;k=e[k].nxt) 96 if(!vis[e[k].to]&&e[k].cap>e[k].flow&&h[e[k].to]==h[u]-e[k].d) 97 { 98 f=dfs(e[k].to,min(x-flow,e[k].cap-e[k].flow)); 99 e[k].flow+=f;e[k^1].flow-=f;flow+=f; 100 if(flow==x) return flow; 101 } 102 return flow; 103 } 104 void augment() 105 { 106 int f; 107 while(1) 108 { 109 memset(vis,0,sizeof(int)*(n+1)); 110 f=dfs(S,0x3f3f3f3f); 111 if(!f) break; 112 flow+=f;cost+=f*h[S]; 113 } 114 } 115 void solve() 116 { 117 flow=cost=0; 118 memset(h,0,sizeof(int)*(n+1)); 119 if(!spfa()) return; 120 do 121 { 122 update(); 123 augment(); 124 }while(dij()); 125 } 126 void me(int a,int b,int c,int d) 127 { 128 e[++ne].to=b;e[ne].nxt=f1[a];f1[a]=ne;e[ne].from=a; 129 e[ne].cap=c;e[ne].d=d; 130 e[++ne].to=a;e[ne].nxt=f1[b];f1[b]=ne;e[ne].from=b; 131 e[ne].cap=0;e[ne].d=-d; 132 } 133 134 } 135 136 int n,m,S,T; 137 int main() 138 { 139 int i,a,b,c,d; 140 scanf("%d%d%d%d",&n,&m,&S,&T);F::S=S;F::T=T;F::n=n; 141 for(i=1;i<=m;i++) 142 { 143 scanf("%d%d%d%d",&a,&b,&c,&d); 144 F::me(a,b,c,d); 145 } 146 F::solve(); 147 printf("%d %d",F::flow,F::cost); 148 return 0; 149 }
单路增广版本：
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<vector> 5 #include<queue> 6 #include<ext/pb_ds/assoc_container.hpp> 7 #include<ext/pb_ds/priority_queue.hpp> 8 using namespace std; 9 #define fi first 10 #define se second 11 #define mp make_pair 12 #define pb push_back 13 typedef long long ll; 14 typedef unsigned long long ull; 15 typedef pair<int,int> pii; 16 namespace F 17 { 18 19 struct E 20 { 21 int to,nxt,d,from,cap,flow; 22 }e[101000]; 23 int f1[5010],ne=1; 24 int S,T,n; 25 bool inq[5010],*vis=inq; 26 int d[5010]; 27 int h[5010]; 28 int pre[5010]; 29 int flow,cost; 30 typedef __gnu_pbds::priority_queue<pii,greater<pii> > pq; 31 pq::point_iterator it[5010]; 32 //const int INF=0x3f3f3f3f; 33 bool spfa() 34 { 35 int u,k,k1; 36 memset(d,0x3f,sizeof(int)*(n+1)); 37 memset(inq,0,sizeof(bool)*(n+1)); 38 queue<int> q; 39 q.push(T);d[T]=0;inq[T]=1;pre[T]=0; 40 while(!q.empty()) 41 { 42 u=q.front();q.pop(); 43 inq[u]=0; 44 for(k1=f1[u];k1;k1=e[k1].nxt) 45 { 46 k=k1^1; 47 if(e[k].cap>e[k].flow&&d[u]+e[k].d<d[e[k].from]) 48 { 49 d[e[k].from]=d[u]+e[k].d; 50 pre[e[k].from]=k; 51 if(!inq[e[k].from]) 52 { 53 inq[e[k].from]=1; 54 q.push(e[k].from); 55 } 56 } 57 } 58 } 59 return d[S]!=0x3f3f3f3f; 60 } 61 bool dij() 62 { 63 int i,u,k,k1;pii t; 64 memset(d,0x3f,sizeof(int)*(n+1)); 65 memset(vis,0,sizeof(bool)*(n+1)); 66 pq q; 67 for(i=1;i<=n;i++) it[i]=q.end(); 68 it[T]=q.push(mp(0,T));d[T]=0;pre[T]=0; 69 while(!q.empty()) 70 { 71 t=q.top();q.pop(); 72 u=t.se; 73 if(vis[u]) continue; 74 vis[u]=1; 75 for(k1=f1[u];k1;k1=e[k1].nxt) 76 { 77 k=k1^1; 78 if(e[k].cap>e[k].flow&&d[u]+e[k].d+h[u]-h[e[k].from]<d[e[k].from]) 79 { 80 d[e[k].from]=d[u]+e[k].d+h[u]-h[e[k].from]; 81 pre[e[k].from]=k; 82 if(it[e[k].from]!=q.end()) q.modify(it[e[k].from],mp(d[e[k].from],e[k].from)); 83 else it[e[k].from]=q.push(mp(d[e[k].from],e[k].from)); 84 } 85 } 86 } 87 return d[S]!=0x3f3f3f3f; 88 } 89 void update() 90 { 91 for(int i=1;i<=n;i++) h[i]+=d[i]; 92 } 93 void augment() 94 { 95 int k,f=0x3f3f3f3f; 96 for(k=pre[S];k;k=pre[e[k].to]) 97 { 98 f=min(f,e[k].cap-e[k].flow); 99 } 100 for(k=pre[S];k;k=pre[e[k].to]) 101 { 102 e[k].flow+=f;e[k^1].flow-=f; 103 } 104 flow+=f;cost+=f*h[S]; 105 } 106 void solve() 107 { 108 flow=cost=0; 109 memset(h,0,sizeof(int)*(n+1)); 110 if(!spfa()) return; 111 do 112 { 113 update(); 114 augment(); 115 }while(dij()); 116 } 117 void me(int a,int b,int c,int d) 118 { 119 e[++ne].to=b;e[ne].nxt=f1[a];f1[a]=ne;e[ne].from=a; 120 e[ne].cap=c;e[ne].d=d; 121 e[++ne].to=a;e[ne].nxt=f1[b];f1[b]=ne;e[ne].from=b; 122 e[ne].cap=0;e[ne].d=-d; 123 } 124 125 } 126 127 int n,m,S,T; 128 int main() 129 { 130 int i,a,b,c,d; 131 scanf("%d%d%d%d",&n,&m,&S,&T);F::S=S;F::T=T;F::n=n; 132 for(i=1;i<=m;i++) 133 { 134 scanf("%d%d%d%d",&a,&b,&c,&d); 135 F::me(a,b,c,d); 136 } 137 F::solve(); 138 printf("%d %d",F::flow,F::cost); 139 return 0; 140 }
以上是关于费用流(自用，勿看)的主要内容，如果未能解决你的问题，请参考以下文章