Hierarchical Index

00001 //#include<scil/constraints/disjunctive.h>
00002 
00003 class dir_cut : public cons_obj {
00004   private:
00005   var_map<edge>& VMy;
00006   var_map<edge>& VMz;
00007   graph& G;
00008   edge_array<double>& D;
00009   list<row> NZ;
00010 
00011   public:
00012 
00013   dir_cut(graph& G_, node_array<double>& C, var_map<edge>& VMy_, 
00014           var_map<edge>& VMz_, edge_array<double>& D_, GraphWin& GW) 
00015   : cons_obj(Greater, 1), VMy(VMy_), VMz(VMz_), G(G_), D(D_) {
00016     node_array<double> CD(G);
00017     node u;
00018     forall_nodes(u, G) CD[u]=C[u];
00019     edge e;
00020     forall_edges(e, G) {
00021       if((CD[source(e)]>0) && (CD[target(e)]==0) && (CD[source(e)]>D[e]))
00022         CD[source(e)]=D[e];
00023       if((CD[target(e)]>0) && (CD[source(e)]==0) && (CD[target(e)]>D[e])) 
00024         CD[target(e)]=D[e];
00025     }
00026     forall_edges(e, G) {
00027       //GW.set_color(e, black);
00028       if(coeff_y(e, CD)>0) { 
00029         NZ.append(VMy[e]); 
00030         //GW.set_color(e, red);
00031       }
00032       if(coeff_z(e, CD)>0) { 
00033         NZ.append(VMz[e]); 
00034         //GW.set_color(e, red); 
00035       }
00036     };
00037     //GW.wait();
00038   };
00039 
00040   double coeff_y(edge e, node_array<double>& CD) {
00041     if((CD[source(e)]>0) && (CD[source(e)]<=D[e])) return 1;
00042     return 0;
00043   };
00044 
00045   double coeff_z(edge e, node_array<double>& CD) {
00046     if((CD[target(e)]>0) && (CD[target(e)]<=D[e])) return 1;
00047     return 0;
00048   };
00049 
00050   virtual void non_zero_entries(row& r) {
00051     row v;
00052     forall(v, NZ) r+=v;
00053   };
00054 };
00055 
00056 class Broadcast : public sym_constraint {
00057 
00058   private:
00059   graph& G;
00060   graph GD;
00061   node_array<node> GtoGD;
00062   edge_array<edge> REV;
00063   var_map<edge>& VM;
00064   var_map<edge>& VMy;
00065   var_map<edge>& VMz;
00066   var_map<edge> VMD;
00067   GraphWin& GW;
00068   edge_array<double>& C;
00069 
00070   public:
00071 
00072   Broadcast(graph& G_, var_map<edge>& VM_, var_map<edge>& VMy_, 
00073     var_map<edge>& VMz_, GraphWin& GW_, edge_array<double>& C_) 
00074     : G(G_), VM(VM_), VMy(VMy_), VMz(VMz_), GW(GW_), C(C_) {
00075   };
00076 
00077   void init(subproblem& s) {
00078     node_array<double> CD(G, 0);
00079 
00080     GtoGD.init(G);
00081     node u;
00082     forall_nodes(u, G)
00083       GtoGD[u]=GD.new_node();
00084 
00085     edge e,f,g;
00086     edge_array<double> C1(GD, 2*G.number_of_edges(), 0);
00087     REV.init(GD, 2*G.number_of_edges(), e);
00088     forall_edges(e, G) {
00089       f=GD.new_edge(GtoGD[source(e)], GtoGD[target(e)]);
00090       VMD[f]=VMy[e];
00091       C1[f]=-C[e];
00092       g=GD.new_edge(GtoGD[target(e)], GtoGD[source(e)]);
00093       VMD[g]=VMz[e];
00094       C1[g]=-C[e];
00095       REV[g]=f;
00096       REV[f]=g;
00097     };
00098 
00099     GD.sort_edges(C1);
00100     /*
00101     forall_edges(e, G) {
00102       CD[source(e)]=MAX_DOUBLE;
00103       CD[target(e)]=MAX_DOUBLE;
00104       s.add_basic_constraint(new dir_cut(G, CD, VMy, VMz, C, GW));
00105       CD[source(e)]=0;
00106       CD[target(e)]=0;
00107     };
00108     */
00109   };
00110 
00111   void mc_dfs(node u, edge_array<int>& c, edge_array<int>& f, node_array<bool>& r) {
00112     r[u]=true;
00113     edge e;
00114     forall_out_edges(e, u) 
00115       if((!r[target(e)]) && (f[e]!=c[e]))
00116         mc_dfs(target(e), c, f, r);
00117   };
00118 
00119   status separate(subproblem& s) {
00120     edge_array<int> cap(GD);
00121     edge e;
00122 
00123     status st=no_constraint_found;
00124     
00125     node_array<double> SE(GD, 0);
00126 
00127 
00128     forall_edges(e, GD) {
00129       cap[e]=(int) (1000*(s.value(VMD[e])+SE[source(e)]));
00130       SE[source(e)]+=s.value(VMD[e]);
00131     };
00132 
00133     edge_array<int> capr(GD);
00134     forall_edges(e, GD) capr[e]=cap[REV[e]];
00135 
00136     node u, v, w;
00137     int k;
00138     edge_array<int> f(GD);
00139     u=GD.choose_node();
00140     forall_nodes(v, GD) if (u!=v) {
00141       {
00142         k=MAX_FLOW(GD, u, v, cap, f);
00143         // cout<<k<<" flow\n";
00144         node_array<bool> r(GD, false);
00145         mc_dfs(u, cap, f, r);
00146         node_array<double> r1(G);
00147         forall_nodes(w, G) if(r[GtoGD[w]]) r1[w]=MAX_DOUBLE; else r1[w]=0;
00148         cons_obj* c=new dir_cut(G, r1, VMy, VMz, C, GW);
00149         if(c->violated(s)) {
00150           //cout<<k<<" violated\n";
00151           s.add_basic_constraint(c);
00152           st=constraint_found;
00153         };
00154       }
00155       {
00156         k=MAX_FLOW(GD, u, v, capr, f);
00157         // cout<<k<<" flow\n";
00158         node_array<bool> r(GD, false);
00159         mc_dfs(u, capr, f, r);
00160         node_array<double> r1(G);
00161         forall_nodes(w, G) if(!r[GtoGD[w]]) r1[w]=MAX_DOUBLE; else r1[w]=0;
00162         cons_obj* c=new dir_cut(G, r1, VMy, VMz, C, GW);
00163         if(c->violated(s)) {
00164           //cout<<k<<" violated\n";
00165           s.add_basic_constraint(c);
00166           st=constraint_found;
00167         };
00168       }
00169     }
00170     
00171     return st;
00172   }
00173   /*
00174   status close_subproblem(subproblem& s) {
00175     edge e;
00176     char ss[20];
00177     GW.get_window().set_line_width(4);
00178     forall_edges(e, G) {
00179       if(s.value(VMy[e])>0.1) {
00180         sprintf(ss, "%f", s.value(VMy[e]));
00181         GW.set_label(e, string(ss));
00182       }
00183       if(s.value(VMz[e])>0.1) {
00184         sprintf(ss, "%f", s.value(VMz[e]));
00185         GW.set_label(e, string(ss));
00186       }
00187     }
00188     GW.redraw();
00189     forall_edges(e, G) {
00190       if(s.value(VMy[e])>0.1)
00191         GW.get_window().draw_arrow(GW.get_position(source(e)), 
00192                                   GW.get_position(target(e)), green);
00193       if(s.value(VMz[e])>0.1)
00194         GW.get_window().draw_arrow(GW.get_position(target(e)), 
00195                                   GW.get_position(source(e)), green);
00196       if(s.value(VMy[e])>0.99) 
00197         GW.get_window().draw_arrow(GW.get_position(source(e)), 
00198                                   GW.get_position(target(e)), red);
00199       if(s.value(VMz[e])>0.99)
00200         GW.get_window().draw_arrow(GW.get_position(target(e)), 
00201                                   GW.get_position(source(e)), red);
00202     }
00203     //    GW.redraw();
00204     GW.wait();
00205     forall_edges(e, G) GW.set_color(e, black);
00206     GW.redraw();
00207     return continue_work;
00208   }
00209   */
00210 };
00211 
00212 double ComputeBroadcast(graph& G, edge_array<double>& Cost, list<edge>& T, 
00213                       GraphWin& GW) {
00214 
00215   list<edge> MST=MIN_SPANNING_TREE(G, Cost);
00216   node_array<double> MSTC(G,0);
00217   edge e;
00218   forall(e,MST) {
00219     MSTC[source(e)]=leda_max(MSTC[source(e)], Cost[e]);
00220     MSTC[target(e)]=leda_max(MSTC[target(e)], Cost[e]);
00221   };
00222   node u; double mstc=0;
00223   forall_nodes(u, G) mstc+=MSTC[u]*MSTC[u];
00224 
00225   //cout<<mstc<<" cost of mst\n";
00226 
00227   ILP_Problem IP(Optsense_Min);
00228 
00229   edge f; e=nil;
00230   var_map<edge> VM;
00231   var_map<edge> VMy;
00232   var_map<edge> VMz;
00233 
00234   forall_edges(e,G) { 
00235     VM[e]=IP.add_variable(-0.1, 0, 1, Vartype_Integer);
00236     VMy[e]=IP.add_variable(Cost[e]*Cost[e], 0, 1, Vartype_Integer);    
00237     VMz[e]=IP.add_variable(Cost[e]*Cost[e], 0, 1, Vartype_Integer);
00238   }
00239 
00240   forall_edges(e,G) {
00241     row r1;
00242     forall_out_edges(f,source(e)) if(Cost[f]>=Cost[e]-0.001) r1+=VMy[f];
00243     forall_in_edges(f, source(e)) if(Cost[f]>=Cost[e]-0.001) r1+=VMz[f];
00244     IP.add_basic_constraint(r1 >= VM[e]);
00245 
00246     row r2;
00247     forall_out_edges(f,target(e)) if(Cost[f]>=Cost[e]-0.001) r2+=VMy[f];
00248     forall_in_edges(f, target(e)) if(Cost[f]>=Cost[e]-0.001) r2+=VMz[f];
00249     IP.add_basic_constraint(r2 >= VM[e]);
00250   }    
00251 
00252   forall_nodes(u, G) {
00253     row r;
00254     forall_out_edges(e, u) r+=VMy[e];
00255     forall_in_edges(e, u) r+=VMz[e];
00256     IP.add_basic_constraint(r==1);
00257   };
00258 
00259   IP.add_sym_constraint(new SpanTree(G,VM));
00260 
00261   //IP.add_sym_constraint(new Broadcast(G, VM, VMy, VMz, GW, Cost));
00262 
00263   //IP.add_sym_constraint(new disjunctive());
00264 
00265   IP.optimize();
00266   
00267   double optc=0;
00268   forall_edges(e, G) {
00269     if(IP.get_solution(VM[e])>0.5) { 
00270       T.append(e);
00271     };
00272   }
00273   forall_nodes(u, G) MSTC[u]=0;
00274   forall(e,T) {
00275     MSTC[source(e)]=leda_max(MSTC[source(e)], Cost[e]);
00276     MSTC[target(e)]=leda_max(MSTC[target(e)], Cost[e]);
00277   };
00278   forall_nodes(u, G) optc+=MSTC[u]*MSTC[u];
00279 
00280   cout<<"Number of nodes "<<G.number_of_nodes()<<endl;
00281   cout<<"Cost of mst "<<mstc<<"\n";
00282   cout<<"Cost of optimal solution "<<optc<<"\n";
00283 
00284   cout<<"Save "<<(mstc-optc)/mstc*100<<"\n";
00285 
00286 
00287   return (mstc-optc)/mstc*100;
00288 } // END_COMPUTE_BROADCAST
sym_connect.c
