dcmst.cc

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <iostream>
#include <fstream>
#include <list>
#include <scil/scil.h>
#include <scil/constraints/matching.h>
#include <scil/constraints/spantree.h>
#include <string>


#define OUTPUT
#define INPATH ""
#define OUTPATH ""
#define EPS 1e-8

using std::cout;
using std::cerr;
using std::endl;
using std::ifstream;
using std::atoi;

using namespace boost;
using namespace SCIL;

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;

bool readExtTsplibFile(const char*fileName, Graph& G, 
   map<vertex_descriptor, int>& nc, map<vertex_descriptor, int>& node_mapG,
   map<edge_descriptor, double>& edge_mapG) {
   FILE* dcmstFile = fopen(fileName, "r");
   if(dcmstFile == NULL){
      cerr<<"dcmstFile "<<fileName<<" could not be opened."<<endl;
      exit(1);
   }

   const unsigned int maxCharPerLine = 1024;
   char lineBuf[maxCharPerLine + 1];
   bool dimensionFound = false;
   bool formatFound = false;
   bool constraintSectionFound = false;
   bool dataSectionFound = false;
   bool allInteger = true;
   unsigned int nNodes = 0;
   unsigned int nEdges = 0;

   while(fgets(lineBuf,maxCharPerLine,dcmstFile)){
      //the first keyword that has to be found is DIMENSION
      if(strncmp(lineBuf,"DIMENSION",strlen("DIMENSION"))==0){
         if(sscanf(lineBuf,"DIMENSION : %d",&nNodes)!=1){
            cerr<<"Error when reading dimension of problem."<<endl;
            exit(1);
         }
         dimensionFound= true;
      }
      else if(strncmp(lineBuf,"EDGE_DATA_FORMAT",
         strlen("EDGE_DATA_FORMAT"))==0) {
         if(strncmp(lineBuf,"EDGE_DATA_FORMAT : EDGE_WEIGHT_LIST",
            strlen("EDGE_DATA_FORMAT : EDGE_WEIGHT_LIST"))) {
            cerr<<"Invalid EDGE_DATA_FORMAT, must be EDGE_WEIGHT_LIST."<<endl;
            exit(1);
         }
         formatFound= true;
      }
      else if(strncmp(lineBuf,"NODE_DEGREE_CONSTRAINT_SECTION",
         strlen("NODE_DEGREE_CONSTRAINT_SECTION"))==0) {
         constraintSectionFound= true;
         int tempA;
         int tempB;

         while(true) {
            if(fscanf(dcmstFile,"%d %d",&tempA,&tempB)!=2) {
               if(tempA== -1)
                  break;
               else {
                  cerr<<"Error while reading degree constraints for nodes.";
                  cerr<<endl;
                  exit(1);
               }
            }
         }
      }
      else if(strncmp(lineBuf,"EDGE_DATA_SECTION", 
         strlen("EDGE_DATA_SECTION"))==0) {
         dataSectionFound= true;
         int tempA;
         int tempB;
         double tempC;
         nEdges= 0;
         while(true) {
            if(fscanf(dcmstFile,"%d %d %lf",&tempA,&tempB,&tempC)==3)
               ++nEdges;
            else {
               if(tempA== -1)
                  break;
               else {
                  cerr<<"Error while reading edges."<<endl;
                  exit(1);
               }
            }
         }
      }
   }

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

   if(!formatFound){
      cerr<<"Keyword EDGE_DATA_FORMAT not found in file "<<fileName<<endl;
      exit(1);
   }

   if(!constraintSectionFound){
      cerr<<"Keyword NODE_DEGREE_CONSTRAINT_SECTION not found in file ";
      cerr<<fileName<<endl;
      exit(1);
   }

   if(!dataSectionFound){
      cerr<<"Keyword EDGE_DATA_SECTION not found in file "<<fileName<<endl;
      exit(1);
   }

   //rewind the file for getting data
   if(fseek(dcmstFile,0L,SEEK_SET)){
      cerr<<"Error while rewinding file "<<fileName<<"."<<endl;
      exit(1);
   }

   //initialize the nodes
   vertex_descriptor* nodes = new vertex_descriptor[nNodes];
   for(unsigned int i=0; i<nNodes; ++i) {
      nodes[i] = add_vertex(G);
      node_mapG[nodes[i]] = i;
   }

   //initialize the degree constraints of the nodes with nNodes-1
   //(that means no constraint)
   BGL_FORALL_VERTICES(v, G, Graph)
      nc[v] = nNodes-1;

   while(fgets(lineBuf,maxCharPerLine,dcmstFile)) {
      if(strncmp(lineBuf,"NODE_DEGREE_CONSTRAINT_SECTION", 
         strlen("NODE_DEGREE_CONSTRAINT_SECTION"))==0) {
         int tempA;
         int tempB;
         while(true) {
            if(fscanf(dcmstFile,"%d %d",&tempA,&tempB)==2) {
               //save the degree constraints for those nodes that have a real 
               //restriction
               if(0 < tempB && tempB < (int)nNodes-1)
                  nc[nodes[tempA]]= tempB;
            }
            else {
               if(tempA==-1)
                  break;
               else {
                  cerr<<"Error while reading degree constraints for nodes.";
                  cerr<<endl;
                  exit(1);
               }
            }
         }
      }
      else if(strncmp(lineBuf,"EDGE_DATA_SECTION", 
         strlen("EDGE_DATA_SECTION"))==0) {
         int tempA;
         int tempB;
         double tempC;
         int i=0;
         while(true) {
            if(fscanf(dcmstFile,"%d %d %lf",&tempA,&tempB,&tempC)==3) {
               if(fabs(tempC - floor(tempC)) >= EPS)
                  allInteger = false;
               else
                  tempC = round(tempC); //make it really integer
               //set the edges and store their costs in the temporary array
               edge_descriptor new_edge;
               new_edge = (add_edge(nodes[tempA], nodes[tempB], G)).first;
               edge_mapG[new_edge] = tempC;
            }
            else {
               if(tempA== -1)
                  break;
               else {
                  cerr<<"Error while reading edges."<<endl;
                  exit(1);
               }
            }
         }
      }
   }

   //we remove trivial node degree constraints
   BGL_FORALL_VERTICES(v, G, Graph) {
      if( (int)out_degree(v, G) <= nc[v] )
         nc[v] = nNodes-1;
   }

   delete[]nodes;
   if(fclose(dcmstFile)){
      cerr<<"error in closing file "<<fileName<<endl;
      exit(1);
   }
   return true;
}

/*------------------------------------------------------------------*/
struct edgeCostsPair{
   const edge_descriptor e;
   const double c;
   edgeCostsPair(const edge_descriptor _e, const double _c) : e(_e), c(_c) {}
   bool operator<(const edgeCostsPair& a) const{
      return c < a.c;
   }
};

class dcmst_heuristic : public primal_heuristic {
   private:
      Graph& G;
      var_map<edge_descriptor>& VM;
      map<vertex_descriptor, int>& dc; //degree constraints

   public:
      dcmst_heuristic(Graph& G_, var_map<edge_descriptor>& VM_, 
         map<vertex_descriptor, int>& dc_): G(G_), VM(VM_), dc(dc_) { }

      int heuristic(subproblem& S, solution& newSol) {
         /* extended Kruskal heuristic */
         map<vertex_descriptor, int> act_degree;
         BGL_FORALL_VERTICES(v, G, Graph)
            act_degree[v] = 0;

         std::list<edgeCostsPair> edgesPriority;

         newSol.save_solution(S);
         BGL_FORALL_EDGES(e, G, Graph) {
            newSol.set_value(VM[e], 0.);
            edgesPriority.push_back(edgeCostsPair(e, 1. - S.value(VM[e])));
         }
         edgesPriority.sort();

         unsigned int E_T_counter = 0;
         typedef map<vertex_descriptor,set<vertex_descriptor> >vertexPartition;
         vertexPartition V_T; //= new vertexPartition();
         //Initialize the vertex_partition
         BGL_FORALL_VERTICES(v, G, Graph) {
            //(V_T)[v] = set<vertex_descriptor>;
            (V_T)[v].insert(v);
         }

         int cr = 0;
         for(std::list<edgeCostsPair>::const_iterator it=edgesPriority.begin();
            it != edgesPriority.end(); ++it) {
            const vertex_descriptor s = source(it->e, G);
            const vertex_descriptor t = target(it->e, G);

            if( act_degree[s] + 1 <= dc[s] && act_degree[t] + 1 <= dc[t] 
               && ((V_T)[s]) != ((V_T)[t]) ) {
               newSol.set_value(VM[it->e], 1.);
               ++E_T_counter;
               ++act_degree[s];
               ++act_degree[t];
               //Union sets of vertices s with sets of vertices t
               BGL_FORALL_VERTICES(v, G, Graph) {
                  set<vertex_descriptor> tmp = (V_T)[t];
                  (V_T)[v].insert(tmp.begin(), tmp.end());
                  set<vertex_descriptor>::iterator it;
                  //Update the pointers of the vertices in V_T[t]
                  for(it = tmp.begin(); it != tmp.end(); it++)
                     (V_T)[*it] = (V_T)[s];
               }
            }

            if(E_T_counter == num_vertices(G) - 1)
               break;
         }

         //FIXME
         /*
         BGL_FORALL_VERTICES(v, G, Graph)
            delete (*V_T)[v];
            */
         //delete V_T;

         if(E_T_counter != num_vertices(G) - 1) {
            //cout << "extendedKruskal could not find a dcst" << endl;
            return 0;
         }

         return 1;
      }
};

/*------------------------------------------------------------------*/


int main(int argc, char** argv)
{

   //initialize the input and output file names
   if( argc < 2 ) {
      cerr << "Please give the name of an input file!" << endl;
      return 1;
   }

   char* pch = strrchr(argv[argc - 1],'/');
   std::string outname;
   if( pch != NULL ) {
      outname += pch;
      outname.erase(0, 1);
   }
   else
      outname += argv[argc - 1];

   std::string logname = OUTPATH + outname;
   outname = OUTPATH + outname + ".res";
   //------------------------------------------------------

   //initialize the graph
   Graph G;
   map<vertex_descriptor, int> node_mapG;
   map<edge_descriptor, double> edge_mapG;
   map<vertex_descriptor, int> nc;
   if( !readExtTsplibFile(argv[1], G, nc, node_mapG, edge_mapG) )
      return 1;
   //------------------------------------------------------

   //initialize the ILP
   ILP_Problem IP(Optsense_Min);
   var_map<edge_descriptor> VM;
   std::list<var> vl;
   var va;

   BGL_FORALL_EDGES(e, G, Graph) {
      vl.push_back(IP.add_binary_variable(edge_mapG[e]));
      VM[e] = vl.back();
   }
   //------------------------------------------------------

   //add the symbolic constraints
   IP.add_sym_constraint(new SpanTree<Graph>(G, VM));

   IP.add_sym_constraint(new MATCHING<Graph>(G, nc, VM, false));
   //------------------------------------------------------
   IP.set_primal_heuristic(new dcmst_heuristic(G, VM, nc));
   //------------------------------------------------------

   //solve the problem
   IP.optimize();
   //------------------------------------------------------

   //write the solution and running time to the result file 
#ifdef OUTPUT
   ofstream outFile;
   outFile.open(outname.c_str());
   if( outFile.fail() ) {
      cerr << "File " << outname << " could not be opened." << endl;
      exit(1);
   }

   //output the solution
   outFile << "\nName: " << argv[argc - 1] << endl;

   int base = 1;
   double sum = 0.;
   solution sol = IP.get_solution();

   BGL_FORALL_EDGES(e, G, Graph) {
      sum += double(edge_mapG[e]) * sol.value(VM[e]);
      if( sol.value(VM[e]) > 1. - EPS ) {
         if( node_mapG[source(e, G)] < node_mapG[target(e, G)] ) {
            outFile << node_mapG[source(e, G)] + base;
            outFile << " -- " << node_mapG[target(e, G)] + base << endl;
         } else {
            outFile << node_mapG[target(e, G)] + base;
            outFile << " -- " << node_mapG[source(e, G)] + base << endl;
         }
      }
   }
   outFile << endl;
   outFile.precision(2);
   outFile << fixed << "solution value: " << sum << endl;

#endif //OUTPUT
   //------------------------------------------------------
   return 0;
}

---