maxstableset.cc

#include <boost/foreach.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/depth_first_search.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/iteration_macros.hpp>
#include <deque>
#include <fstream>
#include <iostream>
#include <scil/scil.h>
#include <scil/constraints/stableset.h>
#include <string>
#include <queue>

#define INPATH ""

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;

typedef adjacency_list<vecS, vecS, bidirectionalS > diGraph;
typedef boost::graph_traits<diGraph> diGraphTraits;
typedef diGraphTraits::vertex_descriptor divertex_descriptor;
typedef diGraphTraits::edge_descriptor diedge_descriptor;


class dfs_back_visitor : public default_dfs_visitor {
   public:
      dfs_back_visitor (std::list<diedge_descriptor>* _back) : back(_back) {}
      void back_edge(diedge_descriptor e, const diGraph &G) const {
         back->push_back(e);
      }
      std::list<diedge_descriptor> *back;
};


bool readgraphComplement(const char* infile, Graph& uG, 
   map<vertex_descriptor, double>& nc, map<vertex_descriptor, int>& node_mapG) {

   diGraph G;
   std::string path(INPATH);
   path += infile;
   ifstream file(path.c_str());
   if( file.fail() ) {
      cerr << "%%%%%%%% Could not open file " << infile << endl;
      cerr << "%%%%%%%% tried " << path << endl;
      return false;
   }

   char buf[1024];
   while( file.peek() != 'p' ) { //search for 'p' indicating the size of the graph
      file.getline(buf, 1024);
   }

   //read the size of the graph and the edge data
   char c;
   int s, t, numnodes, numedges;
   double u;
   file >> c; //skip the marker 'p'
   file >> buf; //skip the format string
   file >> numnodes; //read the number of nodes in the graph
   file >> numedges; //read the number of edges in the graph

   /* Create a complete graph with two edges (v1,v2) and (v2,v1) between all 
      pairs of nodes v1 and v2 */
   for(int i = 0; i < numnodes; i++)
      add_vertex(G);
   diGraphTraits::vertex_iterator it, iit, it_end;
   for(tie(it, it_end) = vertices(G); it != it_end; it++)
      for( iit = it+1; iit != it_end; iit++) {
         add_edge(*it, *iit, G);
         add_edge(*iit, *it, G);
      }

   std::vector<divertex_descriptor> an;
   BGL_FORALL_VERTICES(v, G, diGraph) {
      nc[v] = 1.0;
      an.push_back(v);
   }

   while((file >> c) && file.good()) {
      if( c == 'e' ) {
         file >> s >> t;

         add_edge(an[s-1], an[t-1], G);
         add_edge(an[t-1], an[s-1], G);
      } 
      else if( c == 'n' ) {
         file >> s >> u;
         nc[an[s-1]] = u;
      }
      else
         file.getline(buf, 1024);
   }

   //Find all parallel edges and delete them
   std::set<diedge_descriptor> parallel_edges;
   BGL_FORALL_EDGES(e, G, diGraph) {
      diGraphTraits::in_edge_iterator ie_it, ie_it_end;
      for(tie(ie_it, ie_it_end) = in_edges(target(e ,G), G);
         ie_it != ie_it_end && *ie_it != e; ie_it++) {
         if (source(*ie_it,G) == source(e, G)) { //parallel edge exists
            parallel_edges.insert(e);
            parallel_edges.insert(*ie_it);
         }
      }
   }
   diedge_descriptor e;
   BOOST_FOREACH(e, parallel_edges)
      remove_edge(e,G);

   //Make input-graph acyclic by deleting the back edges discoverd with dfs
   std::list<diedge_descriptor> back_edges;
   dfs_back_visitor vis(&back_edges);
   depth_first_search(G, visitor(vis));
   BOOST_FOREACH(e, back_edges)
      remove_edge(e, G);

   /* Transform bidirectional Graph G in undirected Graph uG
      Since both Graphs use vecS as VertexList vertex_descriptor is the same
      type as divertex_descriptor */
   BGL_FORALL_VERTICES(v, G, diGraph)
      add_vertex(uG);
   BGL_FORALL_EDGES(e, G, diGraph)
      add_edge(source(e, G), target(e, G), uG);

   vertex_descriptor v;
   int i = 1;
   BOOST_FOREACH(v, an) {
      node_mapG[v] = i++;
         }
   return true;
}

//TODO: Make the heuristic work.
class stableset_heuristic : public primal_heuristic {
   private:
      Graph& G;
      var_map<vertex_descriptor>& VM;
      map<vertex_descriptor, double> nc;
      //compare function for priority queue
      template <typename T>
      struct greater_second : std::binary_function<T, T, bool> {
         inline bool operator() (const T& lhs, const T& rhs) {
            return (lhs.second > rhs.second);
         }
      };

   public:
      stableset_heuristic(Graph& G_, var_map<vertex_descriptor>& VM_,
         map<vertex_descriptor, double> nc_) : G(G_), VM(VM_), nc(nc_) { }

      int heuristic(subproblem& S, solution& newSol) {
         cerr << "primal heuristic" << endl;
         newSol.save_solution(S);
         map<vertex_descriptor, bool> seen;
         BGL_FORALL_VERTICES(v, G, Graph)
            seen[v] = false;
         typedef std::pair<vertex_descriptor, double> entry_t;
         typedef std::priority_queue<entry_t, std::deque<entry_t>, 
            greater_second<entry_t> > queue_t;
         queue_t prio;
         BGL_FORALL_VERTICES(v, G, Graph) {
            newSol.round_to_integer(VM[v]);
            prio.push(entry_t(v, out_degree(v, G)));
         }

         bool stop = false;
         while (!prio.empty()) {
            stop = false;
            if( newSol.value(VM[(prio.top()).first]) < .5 ) { 
               GraphTraits::adjacency_iterator it, it_end;
               for(tie(it,it_end) = adjacent_vertices((prio.top()).first,G);
                  it != it_end; it++) {
                  if( newSol.value(VM[*it]) > .5 ) {
                     stop = true;
                     break;
                  }
               }
               if( !stop )
                  newSol.set_value(VM[(prio.top()).first], 1.);
            }
            prio.pop();
         }
         return 1;
      }
};


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

   //initialize the graph
   Graph G;
   map<vertex_descriptor, int> node_mapG;
   map<edge_descriptor, double> edge_mapG;

   map<vertex_descriptor, double> nc;
   if( !readgraphComplement(argv[1], G, nc, node_mapG) )
      return 1;

   //initialize the ILP
   ILP_Problem IP(Optsense_Max);
   vertex_descriptor v;
   var_map<vertex_descriptor> VM;
   BGL_FORALL_VERTICES(v, G, Graph) 
      VM[v]=IP.add_binary_variable(nc[v]);

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

   //set the primal heuristic
   //TODO: Make the heuristic work.
//   IP.set_primal_heuristic(new stableset_heuristic(G, VM, nc));

   //solve the problem
   IP.optimize();

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

   return 0;
}

---