perfectBmatching.cc

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <scil/constraints/matching.h>
#include <scil/scil.h>

#include<iostream>
#include<fstream>
#include<string>

#define INPATH ""

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

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;

struct coord {
   double x, y;
};

//this function reads a graph in TSPLIB format with euclidian edge weights
//if writeout is true the graph is output in edge format
bool readgraphFromTSPLIB(bool writeout, int base, const char* infile, Graph& G, 
                         map<vertex_descriptor, int>& nc,
                         map<vertex_descriptor, int>& node_mapG,
                         map<edge_descriptor, double>& edge_mapG) {
   string path(INPATH);
   path += infile;
   ifstream file(path.c_str());
   if( file.fail() ) {
      cerr << "%%%%%%%% Could not open file " << infile << endl;
      cerr << "%%%%%%%% tried " << path << endl;
      exit(1);
   }
   double sign = 1.;
   if( base == -1 ) {
      base = 1;
      sign = -1.;
   }

   int count = 0;
   std::string in;
   do {
      getline(file, in);
   } while( (in.compare(0, 8, "DIMENSION", 0, 8) != 0) && (count++ < 100) ) ;
   int pos = in.find_last_of(' ');
   int numnodes = atoi(in.substr(pos+1, in.length()-pos-1).c_str());

   do {
      getline(file, in);
   } while( (in.compare(0, 9, "NODE_COORD_SECTION", 0, 9) != 0) && (count++ < 100) ) ;

   vertex_descriptor* nodelist = new vertex_descriptor[numnodes];

   for( int i = 0; i < numnodes; i++) {
      nodelist[i] = add_vertex(G);
      node_mapG[nodelist[i]] = i;
   }

   map<vertex_descriptor, coord> coords;

   BGL_FORALL_VERTICES(v, G, Graph)
      nc[v] = 0;
   int s;
   for( int i = 0; i < numnodes; i++) {
      file >> s >> coords[nodelist[i]].x >> coords[nodelist[i]].y;
      nc[nodelist[i]] = 1;
   }

   double xd, yd;
   int dist;
   for( int i = 0; i < numnodes; i++ ) {
      for( int j = i + 1; j < numnodes; j++ ) {
         xd = coords[nodelist[i]].x - coords[nodelist[j]].x;
         yd = coords[nodelist[i]].y - coords[nodelist[j]].y;
         dist = int(sqrt(xd * xd + yd * yd) + .5);
         edge_mapG[(add_edge(nodelist[i], nodelist[j], G)).first] = dist;
      }
   }

   if(writeout) {
      cout << num_vertices(G) << " " << num_edges(G) << endl;
      BGL_FORALL_EDGES(e, G, Graph) {
         cout << node_mapG[source(e, G)] << " " << node_mapG[target(e, G)];
         cout << " " << edge_mapG[e] << endl;
      }
      return false;
   }

   delete[]nodelist;
   return true;
}

bool readgraph(int base, const char* infile, Graph& G, 
               map<vertex_descriptor, int>& nc,
               map<vertex_descriptor, int>& node_mapG,
               map<edge_descriptor, double>& edge_mapG) {
   string path(INPATH);
   path += infile;
   ifstream file(path.c_str());
   if( file.fail() ) {
      cerr << "%%%%%%%% Could not open file " << infile << endl;
      cerr << "%%%%%%%% tried " << path << endl;
      exit(1);
   }
   double sign = 1.;
   if( base == -1 ) {
      base = 1;
      sign = -1.;
   }

   int s, t, numnodes, numedges;
   double u;
   file >> numnodes; //read the number of nodes in the graph
   file >> numedges; //read the number of edges in the graph

   vertex_descriptor* nodelist = new vertex_descriptor[numnodes];

   for( int i = 0; i < numnodes; i++) {
      nodelist[i] = add_vertex(G);
      node_mapG[nodelist[i]] = i;
   }

   BGL_FORALL_VERTICES(v, G, Graph)
      nc[v] = 0;
   for( int i = 0; i < numnodes; i++) {
      //file >> s;
      nc[nodelist[i]] = 1;
   }

   while((file >> s >> t >> u) && file.good()) {
      edge_mapG[(add_edge(nodelist[s-base],nodelist[t-base],G)).first]=sign*u;
   }

   delete[]nodelist;
   return true;
}


int main(int argc, char** argv)
{
   if( argc < 2 ) {
      cerr << "Please give the name of an input file!" << endl;
      return 1;
   }
   cerr << "Start Matching-Example";

   int base = 1;
   int filearg = 1;
   bool writeout = false;

   if( (argc == 3) ) {
      cerr << argv[1] << endl;
      if( argv[1][0] == 'p' ) {
      writeout = true;
      filearg = 2;
   }
   else 
      return 1;
   }

   //initialize the graph
   Graph G;
   map<vertex_descriptor, int> node_mapG;
   map<edge_descriptor, double> edge_mapG;
   map<vertex_descriptor, int> nc;
   //if( !readgraph(base, argv[1], G, nc) )
   if( !readgraphFromTSPLIB(writeout, base, argv[filearg], G, nc, node_mapG,
                            edge_mapG) )
      return 1;

   //initialize the ILP
   ILP_Problem IP(Optsense_Min);
   edge_descriptor e;
   var_map<edge_descriptor> VM;
   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 constraint
   IP.add_sym_constraint(new MATCHING<Graph>(G, nc, VM));

   //solve the problem
   IP.optimize();

   //output the solution
   solution sol =  IP.get_solution();
   double sum = 0.;
   cout << "\nName: " << argv[1] << endl;
   BGL_FORALL_EDGES(e, G, Graph) {
      sum += double(edge_mapG[e]) * sol.value(VM[e]);
      if( sol.value(VM[e]) > .9 ) {
         cout << node_mapG[source(e, G)] + base << " -- ";
         cout << node_mapG[target(e, G)] + base << endl;
      }
   }
   cout << endl;
   cout.precision(2);
   cout << fixed << "solution value: " << sum << endl;

   return 0;
}

---