tsp.cc

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <boost/property_map/property_map.hpp>
#include <fstream>
#include <iostream>
#include <scil/scil.h>
#include <scil/constraints/tour.h>
#include <string>

using namespace boost;
using namespace SCIL;
using namespace std;

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


void readTsplibFile(Graph& G, const char*fileName,  map<vertex_descriptor, double>& xp, map<vertex_descriptor, double>& yp)
{
//works on complete graphs!

   FILE*tspFile= fopen(fileName,"r");

   if(tspFile==NULL){
      std::cerr<<"TSPLIB file "<<fileName<<" could not be opened."<<std::endl;
      exit(1);
   }

   const int maxCharPerLine= 1024;
   char lineBuf[maxCharPerLine+1];
   bool dimensionFound= false;
   bool typeFound= false;
   bool coordSectionFound= false;
   int nNodes_;

   while(fgets(lineBuf,maxCharPerLine,tspFile)){
      if(strncmp(lineBuf,"DIMENSION",strlen("DIMENSION"))==0){
         if(sscanf(lineBuf,"DIMENSION : %d",&nNodes_)!=1){
            std::cerr<<"Error when reading dimension of problem."<<std::endl;
            exit(1);
         }
         dimensionFound= true;
      }
      else if(strncmp(lineBuf,"EDGE_WEIGHT_TYPE",strlen("EDGE_WEIGHT_TYPE"))
            ==0){
         if(strncmp(lineBuf,"EDGE_WEIGHT_TYPE : EUC_2D",
                  strlen("EDGE_WEIGHT_TYPE : EUC_2D"))&&
               strncmp(lineBuf,"EDGE_WEIGHT_TYPE: EUC_2D",
                  strlen("EDGE_WEIGHT_TYPE: EUC_2D"))){
            std::cerr<<"Invalid EDGE_WEIGHT_TYPE, must be EUC_2D."<<std::endl;
            exit(1);
         }
         typeFound= true;
      }
      else if(strncmp(lineBuf,"NODE_COORD_SECTION",strlen("NODE_COORD_SECTION"))
            ==0){
         coordSectionFound= true;
         break;
      }
   }

   if(!typeFound){
      std::cerr<<"Keyword EDGE_WEIGHT_TYPE not found in file "<<fileName<<".";
      std::cerr<<std::endl;
      exit(1);
   }

   if(!dimensionFound){
      std::cerr<<"Keyword DIMENSION not found in file "<<fileName<<".";
      std::cerr<<std::endl;
      exit(1);
   }

   if(!coordSectionFound){
      std::cerr<<"Keyword  NODE_COORD_SECTION not found in file "<<fileName;
      std::cerr<<"."<<std::endl;
      exit(1);
   }

   vertex_descriptor v;
   for(int i = 0; i < nNodes_; i++){
     v=add_vertex(G);
     xp[v]=0;
     yp[v]=0;
   }

   int nodeNumber;
   float x,y;

   BGL_FORALL_VERTICES(v, G, Graph) {
      if(fscanf(tspFile,"%d %f %f",&nodeNumber,&x,&y)!=3){
         std::cerr<<"Error while reading coordinates of node "<<nodeNumber<<"."<<std::endl;
         exit(1);
      } else
      xp[v]=x;
      yp[v]=y;
   }

   if(fclose(tspFile)){
      std::cerr<<"error in closing file "<<fileName<<"."<<std::endl;
      exit(1);
   }

   return;
}



int main(int argc, char ** argv) {
   Graph G;
   map<vertex_descriptor, double> xp,yp;
   readTsplibFile(G, argv[1], xp, yp);

   //calculate costs
   map<edge_descriptor, double> costs;
   vertex_iterator it1,it2;
   for(it1 = vertices(G).first; it1 != vertices(G).second; it1++){
      it2 = it1;
      it2++;
      for(;it2 != vertices(G).second; it2++){
         costs[add_edge(*it1, *it2, G).first] = sqrt((xp[*it1]-xp[*it2])*(xp[*it1]-xp[*it2])+(yp[*it1]-yp[*it2])*(yp[*it1]-yp[*it2]));
      }
   }

  ILP_Problem IP(Optsense_Min);

  var_map<edge_descriptor> VM;
  BGL_FORALL_EDGES(e,G, Graph) { 
    VM[e]=IP.add_binary_variable(costs[e]);
  }
  IP.add_sym_constraint(new TOUR<Graph>(G,VM));
  IP.optimize();
 
 
  //examine the solution
   cout<<"Optimal Solution Value: "<<IP.get_optimum()<<endl;
   cout<<"Tour: "<<endl;
   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();
   unsigned 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);
      }
   }
}

---