费用流(自用，勿看)

注意：这是一篇个人学习笔记，如果有人因为某些原因点了进来并且要看一下，请一定谨慎地阅读，因为可能存在各种奇怪的错误，如果有人发现错误请指出谢谢！

最小费用最大流

https://www.luogu.org/problemnew/show/P3381

注意：以下算法不能处理费用带负环的图（如果我没搞错的话，应该只要原图（不考虑退流用的反向边）没有负环就行）；下面zkw博客里有关于负环的处理方法，先咕了

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 }

View Code

Primal-Dual 原始对偶算法（费用流）

zkw原文

先看这个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(d[0])*(n+1));
 31     memset(inq,0,sizeof(inq[0])*(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(d[0])*(n+1));
 58     memset(vis,0,sizeof(vis[0])*(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(vis[0])*(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(h[0])*(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 }

View Code

以上版本改了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(d[0])*(n+1));
 36     memset(inq,0,sizeof(inq[0])*(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(d[0])*(n+1));
 63     memset(vis,0,sizeof(vis[0])*(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(vis[0])*(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(h[0])*(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 }

View Code

单路增广版本：

  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(d[0])*(n+1));
 37     memset(inq,0,sizeof(inq[0])*(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(d[0])*(n+1));
 65     memset(vis,0,sizeof(vis[0])*(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(h[0])*(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 }

View Code