Improve算法

  1 #include <bits/stdc++.h>
  2 #define maxn 100000
  3 #define maxm 2000000
  4 #define eps 1e-7
  5 #define INF 1<<30
  6 using namespace std;
  7 
  8 struct Edge{
  9     int u;          //u表示边的起点
 10     int v;          //v表示边的终点
 11     int next;       //next表示下一条边编号
 12     double w;       //w表示边的容量
 13 }edge[maxm];        //网络流建模时用到的边
 14 
 15 int e;                  //边的计数器
 16 int S;                  //超级源点
 17 int T;                  //超级汇点
 18 int first[maxn];        //图邻接表头结点
 19 int d[maxn];            //层次图距离标记
 20 int work[maxn];         //dinic优化
 21 int q[maxn];            //bfs队列
 22 int deg[maxn];          //每个点的度数
 23 set<int>::iterator iter;
 24 set<int> A;             //A集合中的点
 25 vector<int> G[maxn];    //原始图上的边
 26 
 27 bool f1[maxn];          //残量网络与S连通的点集
 28 bool f2[maxn];          //残量网络与T连通的点集
 29 
 30 int n;      //图上的点数
 31 int m;      //图上的边数
 32 int k;      //A集合中点的数量
 33 int vol_a;  //A集合的容量vol(A)
 34 int cut_a;  //跨越A和V-A的边数
 35 double f;       //f(A)=vol(A)/vol(V-A)
 36 double alpha;   //与网络流边容量相关的参数α
 37 /*图的初始化*/
 38 void init(){
 39     e = 0;
 40     memset(first,-1,sizeof(first));
 41 }
 42 
 43 /*加边*/
 44 void add_edge(int a,int b,double c){
 45     //printf("add_edge:from %d to %d,val = %.4f\n",a,b,c);
 46     edge[e].u = a;
 47     edge[e].v = b;
 48     edge[e].next = first[a];
 49     edge[e].w = c;
 50     first[a] = e++;
 51 
 52     edge[e].u = b;
 53     edge[e].v = a;
 54     edge[e].next = first[b];
 55     edge[e].w = 0;
 56     first[b] = e++;
 57 }
 58 
 59 /*bfs构造层次图*/
 60 int bfs(){
 61     int rear = 0;
 62     memset(d,-1,sizeof(d));
 63     d[S] = 0;q[rear++] = S;
 64     for(int i = 0;i < rear;i++){
 65         for(int j = first[q[i]];j != -1;j = edge[j].next){
 66             int to = edge[j].v;
 67             double val = edge[j].w;
 68             if(abs(val) > eps && d[to] == -1){
 69                 d[to] = d[q[i]]+1;
 70                 q[rear++] = to;
 71                 if(to == T) return 1;
 72             }
 73         }
 74     }
 75     return 0;
 76 }
 77 
 78 /*dfs计算阻塞流*/
 79 double dfs(int cur,double a){
 80     if(cur == T)    return a;
 81     for(int &i = work[cur];i != -1;i = edge[i].next){
 82         int to = edge[i].v;
 83         double val = edge[i].w;
 84         if(abs(val) > eps && d[to] == d[cur]+1){
 85             if(double t = dfs(to,min(a,val))){
 86                 edge[i].w -= t;edge[i^1].w += t;
 87                 return t;
 88             }
 89         }
 90     }
 91     return 0;
 92 }
 93 
 94 void dfs1(int u){
 95     f1[u] = 1;
 96     for(int i = first[u];i != -1;i = edge[i].next){
 97         int to = edge[i].v;
 98         double val = edge[i].w;
 99         if(!f1[to] && val > eps)    dfs1(to);
100     }
101 }
102 
103 void dfs2(int u){
104     f2[u] = 1;
105     for(int i = first[u];i != -1;i = edge[i].next){
106         int to = edge[i].v;
107         double val = edge[i^1].w;
108         if(!f2[to] && val > eps)    dfs2(to);
109     }
110 }
111 
112 void build(double alpha){
113     init();
114     for(int i = 1;i <= n;i++){
115         if(A.find(i) != A.end()){
116             add_edge(S,i,deg[i]);
117         }else{
118             add_edge(i,T,f*deg[i]);
119         }
120     }
121     for(int u = 1;u <= n;u++){
122         for(int i = 0;i < G[u].size();i++){
123             int v = G[u][i];
124             add_edge(u,v,1.0/alpha);
125             add_edge(v,u,1.0/alpha);
126         }
127     }
128 }
129 
130 double dinic(){
131     double ans = 0,t;
132     while(bfs()){
133         memcpy(work,first,sizeof(first));
134         while(true){
135             t = dfs(S,INF);
136             if(t < eps) break;
137             ans += t;
138         }
139     }
140     return ans;
141 }
142 
143 double min_cut(){
144     memset(f1,0,sizeof(f1));
145     memset(f2,0,sizeof(f2));
146     dfs1(S);
147     dfs2(T);
148 
149     int cnt = 0;    //表示跨越集合的边数
150 
151     //printf("min_cut edge:\n");
152     for(int i = 0;i < e;i += 2){
153         if(f1[edge[i].u] && f2[edge[i].v] && edge[i].w < eps){
154             //printf("from %d to %d\n",edge[i].u,edge[i].v);
155             if(edge[i].u == S || edge[i].u == T)    continue;
156             if(edge[i].v == S || edge[i].v == T)    continue;
157             cnt++;
158         }
159     }
160 
161     int vol1 = 0,vol2 = 0;
162     for(int i = S+1;i < T;i++){
163         if(f1[i]){
164             vol1 += deg[i];
165             //printf("flag1:%d\n",i);
166         }
167         if(f2[i]){
168             if(A.find(i) != A.end())  vol2 += deg[i];
169         }
170     }
171 
172     printf("cnt = %d,vol1 = %d,total-vol1 = %d\n",cnt,vol1,2*m-vol1);
173     if(cnt == 0)   return 0;
174     double conductance = cnt*1.0/min(vol1,2*m-vol1);
175     printf("conductance = %.5f\n",conductance);
176 
177     double QAS = cnt*1.0/(vol1-f*vol2);
178     return QAS;
179 }
180 
181 int main()
182 {
183     freopen("wing_nodal.txt","r",stdin);
184     freopen("output.txt","w",stdout);
185     scanf("%d%d%d",&n,&m,&k);
186     S = 0,T = n+1;
187     vol_a = cut_a = 0;
188     for(int i = 0;i < m;i++){
189         int a,b;
190         scanf("%d%d",&a,&b);
191         G[a].push_back(b);
192         G[b].push_back(a);
193         deg[a]++;deg[b]++;
194     }
195     for(int i = 0;i < k;i++){
196         int a;
197         scanf("%d",&a);
198         A.insert(a);
199         vol_a += deg[a];
200     }
201     for(iter = A.begin();iter != A.end();iter++){
202         int u = *iter;
203         for(int i = 0;i < G[u].size();i++){
204             int v = G[u][i];
205             if(A.find(v) != A.end())    continue;
206             cut_a++;
207         }
208     }
209     f = vol_a*1.0/(2*m-vol_a);
210     alpha = cut_a*1.0/vol_a;
211     double best_alpha = alpha;
212     printf("init:cut_a = %d,vol_a = %d\n",cut_a,vol_a);
213     printf("before algorithm,conductance = %.5f\n",alpha);
214     clock_t time_begin,time_end;
215     time_begin = clock();
216     while(true){
217         build(alpha);
218         double maxflow = dinic();
219         double QAS = min_cut();
220         printf("maxflow = %.5f\n",maxflow);
221         printf("QAS = %.5f\n",QAS);
222         if(abs(QAS) < eps){
223             break;
224         }else{
225             best_alpha = alpha;
226         }
227         if(QAS < alpha){
228             alpha = QAS;
229         }else{
230             break;
231         }
232     }
233     puts("------------------result-----------------");
234     printf("best alpha = %.5f\n",best_alpha);
235     build(best_alpha);
236     dinic();
237     min_cut();
238 
239     time_end = clock();
240     double delta = (double)(time_end - time_begin) / CLOCKS_PER_SEC;
241     printf("\n,delta_time = %.5f s\n",delta);
242     return 0;
243 }