SCL - MCMF based on Primal Dual

2017-03-19 16:35:01

最小费用流  之  原始对偶 (Primal-Dual) 算法

这个算法在搞竞赛时并没有学... 但是在申研时和通通问我的题目中都涉及到了，决定学习下、写一发。

学习推荐博客：https://artofproblemsolving.com/community/c1368h1020435

 

  1 #include <stdio.h>
  2 #include <string.h>
  3 #include <deque>
  4 #include <vector>
  5 #include <algorithm>
  6 using namespace std;
  7 
  8 const int MAXN = 1010;
  9 const int MAXM = 2010;
 10 const int INF = 0x3F3F3F3F;
 11 
 12 struct edge{
 13     int v,nxt,cost,cp;
 14     edge() {}
 15     edge(int a,int b,int c,int d) : v(a),nxt(b),cost(c),cp(d) {}
 16 };
 17 
 18 struct MCMF_Primal_Dual{
 19     edge e[MAXM];
 20     int N;
 21     int sou,sin,cost,piS;
 22     int first[MAXN],ecnt;
 23     int link[MAXN];
 24     bool vis[MAXN];
 25 
 26     void init(int n,int a,int b){
 27         N = n;
 28         sou = a;
 29         sin = b;
 30         memset(first,-1,sizeof(first));
 31         ecnt = ecnt = cost = piS = 0;
 32     }
 33     void add_edge(int u,int v,int fee,int cap){
 34         e[ecnt] = edge(v,first[u],fee,cap);
 35         first[u] = ecnt++;
 36         e[ecnt] = edge(u,first[v],-fee,0);
 37         first[v] = ecnt++;
 38     }
 39     int aug(int p,int minf){
 40         if(p == sin) return cost += piS * minf,minf;
 41         vis[p] = true;
 42         int tmp = minf;
 43         for(int i = first[p]; ~i; i = e[i].nxt){
 44             if(e[i].cp && !e[i].cost && !vis[e[i].v]){
 45                 int d = aug(e[i].v,tmp < e[i].cp ? tmp : e[i].cp);
 46                 e[i].cp -= d,e[i ^ 1].cp += d,tmp -= d;
 47                 if(!tmp) return minf;
 48             }
 49         }
 50         return minf - tmp;
 51     }
 52     bool modlabel(){
 53         static int D[MAXN];
 54         memset(D,0x3f,sizeof(D));
 55         D[sin] = 0;
 56         static deque<int> Q;
 57         Q.push_back(sin);
 58         while(Q.size()){
 59             int dt,now = Q.front(); Q.pop_front();
 60             for(int i = first[now]; ~i; i = e[i].nxt){
 61                 if(e[i ^ 1].cp && (dt = D[now] - e[i].cost) < D[e[i].v])
 62                     (D[e[i].v] = dt) <= D[Q.size() ? Q.front() : 0] ? 
 63                         Q.push_front(e[i].v) : Q.push_back(e[i].v);
 64             }
 65         }
 66         for(int i = 0; i <= N; ++i)
 67             for(int j = first[i]; ~j; j = e[j].nxt)
 68                 e[j].cost += D[e[j].v] - D[i];
 69         piS += D[sou];
 70         return D[sou] < INF;
 71     }
 72     pair<int,int> solve(){
 73         int sum_flow = 0;
 74         while(modlabel()){
 75             while(1){
 76                 memset(vis,false,sizeof(vis));
 77                 int f = aug(sou,INF);
 78                 if(!f) break;
 79                 sum_flow += f;
 80             }
 81         }
 82         return make_pair(cost,sum_flow);
 83     }
 84     int Dfs(int p,int minf){
 85         if(p == sin || minf == 0) return minf;
 86         vis[p] = true;
 87         for(int i = first[p]; ~i; i = e[i].nxt){
 88             if(e[i ^ 1].cp && !vis[e[i].v]){
 89                 int d = Dfs(e[i].v,min(minf,e[i ^ 1].cp));
 90                 if(d > 0){
 91                     e[i ^ 1].cp -= d;
 92                     link[p] = e[i].v;
 93                     // e[i].cp += d; 这句可加可不加，如果需要从残余网络复原，则加
 94                     return d;
 95                 }
 96             }
 97         }
 98         return 0;
 99     }
100     vector<vector<int> > find_path(){
101         vector<vector<int> > res;
102         while(1){
103             memset(vis,false,sizeof(vis));
104             int flow = Dfs(sou,INF);
105             if(!flow) break;
106             vector<int> tmp(0);
107             int p = sou;
108             while(p != sin){
109                 tmp.push_back(p);
110                 p = link[p];
111             }
112             tmp.push_back(sin);
113             tmp.push_back(flow);
114             res.push_back(tmp);
115         }
116         return res;
117     }
118 }SSR;
119 
120 int main(){
121     int n = 5;
122     SSR.init(n,0,n);
123     SSR.add_edge(0,1,0,2);
124     SSR.add_edge(1,2,1,1); SSR.add_edge(2,5,1,1);
125     SSR.add_edge(1,3,1,1); SSR.add_edge(3,5,1,1);
126     SSR.add_edge(1,4,2,1); SSR.add_edge(4,5,2,1);
127     pair<int,int> res = SSR.solve();
128     vector<vector<int> > res_path = SSR.find_path();
129     for(int i = 0; i < res_path.size(); ++i){
130         vector<int> &now = res_path[i];
131         for(int j = 0; j < now.size(); ++j) printf("%d ",now[j]);
132         puts("");
133     }
134     printf("%d %d\n",res.first,res.second);
135     return 0;
136 }

 

posted @ 2017-03-19 16:39  Naturain  阅读(255)  评论(0)    收藏  举报