LocalImprove算法

  1 #include <bits/stdc++.h>
  2 #define maxn 100000
  3 #define maxm 3000000
  4 #define eps 1e-6
  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> A;             //A集合中的点
 24 set<int> Bs;            //Bs集合中的点
 25 vector<int> G[maxn];    //原始图上的边
 26 
 27 int n;                  //图上的点数
 28 int m;                  //图上的边数
 29 int k;                  //A集合中点的数量
 30 int vol_a;              //A集合的容量vol(A)
 31 double epsilon_sigma;   //必要的参数ε(σ)即ε,文中定义ε=1/(3*(1/σ-1))
 32 double sigma;           //必要的参数σ,可推导出σ=3ε/(3ε+1)
 33 
 34 /*图的初始化*/
 35 void init(){
 36     e = 0;
 37     memset(first,-1,sizeof(first));
 38 }
 39 
 40 /*加边*/
 41 void add_edge(int a,int b,double c){
 42     //printf("add_edge:from %d to %d,val = %.4f\n",a,b,c);
 43     edge[e].u = a;
 44     edge[e].v = b;
 45     edge[e].next = first[a];
 46     edge[e].w = c;
 47     first[a] = e++;
 48 
 49     edge[e].u = b;
 50     edge[e].v = a;
 51     edge[e].next = first[b];
 52     edge[e].w = 0;
 53     first[b] = e++;
 54 }
 55 
 56 /*bfs构造层次图*/
 57 int bfs(){
 58     int rear = 0;
 59     memset(d,-1,sizeof(d));
 60     d[S] = 0;q[rear++] = S;
 61     for(int i = 0;i < rear;i++){
 62         for(int j = first[q[i]];j != -1;j = edge[j].next){
 63             int to = edge[j].v;
 64             double val = edge[j].w;
 65             if(abs(val) > eps && d[to] == -1){
 66                 d[to] = d[q[i]]+1;
 67                 q[rear++] = to;
 68                 if(to == T) return 1;
 69             }
 70         }
 71     }
 72     return 0;
 73 }
 74 
 75 /*dfs计算阻塞流*/
 76 double dfs(int cur,double a){
 77     if(cur == T)    return a;
 78     for(int &i = work[cur];i != -1;i = edge[i].next){
 79         int to = edge[i].v;
 80         double val = edge[i].w;
 81         if(abs(val) > eps && d[to] == d[cur]+1){
 82             if(double t = dfs(to,min(a,val))){
 83                 edge[i].w -= t;edge[i^1].w += t;
 84                 return t;
 85             }
 86         }
 87     }
 88     return 0;
 89 }
 90 
 91 
 92 bool f1[maxn],f2[maxn];
 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 //min_cut输出最小割，并给出当前划分的导率conductance
113 double min_cut(){
114     /*for(int i = 0;i < e;i += 2){
115         printf("edge:from %d to %d,val = %.4f\n",edge[i].u,edge[i].v,edge[i].w);
116     }*/
117     memset(f1,0,sizeof(f1));
118     memset(f2,0,sizeof(f2));
119     dfs1(S);
120     dfs2(T);
121 
122     int cnt = 0;    //表示跨越集合的边数
123 
124     //printf("min_cut edge:\n");
125     for(int i = 0;i < e;i += 2){
126         if(f1[edge[i].u] && f2[edge[i].v] && edge[i].w < eps){
127             //printf("from %d to %d\n",edge[i].u,edge[i].v);
128             if(edge[i].u == S || edge[i].u == T)    continue;
129             if(edge[i].v == S || edge[i].v == T)    continue;
130             cnt++;
131         }
132     }
133     printf("\n");
134 
135     int vol1 = 0;
136     //set<int> S_new;//新的划分中S中的点
137     for(int i = S+1;i < T;i++){
138         if(f1[i]){
139             vol1 += deg[i];
140             //printf("flag1:%d\n",i);
141         }
142     }
143     /*set<int>::iterator iter;
144     for(iter = S_new.begin();iter != S_new.end();iter++){
145         int u = *iter;
146         for(int i = 0;i < G[u].size();i++){
147             int v = G[u][i];
148             if(S_new.find(v) != S_new.end())    continue;
149             cnt++;
150         }
151     }*/
152     printf("cnt = %d,vol1 = %d,total-vol1 = %d\n",cnt,vol1,2*m-vol1);
153     if(cnt == 0){
154         return 0;
155     }
156     double conductance = cnt*1.0/min(vol1,2*m-vol1);
157     printf("conductance = %.5f\n",conductance);
158     return conductance;
159 }
160 
161 double LocalFlow(double alpha){
162 
163     init();
164     Bs.clear();
165     double ans = 0;
166     printf("alpha=%.5f,vol_a=%d,epsilon_sigma=%.5f,sigma=%.5f\n",alpha,vol_a,epsilon_sigma,sigma);
167     int max_number_of_iterations = (int)(5.0/alpha*log(3*vol_a/sigma)/log(2));
168     printf("max_number_of_iterations=%d\n",max_number_of_iterations);
169     int number_of_iterations = 0;
170 
171     set<int>::iterator iter;
172 
173     //构造local graph G‘‘=G‘<Bs>
174     //edges(s,v) for all v ∈ A
175     for(iter = A.begin();iter != A.end();iter++){
176         int v = *iter;
177         add_edge(S,v,deg[v]);
178     }
179 
180     set<int> old_to_T_set;
181     set<int> old_A_union_Bs;
182 
183     while(++number_of_iterations < max_number_of_iterations){
184         set<int> A_union_Bs;
185         set_union(A.begin(),A.end(),Bs.begin(),Bs.end(),inserter(A_union_Bs, A_union_Bs.begin()));
186 
187         set<int> to_T_set;
188         set<int> N_A_union_Bs;
189 
190         for(iter = A_union_Bs.begin();iter != A_union_Bs.end();iter++){
191             int u = *iter;
192             for(int i = 0;i < G[u].size();i++){
193                 int v = G[u][i];
194                 if(A_union_Bs.find(v) != A_union_Bs.end())  continue;
195                 N_A_union_Bs.insert(v);
196             }
197         }
198         set_union(N_A_union_Bs.begin(),N_A_union_Bs.end(),Bs.begin(),Bs.end(),inserter(to_T_set, to_T_set.begin()));
199         //edges(v,t) for all v ∈ Bs union N(A union Bs)
200         for(iter = to_T_set.begin();iter != to_T_set.end();iter++){
201             int v = *iter;
202             if(old_to_T_set.find(v) != old_to_T_set.end())  continue;
203             add_edge(v,T,epsilon_sigma*deg[v]);
204         }
205         //edges(u,v) for all (v ∈ A union Bs) and (v ∈ V)
206         for(iter = A_union_Bs.begin();iter != A_union_Bs.end();iter++){
207             int u = *iter;
208             if(old_A_union_Bs.find(u) != old_A_union_Bs.end())  continue;
209             for(int i = 0;i < G[u].size();i++){
210                 int v = G[u][i];
211                 if(old_A_union_Bs.find(v) != old_A_union_Bs.end())  continue;
212                 add_edge(u,v,1/alpha);
213                 add_edge(v,u,1/alpha);
214             }
215         }
216 
217         /*
218         因为每次增广只需要在局部的图上进行，
219         因此我们维护两个集合old_to_T_set和old_A_union_Bs，
220         这样每次寻找阻塞流之前通过对比，
221         我们就可以把Bs集合更新后需要加入局部图的新边在图上建边即可。
222         */
223         old_to_T_set = to_T_set;
224         old_A_union_Bs = A_union_Bs;
225 
226         //bfs确定层次图
227         if(!bfs())  break;
228         //寻找阻塞流
229         memcpy(work,first,sizeof(first));
230         while(true){
231             double t = dfs(S,INF);
232             if(abs(t) < eps)    break;
233             ans += t;
234         }
235         //寻找A union Bs中指向超级汇点t中已经饱和的边，生成新的Bs集合
236         for(iter = N_A_union_Bs.begin();iter != N_A_union_Bs.end();iter++){
237             int u = *iter;
238             //printf("u = %d\n",u);
239             for(int i = first[u];i != -1;i = edge[i].next){
240                 int to = edge[i].v;
241                 double val = edge[i].w;
242                 if(to != T)  continue;
243                 //printf("u = %d,val = %.3f\n",u,val);
244                 if(val < eps){
245                     Bs.insert(u);
246                     break;
247                 }
248             }
249         }
250         //puts("----------------------");
251     }
252     /*printf("number_of_iterations=%d\n",number_of_iterations);
253     for(int i = 0;i < e;i += 2){
254         printf("edge:from %d to %d,val = %.4f\n",edge[i].u,edge[i].v,edge[i].w);
255     }*/
256     printf("LocalFlow(%.5f)=%.5f\n",alpha,ans);
257     return ans;
258 }
259 
260 
261 int main()
262 {
263     freopen("wing_nodal.txt","r",stdin);
264     //freopen("output.txt","w",stdout);
265 
266     scanf("%d%d%d%lf",&n,&m,&k,&epsilon_sigma);
267     sigma = 3*epsilon_sigma/(3*epsilon_sigma+1);
268     S = 0,T = n+1;
269     vol_a = 0;
270 
271     for(int i = 0;i < m;i++){
272         int a,b;
273         scanf("%d%d",&a,&b);
274         G[a].push_back(b);
275         G[b].push_back(a);
276         deg[a]++;deg[b]++;
277     }
278     for(int i = 0;i < k;i++){
279         int a;
280         scanf("%d",&a);
281         A.insert(a);
282         vol_a += deg[a];
283     }
284     set<int>::iterator iter;
285     int cut_a = 0;//统计E(A,V-A),即跨越A与V-A的边的数量
286     for(iter = A.begin();iter != A.end();iter++){
287         int u = *iter;
288         for(int i = 0;i < G[u].size();i++){
289             int v = G[u][i];
290             if(A.find(v) != A.end())    continue;
291             cut_a++;
292         }
293     }
294     double conductance = cut_a*1.0/min(vol_a,2*m-vol_a);
295     printf("init:cut_a = %d,vol_a = %d,vol_a/vol_v=%.5f\n",cut_a,vol_a,vol_a*1.0/(2*m));
296     printf("before algorithm,conductance = %.5f\n",conductance);
297 
298     clock_t time_begin,time_end;
299     time_begin = clock();
300     double alpha_max = 1.0,alpha_min = 0.0;
301     while(alpha_max - alpha_min > 1e-5){
302         double alpha_mid = (alpha_max+alpha_min)*0.5;
303         double max_flow = LocalFlow(alpha_mid);
304         conductance = min_cut();
305         if(abs(conductance) < eps){
306             alpha_min = alpha_mid;
307         }else{
308             alpha_max = alpha_mid;
309         }
310     }
311     puts("------------------result-----------------");
312     printf("best alpha=%.5f\n",alpha_max);
313     LocalFlow(alpha_max);
314     min_cut();
315     time_end = clock();
316     double delta = (double)(time_end - time_begin) / CLOCKS_PER_SEC;
317     printf("\n,delta_time = %.5f s\n",delta);
318     return 0;
319 }