Hierarchical Index

00001 #include<scil/scil.h>
00002 #include<scil/constraints/atour.h>
00003 #include<scil/constraints/path.h>
00004 //#include<scil/constraints/disjunctive.h>
00005 #include<LEDA/tuple.h>
00006 #include<LEDA/stream.h>
00007 #include<LEDA/string.h>
00008 #include<fstream.h>
00009 #include<LEDA/graph_alg.h>
00010 #include<LEDA/p_queue.h>
00011 
00012 #define MULT 2
00013 
00014 using namespace SCIL;
00015 using namespace LEDA;
00016 
00017 typedef two_tuple<node, node> prec;
00018 
00019 LEDA::string read_next_line(istream& is) {
00020   LEDA::string s="#";
00021   while((!is.eof()) && ((s.head(1)=="#") || (s.length()==0) || (s.head(1)==" "))) {
00022     s.read_line(is);
00023   };
00024   if(is.eof()) s="";
00025   return s;
00026 };
00027 
00028 class GATOUR : public sym_constraint {
00029   GRAPH<int,int>& G;
00030   var_map<edge>& VM;
00031   row_map<node>& NM;
00032   int n;
00033   list<node>* L;
00034   list<cons_obj*> CL;
00035   node sour,targ;
00036 
00037 public: 
00038   int addition(int c) {
00039     int d=MULT*10;
00040     if(c<d) return 10;
00041     d*=MULT*10;
00042     if(c<d) return 8;
00043     d*=MULT*10;
00044     if(c<d) return 6;
00045     d*=MULT*10;
00046     if(c<d) return 4;
00047     return 2;
00048   };
00049 
00050   GATOUR(GRAPH<int,int>& G_, int n_, var_map<edge>& VM_, row_map<node>& NM_) : G(G_), VM(VM_), NM(NM_) {
00051     n=n_;
00052 
00053     L=new list<node>[n];
00054 
00055     node v;
00056     forall_nodes(v,G) {
00057       if(G[v]==0) {
00058         sour=v;
00059       } else {
00060         if(G[v]==n*n) targ=v;
00061         else
00062           L[G[v]/n].append(v);
00063       }
00064     };
00065   };
00066 
00067   void MYDFS(graph& H, node s, node_array<bool>& R,edge_array<int>& E, edge_array<int>& F) {
00068     R[s]=true;
00069     edge e;
00070     forall_out_edges(e,s) if((F[e]!=E[e]) && (R[target(e)]==false))
00071       MYDFS(H,target(e),R,E,F);
00072     forall_in_edges(e,s) if((F[e]!=0) && (R[source(e)]==false))
00073       MYDFS(H,source(e),R,E,F);
00074   };
00075 
00076   status feasible(solution& s) {
00077     return feasible_solution;
00078 
00079     //return infeasible_solution;
00080   };
00081 
00082   void CDFS(subproblem& s, GRAPH<int,int>& G, node u, int* v, int nv, int n, int c, int& m, list<edge>& L1, list<edge>& L2) {
00083     edge e;
00084     forall_out_edges(e, u) {
00085       if((target(e)==targ) && (nv==n)) {
00086         if(c+G[e]<m) {
00087           m=c+G[e];
00088           L2.clear();
00089           edge f;
00090           forall(f,L1) L2.append(f);
00091           L2.append(e);
00092         }
00093       } else {
00094         if((s.value(VM[e])>0.001) && (v[G[target(e)]/n]==0)) {
00095           v[G[target(e)]/n]=1;
00096           L1.append(e);
00097           CDFS(s, G, target(e), v, nv+1, n, c+G[e], m, L1, L2);
00098           L1.pop_back();
00099           v[G[target(e)]/n]=0;
00100         }
00101       }
00102     };
00103   };
00104 
00105   status separate(subproblem& s) {
00106     int ccount=0;
00107     {
00108       bool found_solution;
00109       int v[n];
00110       edge e,e1,f;
00111       forall_out_edges(f,sour) if(s.value(VM[f])!=0) {
00112         found_solution=true;
00113         solution SOL;
00114         for(int i=0; i<n; i++) v[i]=0;
00115         int vi=0;
00116         int d=0;
00117         node w=sour;
00118         while(vi<n-1) {
00119           vi++;
00120           e1=nil;
00121           if(w!=sour) {
00122             forall_out_edges(e,w) if(target(e)!=targ) 
00123               if((v[G[target(e)]%n]==0) && ((e1==nil) || (1000*s.value(VM[e])-G[e]>1000*s.value(VM[e1])-G[e1]))) 
00124                 e1=e;
00125           } else e1=f;
00126           //assert(e1!=nil);
00127           if(e1==nil) { found_solution=false; break; }
00128           if(vi==1) v[G[target(e1)]/n]=1;
00129           v[G[target(e1)]%n]=1;
00130           w=target(e1);
00131           SOL.set_value(VM[e1],1);
00132           d+=G[e1];
00133         }
00134         forall_out_edges(e,w) if(target(e)==targ) {
00135           SOL.set_value(VM[e],1);
00136         }
00137         //cout<<d<<" sol\n";
00138         if(found_solution) s.save_solution(SOL);
00139       }
00140     }
00141 /*
00142     {
00143       int v[n];
00144       for(int i=0; i<n; i++) v[i]=0;
00145       int d=10000000;
00146       list<edge> L1, L2;
00147       CDFS(s,G,sour,v,0,n,0,d, L1, L2);
00148       cout<<d<<endl;
00149       if(d<10000000) {
00150         solution SOL;
00151         edge e;
00152         forall(e,L2) {
00153           cout<<G[e]<<" ";
00154           SOL.set_value(VM[e],1);
00155         };
00156         cout<<endl;
00157         s.save_solution(SOL);
00158       };
00159     };
00160 */
00161 
00162     return no_constraint_found;
00163 
00164     status sta=no_constraint_found;
00165 
00166     cons_obj* co;
00167     forall(co,CL) if(co->violated(s)) {
00168       sta=constraint_found;
00169       s.add_basic_constraint(co);
00170       ccount++;
00171     }
00172     if(sta==constraint_found) {
00173       cout<<"POOL SEPARATION\n";
00174       return sta;
00175     };
00176 
00177     node v,u;
00178     edge e,f;
00179     node_array<node> GtoH(G);
00180     graph H;
00181     edge_array<int> C(H,2*G.number_of_edges()+4*n*n,0);
00182     var_map<edge> HM;
00183     node HE[n];
00184     forall_nodes(v,G) GtoH[v]=H.new_node();
00185     for(int i=0; i<n; i++) HE[i]=H.new_node();
00186     forall_edges(e,G) if((s.value(VM[e])>0.001) && (target(e)!=targ)) {
00187       f=H.new_edge(GtoH[source(e)], GtoH[target(e)]);
00188       HM[f]=VM[e];
00189       C[f]=G[e];
00190       f=H.new_edge(GtoH[target(e)], GtoH[source(e)]);
00191       HM[f]=VM[e];
00192       C[f]=G[e];
00193     }
00194     forall_in_edges(e,targ) {
00195       f=H.new_edge(GtoH[source(e)],HE[G[source(e)]%n]);
00196       HM[f]=VM[e];
00197       C[f]=G[e];
00198       f=H.new_edge(HE[G[source(e)]%n],GtoH[targ]);
00199       HM[f]=VM[e];
00200       C[f]=G[e];
00201       f=H.new_edge(HE[G[source(e)]%n],GtoH[source(e)]);
00202       HM[f]=VM[e];
00203       C[f]=G[e];
00204       f=H.new_edge(GtoH[targ],HE[G[source(e)]%n]);
00205       HM[f]=VM[e];
00206       C[f]=G[e];
00207     }
00208 
00209     p_queue<int,cons_obj*> PQ;
00210 
00211     edge_array<int> E(H, H.number_of_edges()+2*n*n,0);
00212     edge_array<int> F(H, H.number_of_edges()+2*n*n,0);   
00213     list<edge> EL, FL;
00214     node_array<bool> R(H);
00215 
00216     forall_edges(e, H) {
00217       E[e]=leda_max((int) (10000.0*s.value(HM[e]))+addition(C[e]),0);
00218     }
00219 /*
00220     forall_nodes(u,G) if(s.value(NM[u])>0.6) {
00221       forall_nodes(v,G) if((v>u) && (s.value(NM[v])>0.6)) {
00222         cout<<"test\n";
00223         int d=MAX_FLOW(H,GtoH[u],GtoH[v],E,F);
00224         if(d<20000*(s.value(NM[u])+s.value(NM[v])-1)-100) {
00225           cout<<"cut\n";
00226           forall_nodes(v, H) R[v]=false;
00227           MYDFS(H, GtoH[u], R, E, F);
00228           row r;
00229           forall_edges(e,G) 
00230             if(R[GtoH[source(e)]]!=R[GtoH[target(e)]]) r+=VM[e];
00231           s.add_basic_constraint(r>=NM[u]+NM[v]);
00232           ccount++;
00233           sta=constraint_found;
00234         };
00235       };
00236     };
00237 */
00238 
00239     edge ts=H.new_edge(GtoH[targ],GtoH[sour]);
00240     E[ts]=20000;
00241     edge ts1=H.new_edge(GtoH[sour],GtoH[targ]);
00242     E[ts1]=20000;
00243     for(int i=0; i<n-1; i++) {
00244       forall(v, L[i]) {
00245         EL.append(H.new_edge(HE[i], GtoH[v]));
00246       }
00247       forall(e, EL) E[e]=20000;
00248 
00249       for(int j=i+1; j<n; j++) {
00250         forall(v, L[j]) {
00251           FL.append(H.new_edge(GtoH[v],HE[j]));
00252         }
00253         forall(e, FL) E[e]=20000;
00254 
00255         int d=MAX_FLOW(H,HE[i],HE[j],E,F);
00256 
00257         if(d<=18000) {
00258           //cout<<d<<" found constraint\n";
00259 
00260           forall_nodes(v, H) R[v]=false;
00261           MYDFS(H, HE[i], R, E, F);
00262           row r;
00263           forall_edges(e,G) 
00264             if((R[GtoH[source(e)]]) && (!R[GtoH[target(e)]])) {
00265               r+=VM[e];
00266             }
00267           s.add_basic_constraint(r>=1);
00268           ccount++;
00269           //if((20000-d)/r.size()>0) PQ.insert(-(20000-d)/r.size(),r>=2);
00270           //else cout<<d<<" "<<r.size()<<" denied\n";
00271           sta=constraint_found;
00272         };
00273 
00274         forall(e,FL) H.del_edge(e);
00275         FL.clear();
00276       }
00277       forall(e,EL) H.del_edge(e);
00278       EL.clear();
00279     };     
00280     if(sta==constraint_found) cout<<"YES\n"; else cout<<"NO\n";
00281 /*
00282     h_array<var, int> M;
00283     cons_obj* c;
00284     row_entry re;
00285     bool b;
00286     int i=0;
00287     while(!PQ.empty()) {
00288       c=PQ[PQ.find_min()];
00289       if(!c->violated(s)) cout<<"ERROR\n";
00290       //cout<<c->slack(s)<<endl;
00291       PQ.del_min();
00292       row r;
00293       c->non_zero_entries(r);
00294       //cout<<r.size()<<endl;
00295       b=true;
00296       //forall(re,r) if(M[re.Var]>=3) b=false;
00297       if(i<1500) {
00298         //forall(re,r) M[re.Var]++;
00299         s.add_basic_constraint(c);
00300         //CL.append(new cons_obj(*c));
00301         count++;
00302       } else delete c;
00303     }
00304 */
00305     cout<<ccount<<" constraints separated\n";
00306 
00307     return sta;
00308   };
00309 };
00310 
00311 
00312 int main(int argc, char** argv) {
00313 
00314   // READ_GRAPH
00315 
00316   ifstream is(argv[1]);
00317 
00318   GRAPH<int,int> G;
00319 
00320   LEDA::string_istream si(read_next_line(is));
00321   int n,i,j,k,l;
00322 
00323   si>>n;
00324   node N[n][n];
00325   LEDA::string Names[n];
00326 
00327   cout<<n<<" nodes in the graph\n";
00328 
00329   for(int i=0; i<n; i++) Names[i]=read_next_line(is);
00330 
00331   node s=G.new_node(0);
00332   node t=G.new_node(n*n);
00333 
00334   for(int i=0; i<n; i++)
00335     for(int j=0; j<n; j++)
00336       if(i!=j) { 
00337         N[i][j]=G.new_node(n*i+j);
00338         G.new_edge(N[i][j],t,0);
00339       }
00340 
00341   while(!is.eof()) {
00342     string_istream si(read_next_line(is));
00343     si>>i;
00344     if(!si.eof()) {
00345       si>>j; si>>k; 
00346       leda_string ls;
00347       char ch;
00348       ls.read_line(si);
00349       l=0;
00350       bool neg=false;
00351       for(int in=0; in<ls.length(); in++) {
00352         ch=ls[in];
00353         if(ch=='-') neg=true;
00354         if(ch==',') l*=MULT;
00355         if(ch=='0') l=10*l  ;
00356         if(ch=='1') l=10*l+1;
00357         if(ch=='2') l=10*l+2;
00358         if(ch=='3') l=10*l+3;
00359         if(ch=='4') l=10*l+4;
00360         if(ch=='5') l=10*l+5;
00361         if(ch=='6') l=10*l+6;
00362         if(ch=='7') l=10*l+7;
00363         if(ch=='8') l=10*l+8;
00364         if(ch=='9') l=10*l+9;
00365       }
00366       if(neg) l=-l;
00367       if(i!=-1) 
00368         G.new_edge(N[i][j],N[j][k],l);
00369       else
00370         G.new_edge(s,N[j][k],l);
00371     }
00372   };
00373 
00374   // END_READ_GRAPH
00375 
00376   row_map<node> EM;
00377 
00378   ILP_Problem IP(Optsense_Min);
00379 
00380   edge e,f;
00381   var_map<edge> VM;
00382   row_map<edge> VM1;
00383   forall_edges(e, G) {
00384     VM[e]=IP.add_variable(G[e],0,1, Vartype_Integer);
00385   };
00386 
00387   graph H;
00388   node HN[n+1];
00389   row sr[n];
00390   row tr[n];
00391   row_map<edge> VM2;
00392   for(int i=0; i<=n; i++) HN[i]=H.new_node();
00393   for(int i=0; i<n; i++) {
00394     for(int j=0; j<n; j++) if(i!=j) {
00395       row r;
00396       forall_in_edges(e, N[i][j]) {
00397         r+=VM[e];
00398         if(source(e)==s) sr[i]+=VM[e];
00399       }
00400       f=H.new_edge(HN[i],HN[j]);
00401       //VM1[f]=IP.add_variable(0,0,1, Vartype_Integer);
00402       VM2[f]=r;
00403       EM[N[i][j]]=VM2[f];
00404       //IP.add_basic_constraint(r==VM1[f]);
00405     }
00406   }
00407   forall_in_edges(e,t) { 
00408     tr[G[source(e)]%n]+=VM[e];
00409   }
00410   for(int i=0; i<n; i++) {
00411     f=H.new_edge(HN[i],HN[n]);
00412     //VM1[f]=IP.add_variable(0,0,1,Vartype_Integer);
00413     VM2[f]=tr[i];
00414     //IP.add_basic_constraint(tr[i]==VM1[f]);
00415     f=H.new_edge(HN[n],HN[i]);
00416     //VM1[f]=IP.add_variable(0,0,1,Vartype_Integer);
00417     VM2[f]=sr[i];
00418     //IP.add_basic_constraint(sr[i]==VM1[f]);
00419   };
00420 
00421   cout<<H.number_of_nodes()<<" nodes in H\n";
00422 
00423   IP.add_sym_constraint(new ATOUR(H, VM2));
00424 
00425   IP.add_sym_constraint(new GATOUR(G, n, VM, EM));
00426 
00427   node u,v;
00428   forall_nodes(u,G) {
00429     row r;
00430     if(u!=s) forall_in_edges(e,u) r-=VM[e];
00431     if(u!=t) forall_out_edges(e,u) r+=VM[e];
00432     if(u==s) IP.add_basic_constraint(r==1);
00433     else {
00434       if(u==t) IP.add_basic_constraint(r==-1);
00435       else IP.add_basic_constraint(r==0);
00436     };
00437   };
00438 /*
00439   for(int i=0; i<n; i++) {
00440     row r1;
00441     for(int j=0; j<n; j++) if(j!=i)
00442       forall_out_edges(e, N[i][j])
00443         r1+=VM[e];
00444     forall_in_edges(e,t) if(G[source(e)]%n==i) r1+=VM[e];
00445     IP.add_basic_constraint(r1==1);
00446   }
00447 */
00448   IP.optimize();
00449 
00450   int d=0;
00451   forall_edges(e, G) if(IP.get_solution(VM[e])>=0.5) {
00452     cout<<"("<<G[source(e)]/n<<","<<G[source(e)]%n<<")-("
00453              <<G[target(e)]/n<<","<<G[target(e)]%n<<"):"<<G[e]<<" edge\n";
00454     d+=G[e];
00455   } else G.del_edge(e);
00456 
00457   u=s;
00458   do {
00459     u=G.opposite(u,G.first_adj_edge(u));
00460     if(u!=t) cout<<G[u]/n<<" "<<Names[G[u]/n]<<endl;
00461   } while(u!=t);
00462   cout<<"\n\n";
00463 
00464   node_array<int> NN(H);
00465   for(int i=0; i<=n; i++) NN[HN[i]]=i;
00466   forall_edges(e,H) {
00467     double d=IP.get_solution(VM2[e]);
00468     //assert(d<=1.001);
00469     if(d>=0.5) cout<<d<<" "<<NN[source(e)]<<" "<<NN[target(e)]<<" edge\n";
00470   };
00471 
00472   cerr<<argv[1]<<" "<<d<<endl;
00473 
00474   exit(0);
00475 
00476   return (int) (d+0.5);
00477 };
00478
atsp2.c
