費用流(自用，勿看)

本文轉載自查看原文 2018-09-12 22:12 847 圖論--網絡流--費用流

注意：這是一篇個人學習筆記，如果有人因為某些原因點了進來並且要看一下，請一定謹慎地閱讀，因為可能存在各種奇怪的錯誤，如果有人發現錯誤請指出謝謝！

最小費用最大流

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

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 網絡流學習--費用流貪心相關/模擬網絡流、費用流細節梳理/模板（貪心，模擬費用流，棧）【LOJ#3097】[SNOI2019]通信（費用流）網絡流--最小費用最大流（理解）【BZOJ4842】[Neerc2016]Delight for a Cat 線性規划+費用流 [bzoj3876][AHOI2014]支線劇情——上下界費用流 MySql的基操勿六已刪除，勿加！！！ Word 腳本（自用）京東推廣鏈接（自用）