Hierarchical Index

00001 
00002 /*
00003 void MAX_DIJKSTRA(graph& G, node s, edge_array<int>& C, node_array<int>& D, node_array<edge>& P) {
00004   node_pq<int> PQ(G);
00005   P[s]=nil;
00006   PQ.insert(s, -MAX_INT);
00007   node w;
00008   D[w]=MAX_INT;
00009   edge e;
00010   while(!PQ.empty()) {
00011     w=PQ.del_min();
00012     forall_out_edges(e, w) if(C[e]!=0) {
00013       if(D[target(e)]<leda_min(D[w], C[e])) {
00014         if(D[target(e)]==0) PQ.insert(target(e), -leda_min(D[w], C[e]));
00015         else PQ.decrease_p(target(e), -leda_min(D[w], C[e]));
00016         D[target(e)]=leda_min(D[w], C[e]);
00017       }
00018     }
00019   }
00020 }
00021 
00022 void MIN_ST_CUT_SPECIAL_SUB(graph& G, node& s, node& t,
00023                             edge_array<int>& C, edge_array<int>& F,
00024                             list<node>& CUT) {
00025   node_array<int> RCN(G, 0);
00026   edge_array<int> RC(G, 0);
00027   edge e;
00028   forall_rev_edges(e, G) {
00029     RCN[source(e)]+=C[e]-F[e];
00030     RC[e]=RCN[source(e)];
00031   };
00032   
00033   node_array<int> D(G, 0);
00034   node_array<edge> Pred(G);
00035   MAX_DIJKSTRA(G, s, RC, D, Pred);
00036 
00037   if(D[t]!=0) {
00038     node w=t;
00039     while(w!=s) {
00040       F[Pred[w]]+=D[t];
00041       w=source(Pred[w]);
00042     }
00043     MIN_ST_CUT_SPECIAL_SUB(G, s, t, C, F, CUT);
00044   } else {
00045     node u;
00046     forall_nodes(u, G) if (D[u]!=0) { CUT.append(u); }
00047   };
00048 };
00049 */
00050 void mdfs(graph& G, node u, edge_array<double>& U, edge_array<double>& U1, 
00051           node_array<bool>& R, list<node>& C) {
00052   C.append(u);
00053   R[u]=true;
00054   edge e;
00055   forall_out_edges(e, u) if((U[e]-U1[e]>0.5) && (!R[target(e)]))
00056     mdfs(G, target(e), U, U1, R, C);
00057   //forall_in_edges(e, u) if((U1[e]>0.5) && (!R[source(e)]))
00058   //  mdfs(G, source(e), U, U1, R, C);
00059 };
00060 
00061 void MIN_ST_CUT_SPECIAL(graph& G, edge_array<double>& D, node s, node t,
00062                         edge_array<int>& C, list<node>& CUT) {
00063 /*
00064   G.sort_edges(D);
00065   edge e;
00066   forall_edges(e, G) cout<<D[e]<<" ";
00067   cout<<"\nIncreasing?\n";
00068   edge_array<int> F(G, 0);
00069   MIN_ST_CUT_SPECIAL_SUB(G, s, t, C, F, CUT);
00070 */
00071 
00072   ILP_Problem ILP(Optsense_Max);
00073   edge e,f;
00074   double d;
00075   var_map<edge> VMM;
00076   forall_edges(e, G) {
00077     if(source(e)==s) d=1; else d=0;
00078     VMM[e]=ILP.add_variable(d, 0, 1000, Vartype_Float);
00079   };
00080   forall_edges(e, G) {
00081     row r; d=0;
00082     forall_out_edges(f, source(e)) if(D[f]>=D[e]) {
00083       r+=VMM[f]; d+=leda_max(C[f],0);
00084     }
00085     ILP.add_basic_constraint(r<=d);
00086   }
00087 
00088   node v;
00089   forall_nodes(v, G) if((v!=s) && (v!=t)) {
00090     row r;
00091     forall_out_edges(e, v) r+=VMM[e];
00092     forall_in_edges(e, v) r-=VMM[e];
00093     ILP.add_basic_constraint(r==0);
00094   }
00095 
00096   ILP.optimize();
00097 
00098   edge_array<double> U(G);
00099   edge_array<double> U1(G);
00100   forall_edges(e, G) {
00101     row r; d=0;
00102     forall_out_edges(f, source(e)) if(D[f]>=D[e]) {
00103       U1[e]+=ILP.get_solution(VMM[e]); U[e]+=C[e];
00104     }
00105   }
00106 
00107   node_array<bool> R(G, false);
00108   mdfs(G, s, U, U1, R, CUT);
00109 };
00110 
00111 
00112 class dir_cut2 : public cons_obj {
00113   private:
00114   var_map<edge>& VMy;
00115   graph& G;
00116   edge_array<double>& D;
00117   list<row> NZ;
00118 
00119   public:
00120 
00121   dir_cut2(graph& G_, node_array<double>& C, var_map<edge>& VMy_, 
00122            edge_array<double>& D_, GraphWin& GW) 
00123     : cons_obj(Greater, 1), VMy(VMy_), G(G_), D(D_) {
00124     node_array<double> CD(G);
00125     node u;
00126     forall_nodes(u, G) CD[u]=C[u];
00127     edge e;
00128     forall_edges(e, G) {
00129       if((CD[source(e)]>0) && (CD[target(e)]==0) && (CD[source(e)]>D[e]))
00130         CD[source(e)]=D[e];
00131     }
00132     forall_edges(e, G) {
00133       //GW.set_color(e, black);
00134       if(coeff(e, CD)>0) { 
00135         NZ.append(VMy[e]); 
00136         //GW.set_color(e, red);
00137       }
00138     };
00139     //GW.wait();
00140   };
00141 
00142   double coeff(edge e, node_array<double>& CD) {
00143     if((CD[source(e)]>0) && (CD[source(e)]<=D[e])) return 1;
00144     return 0;
00145   };
00146 
00147   virtual void non_zero_entries(row& r) {
00148     row v;
00149     forall(v, NZ) r+=v;
00150   };
00151 };
00152 
00153 
00154 class Broadcast2 : public sym_constraint {
00155 
00156   private:
00157   graph& G;
00158   var_map<edge>& VM;
00159   var_map<edge>& VMy;
00160   GraphWin& GW;
00161   edge_array<double>& C;
00162 
00163   public:
00164 
00165   Broadcast2(graph& G_, var_map<edge>& VM_, var_map<edge>& VMy_, 
00166              GraphWin& GW_, edge_array<double>& C_) 
00167     : G(G_), VM(VM_), VMy(VMy_), GW(GW_), C(C_) {
00168   };
00169 
00170   void init(subproblem& s) {
00171   };
00172 
00173   void mc_dfs(node u, edge_array<int>& c, edge_array<int>& f, node_array<bool>& r) {
00174     r[u]=true;
00175     edge e;
00176     forall_out_edges(e, u) 
00177       if((!r[target(e)]) && (f[e]!=c[e]))
00178         mc_dfs(target(e), c, f, r);
00179   };
00180 
00181   status separate(subproblem& s) {
00182     edge_array<int> cap(G);
00183     edge e;
00184 
00185     status st=no_constraint_found;
00186 
00187     node_array<double> SE(G, 0);
00188    
00189     forall_edges(e, G) {
00190       cap[e]=(int) (1000*(s.value(VMy[e])+SE[source(e)]));
00191       SE[source(e)]+=s.value(VMy[e]);
00192     };
00193 
00194     node u, v, w;
00195     int k;
00196     edge_array<int> f(G);
00197     u=G.choose_node();
00198     forall_nodes(v, G) if (u!=v) {
00199       {
00200         k=MAX_FLOW(G, u, v, cap, f);
00201         // cout<<k<<" flow\n";
00202         node_array<bool> r(G, false);
00203         mc_dfs(u, cap, f, r);
00204         node_array<double> r1(G);
00205         forall_nodes(w, G) if(!r[w]) r1[w]=MAX_DOUBLE; else r1[w]=0;
00206         cons_obj* c=new dir_cut2(G, r1, VMy, C, GW);
00207         if(c->violated(s)) {
00208           //cout<<k<<" violated\n";
00209           s.add_basic_constraint(c);
00210           st=constraint_found;
00211         };
00212       }
00213       {
00214         k=MAX_FLOW(G, v, u, cap, f);
00215         // cout<<k<<" flow\n";
00216         node_array<bool> r(G, false);
00217         mc_dfs(v, cap, f, r);
00218         node_array<double> r1(G);
00219         forall_nodes(w, G) if(!r[w]) r1[w]=MAX_DOUBLE; else r1[w]=0;
00220         cons_obj* c=new dir_cut2(G, r1, VMy, C, GW);
00221         if(c->violated(s)) {
00222           //cout<<k<<" violated\n";
00223           s.add_basic_constraint(c);
00224           st=constraint_found;
00225         };
00226       }
00227     }
00228     /*
00229     if(st==constraint_found) return st;
00230 
00231     forall_edges(e, G) {
00232       cap[e]=(int) (1000*(s.value(VMy[e])));
00233     };
00234 
00235     u=G.choose_node();
00236     forall_nodes(v, G) if (u!=v) {
00237       {
00238         list<node> L;
00239         MIN_ST_CUT_SPECIAL(G, C, u, v, cap, L);
00240         node_array<double> r1(G, 0);
00241         if(L.size()==G.number_of_nodes()) cout<<"Fehler\n"; else {
00242           forall(w, L) r1[w]=MAX_DOUBLE;
00243           cons_obj* c=new dir_cut2(G, r1, VMy, C, GW);
00244           if(c->violated(s)) {
00245             //cout<<k<<" violated\n";
00246             s.add_basic_constraint(c);
00247             st=constraint_found;
00248           };
00249         }
00250       }
00251       {
00252         list<node> L;
00253         MIN_ST_CUT_SPECIAL(G, C, v, u, cap, L);
00254         node_array<double> r1(G, 0);
00255         if(L.size()==G.number_of_nodes()) cout<<"Fehler2\n"; else {
00256           forall(w, L) r1[w]=MAX_DOUBLE;
00257           cons_obj* c=new dir_cut2(G, r1, VMy, C, GW);
00258           if(c->violated(s)) {
00259             //cout<<k<<" violated\n";
00260             s.add_basic_constraint(c);
00261             st=constraint_found;
00262           }
00263         };
00264       }
00265     }
00266     */
00267     return st;
00268   }
00269   /*
00270   status close_subproblem(subproblem& s) {
00271     edge e;
00272     char ss[20];
00273     GW.get_window().set_line_width(4);
00274     forall_edges(e, G) {
00275       if(s.value(VMy[e])>0.1) {
00276         sprintf(ss, "%f", s.value(VMy[e]));
00277         GW.set_label(e, string(ss));
00278       }
00279       if(s.value(VMz[e])>0.1) {
00280         sprintf(ss, "%f", s.value(VMz[e]));
00281         GW.set_label(e, string(ss));
00282       }
00283     }
00284     GW.redraw();
00285     forall_edges(e, G) {
00286       if(s.value(VMy[e])>0.1)
00287         GW.get_window().draw_arrow(GW.get_position(source(e)), 
00288                                   GW.get_position(target(e)), green);
00289       if(s.value(VMz[e])>0.1)
00290         GW.get_window().draw_arrow(GW.get_position(target(e)), 
00291                                   GW.get_position(source(e)), green);
00292       if(s.value(VMy[e])>0.99) 
00293         GW.get_window().draw_arrow(GW.get_position(source(e)), 
00294                                   GW.get_position(target(e)), red);
00295       if(s.value(VMz[e])>0.99)
00296         GW.get_window().draw_arrow(GW.get_position(target(e)), 
00297                                   GW.get_position(source(e)), red);
00298     }
00299     //    GW.redraw();
00300     GW.wait();
00301     forall_edges(e, G) GW.set_color(e, black);
00302     GW.redraw();
00303     return continue_work;
00304   }
00305   */
00306 };
00307 
00308 double ComputeBroadcast2(graph& G, edge_array<double>& Cost, list<edge>& T, 
00309                       GraphWin& GW) {
00310 
00311   list<edge> MST=MIN_SPANNING_TREE(G, Cost);
00312   node_array<double> MSTC(G,0);
00313   edge e;
00314   forall(e,MST) {
00315     MSTC[source(e)]=leda_max(MSTC[source(e)], Cost[e]);
00316     MSTC[target(e)]=leda_max(MSTC[target(e)], Cost[e]);
00317   };
00318   node u; double mstc=0;
00319   forall_nodes(u, G) mstc+=MSTC[u]*MSTC[u];
00320 
00321   graph H;
00322   edge f;
00323   node_array<node> GtoH(G);
00324   node_map<point> P(H);
00325   edge_array<edge> HtoG(H,2*G.number_of_edges(),e);
00326 
00327   forall_nodes(u, G) { 
00328     GtoH[u]=H.new_node();
00329     P[GtoH[u]]=GW.get_position(u);
00330   };
00331 
00332   forall_edges(e, G) {
00333     f=H.new_edge(GtoH[source(e)], GtoH[target(e)]);
00334     HtoG[f]=e;
00335     f=H.new_edge(GtoH[target(e)], GtoH[source(e)]);
00336     HtoG[f]=e;
00337   };
00338 
00339   edge_array<double> Cost2(H);
00340   forall_edges(e, H) Cost2[e]=Cost[HtoG[e]];
00341 
00342   cout<<"created copy\n";
00343 
00344   ILP_Problem IP(Optsense_Min);
00345 
00346   e=nil;
00347   var_map<edge> VM;
00348   var_map<edge> VMy;
00349 
00350   forall_edges(e,H) {
00351     VM[e]=IP.add_variable(0, 0, 1, Vartype_Integer);
00352     VMy[e]=IP.add_variable(Cost[HtoG[e]]*Cost[HtoG[e]], 0, 1, Vartype_Integer);    
00353   }
00354 
00355   cout<<"made constr 1\n";
00356 
00357   forall_edges(e,H) {
00358     row r1;
00359     forall_out_edges(f,source(e)) if(Cost[HtoG[f]]>=Cost[HtoG[e]]-0.001) r1+=VMy[f];
00360     IP.add_basic_constraint(r1 >= VM[e]);
00361   }
00362 
00363   cout<<"made cosntr 2\n";
00364 
00365   forall_nodes(u, H) {
00366     row r;
00367     forall_out_edges(e, u) r+=VMy[e];
00368     IP.add_basic_constraint(r==1);
00369   };
00370 
00371   cout<<"mad constr 3\n";
00372 
00373   IP.add_sym_constraint(new StronglyConnected(H,VM));
00374 
00375   IP.add_sym_constraint(new Broadcast2(H, VM, VMy, GW, Cost2));
00376 
00377   IP.optimize();
00378 
00379   /*  
00380   forall_edges(e, H) {
00381     if(IP.get_solution(VM[e])>0.5) {
00382       GW.get_window().draw_arrow(P[source(e)],P[target(e)],red);
00383       T.append(HtoG[e]);
00384     };
00385   }
00386   GW.wait();
00387   */
00388 
00389   list<edge> T1;
00390   double optc=0;
00391   forall_edges(e, H) {
00392     if(IP.get_solution(VM[e])>0.5) { 
00393       T1.append(e);
00394       T.append(HtoG[e]);
00395     };
00396   }
00397   node_array<double> OPTC(H, 0);
00398   forall(e,T1) {
00399     OPTC[source(e)]=leda_max(OPTC[source(e)], Cost2[e]);
00400   };
00401   forall_nodes(u, H) optc+=OPTC[u]*OPTC[u];
00402 
00403   cout<<"Number of nodes "<<G.number_of_nodes()<<endl;
00404   cout<<"Cost of mst "<<mstc<<"\n";
00405   cout<<"Cost of optimal solution "<<optc<<"\n";
00406 
00407   cout<<"Save "<<(mstc-optc)/mstc*100<<"\n";
00408 
00409 
00410   return (mstc-optc)/mstc*100;
00411 } // END_COMPUTE_BROADCAST2
asym_connect.c
