angleTSP.cc

#include<fstream>
#include<scil/scil.h>
#include<scil/constraints/tour.h>
#include<scil/core/monomial.h>
#include<boost/graph/adjacency_list.hpp>
#include<boost/graph/iteration_macros.hpp>

using namespace std;
using namespace SCIL;

   typedef boost::adjacency_list<vecS, vecS, undirectedS> Graph;
   typedef boost::graph_traits<Graph> GraphTraits;
   typedef GraphTraits::vertex_descriptor vertex_descriptor;
   typedef GraphTraits::edge_descriptor edge_descriptor;

void read_instance(char* file_name, Graph& G, std::vector<edge_descriptor>& VE, var_map<edge_descriptor>& VM, std::vector<monomial>& lm, ILP_Problem& IP) {

   ifstream file;
   file.open(file_name);

   int n;
   file >> n; //read the size of the problem

   vertex_descriptor v, w;

   int m = (n*(n-1))/2;

   std::pair<double,double> coord;
   std::vector<pair<double,double> > coordVector;
   for( int i = 0; i < n; i++ ) {
         add_vertex(G);
         file >> coord.first;
         file >> coord.second;
         coordVector.push_back(coord);
   }

   BGL_FORALL_VERTICES(v, G, Graph){
      BGL_FORALL_VERTICES(w, G, Graph){
         if(v!=w){
            if(!edge(v,w,G).second)
               add_edge(v,w,G);
         }
      }
   }

   edge_descriptor e;
   std::vector<var> VI;
   BGL_FORALL_EDGES(e, G, Graph) {
      VM[e] = IP.add_binary_variable(0.);
      VI.push_back(VM[e]);
      VE.push_back(e);
   }

   //compute angle between every pair of two connected edges
   std::pair<double, double> v1,v2,vt,vector1,vector2;
   monomial mon;
   for( int i = 0; i < m; i++ ) {
      for( int j = 0; j < m; j++ ) {
         if( i != j ) {
            if( source(VE[i], G) == source(VE[j], G)) {
               vt = coordVector[source(VE[i], G)];
               v1 = coordVector[target(VE[i], G)];
               v2 = coordVector[target(VE[j], G)];
               vector1 = make_pair(v1.first - vt.first, v1.second - vt.second);
               vector2 = make_pair(vt.first - v2.first, vt.second - v2.second);
            } else if( source(VE[i], G) == target(VE[j], G) ) {
               vt = coordVector[source(VE[i], G)];
               v1 = coordVector[target(VE[i], G)];
               v2 = coordVector[source(VE[j], G)];
               vector1 = make_pair(v1.first - vt.first, v1.second - vt.second);
               vector2 = make_pair(vt.first - v2.first, vt.second - v2.second);
            } else if( target(VE[i], G) == source(VE[j], G) ) {
               vt = coordVector[source(VE[j], G)];
               v1 = coordVector[target(VE[j], G)];
               v2 = coordVector[source(VE[i], G)];
               vector1 = make_pair(v1.first - vt.first, v1.second - vt.second);
               vector2 = make_pair(vt.first - v2.first, vt.second - v2.second);
            } else if( target(VE[i], G) == target(VE[j], G) ) {
               vt = coordVector[target(VE[i], G)];
               v1 = coordVector[source(VE[i], G)];
               v2 = coordVector[source(VE[j], G)];
               vector1 = make_pair(v1.first - vt.first, v1.second - vt.second);
               vector2 = make_pair(vt.first - v2.first, vt.second - v2.second);
            } else
               continue;

            double angle = acos(((vector1.first * vector2.first) + (vector1.second * vector2.second))/(sqrt((vector1.first * vector1.first) + (vector1.second * vector1.second)) * (sqrt((vector2.first * vector2.first) + (vector2.second * vector2.second)))));

            //add monomial for every pair of edges to objective function
            mon = monomial(angle, VI[i], VI[j]);
            std::vector<monomial>::iterator pos;
            if( lm.size()==0 ){
               lm.push_back(mon);
               IP.add_polynomial(mon);
            }
            else {
               pos = lower_bound(lm.begin(), lm.end(), mon);
               if( pos == lm.end() || !(*pos == mon)){
                  lm.insert(pos, mon);
                  IP.add_polynomial(mon);
               }
            }
         }
      }
   }
   cerr << "end read" << endl;
}

int main(int argc, char** argv){
 
   Graph G;
   std::vector<edge_descriptor> VE;
   var_map<edge_descriptor> VM;
   std::vector<monomial> lm;
   ILP_Problem IP(Optsense_Min);

   //read Euclidian TSP instance and create objective function
   read_instance(argv[1], G, VE, VM, lm, IP);

   //add symbolic constraint for Tour
   IP.add_sym_constraint(new TOUR<Graph>(G, VM));

   //solve the problem
   IP.optimize();

   //examine the solution
   edge_descriptor e;
   std::vector<edge_descriptor> v;
   solution s = IP.get_solution();
   BGL_FORALL_EDGES(e, G, Graph)
      if(s.value(VM[e]) > 0.9)
         v.push_back(e);
   std::vector<edge_descriptor>::iterator it = v.begin();
   int pos,start;
   start = source(*it, G);
   pos=target(*it, G);
   cout<<start<<endl;
   while(pos!=start) {
      it++;
      if(it==v.end())
         it=v.begin();
      if(source(*it, G) == pos) {
         cout<<pos<<endl;
         pos=target(*it, G);
      }else if(target(*it, G) == pos) {
         cout<<pos<<endl;
         pos=source(*it, G);
      }
   }
   return 0;
};

---