stableset.cc

#ifndef SCIL_STABLESET_CC
#define SCIL_STABLESET_CC

#define HIDE_THRES 5
#define LEPS 1e-8

using namespace SCIL;
using namespace boost;
using namespace std;
//using std::cout;
//using std::cerr;
//using std::endl;

template <typename Graph>
STABLESET<Graph>::STABLESET(Graph& G_, var_map<vertex_descriptor>& VM_)
   : G(G_), VM(VM_) {
   feps = 1.0e-4;
   BGL_FORALL_VERTICES_T(v, G, Graph)
      if(VM[v].type() != Vartype_Binary) {
         cerr << "stableset.cc: All variables in VM need to be";
         cerr << "Vartype_Binary" << endl;
         exit(1);
      }
};


template <typename Graph>
void STABLESET<Graph>::init(subproblem& S) {
   if(S.configuration("STABLESET_Debug_Out")=="true")
      cerr << "STABLESET::init\n";
   //Check whether Graph is undirected
   typedef typename GraphTraits::directed_category directed_category;
   if (!(is_same<directed_category,undirected_tag>::value
         || is_base_and_derived<undirected_tag, directed_category>::value)) {
      cerr << "kconnected.cc: The graph has to be undirected." << endl;
      exit(1);
   }
   //Graph has to provide out_edges()
   BOOST_CONCEPT_ASSERT((IncidenceGraphConcept<Graph>));
   //Graph has to provide vertices() and edges()
   BOOST_CONCEPT_ASSERT((VertexAndEdgeListGraphConcept<Graph>));
   //Graph has to provide clear_vertex()
   BOOST_CONCEPT_ASSERT((MutableGraphConcept<Graph>));
   //add edge inequalities
   int i = 0;
   BGL_FORALL_EDGES_T(e, G, Graph) {
      row r;
      r += VM[source(e, G)];
      r += VM[target(e, G)];
      S.add_basic_constraint(r<=1., Static);
   }

   return;
};

template <typename Graph>
int STABLESET<Graph>::exactOddCycleSeparation(subproblem& S, int maxnum) {
   if(S.configuration("STABLESET_Debug_Out")=="true")
      cerr << "STABLESET::exactOddCycleSeparation\n";

   Graph G1;
   map<vertex_descriptor, bool> node_mapG1;
   map<edge_descriptor, bool> edge_mapG1;
   map<vertex_descriptor, vertex_descriptor> GtoG1a;
   map<vertex_descriptor, vertex_descriptor> GtoG1b;
   map<vertex_descriptor, int> node_used;
   map<vertex_descriptor, vertex_descriptor> G1toG;
   map<edge_descriptor, edge_descriptor> back;
   BGL_FORALL_VERTICES_T(v, G, Graph)
      node_used[v] = 0;

   /* G1 consists of two copies G' and G" of the original graph G. */
   BGL_FORALL_VERTICES_T(v, G, Graph) {
      vertex_descriptor new_vertex;
      new_vertex = add_vertex(G1);
      GtoG1a[v] =  new_vertex;
      node_mapG1[new_vertex] = false;
      new_vertex = add_vertex(G1);
      GtoG1b[v] = new_vertex;
      node_mapG1[new_vertex] = true;
      G1toG[GtoG1a[v]]=v;
      G1toG[GtoG1b[v]]=v;
   }
   if(S.configuration("STABLESET_Debug_Out")=="true")
      cerr << "exactOddCycleSeparation: nodes added\n";

   /* G' and G" are joined by cross edges v' -- v" and v" -- v'
      for every node v of the original graph G. */
   BGL_FORALL_EDGES_T(e, G, Graph) {
      edge_descriptor new_edge;
      new_edge = (add_edge(GtoG1a[source(e,G)], GtoG1b[target(e,G)], G1)).first;
      back[new_edge] = e;
      edge_mapG1[new_edge] = true;
      new_edge = (add_edge(GtoG1b[source(e,G)], GtoG1a[target(e,G)], G1)).first;
      back[new_edge] = e;
      edge_mapG1[new_edge] = true;
   }

   if(S.configuration("STABLESET_Debug_Out")=="true")
      cerr << "exactOddCycleSeparation: edges added\n";
   //Typedefs for Dijkstra in- and output
   typedef map<edge_descriptor, double> costs_t;
   typedef map<vertex_descriptor, size_t> vertex_index_t;
   typedef map<vertex_descriptor, vertex_descriptor> pred_t;
   typedef map<vertex_descriptor, double> dist_t;

   costs_t costs;
   BGL_FORALL_EDGES_T(e, G1, Graph) {
      double co = 1. - S.value(VM[G1toG[source(e, G1)]]);
      co += S.value(VM[G1toG[target(e, G1)]]);
      costs[e] = ( co > 0. ) ? co : 0.;       
   }
   
   //vertex_index_map is needed to handle vertices stored in a std::list
   vertex_index_t vertex_index_map;
   int i = 0;
   BGL_FORALL_VERTICES_T(v, G, Graph)
      vertex_index_map[v] = i++;

   if(S.configuration("STABLESET_Debug_Out")=="true")
      cerr << "exactOddCycleSeparation: costs initialized\n";

   //Stores predecessor of vertices
   pred_t pred;
   //Stores distance of vertices to source vertex
   dist_t dist;
   int numfound = 0;

   //Property-maps required for Dijkstra
   associative_property_map<costs_t> pm_costs(costs);
   associative_property_map<vertex_index_t> pm_vertex_index(vertex_index_map);
   associative_property_map<pred_t> pm_pred(pred);
   associative_property_map<dist_t> pm_dist(dist);

   vertex_descriptor s, t;

   /* Backup Graphs to restore the nodes, which will be cleared if they are used
      more than HIDE_THRES times */
   Graph G_bak = G, G1_bak = G1;
   /* Shortest paths between two nodes v' and v" in G1 with weight < 1 
      correspond to violated odd cycle inequalities. */
   BGL_FORALL_VERTICES_T(v, G, Graph) {
      if( numfound >= maxnum )
         break;
      /* Hide nodes that have already been used in several constraints. */
      BGL_FORALL_VERTICES_T(s, G, Graph) {
         if( node_used[s] > HIDE_THRES && out_degree(s, G) != 0 ) {
            clear_vertex(s, G);
            clear_vertex(GtoG1a[s], G1);
            clear_vertex(GtoG1b[s], G1);
         }
      }
      // Degree of 0 = "hidden" vertex
      if ( out_degree(v, G) == 0 )
         continue;

      s = GtoG1a[v];
      t = GtoG1b[v];
      if(S.configuration("STABLESET_Debug_Out")=="true")
         cerr << "exactOddCycleSeparation: DIJKSTRA START:\n";
      dijkstra_shortest_paths(G1, s, weight_map(pm_costs).vertex_index_map
         (pm_vertex_index).distance_map(pm_dist).predecessor_map(pm_pred));
      double w = dist[t];

      if(S.configuration("STABLESET_Debug_Out")=="true")
         cerr << "exactOddCycleSeparation: DIJKSTRA w:" << 1 << 1 << w << endl; 
      
      if( (w < 1. - LEPS) && (pred[t] != t) ) {
         // Follow the path starting with pred[t] until s is reached,  
         bool duplicates = false;
         std::list<vertex_descriptor> cycle;
         std::list<double> cyclecosts;
         map<vertex_descriptor, bool> seen;
         BGL_FORALL_VERTICES_T(v, G, Graph)
            seen[v] = false;
         row r;

         vertex_descriptor a, b;
         double z = 0;

         /* Mark the start node. */
         seen[G1toG[t]] = true;
         node_used[G1toG[t]]++;

         while( !duplicates ) {
            edge_descriptor e = (edge(t, pred[t], G1)).first;
            
            /* Determine the nodes in the original graph. */
            a = G1toG[t];
            b = G1toG[pred[t]];

            /* Append the node of the original graph to the cycle list. */
            cycle.push_back(a);
            r += VM[a];
            z += costs[e];
            cyclecosts.push_back(costs[e]);

            if( seen[b] ) 
               duplicates = true;
            else {
               seen[b] = true;
               node_used[b]++;
            }

            t = pred[t];
         }

         /* If b is the start node, we are done.
            Otherwise we have to remove nodes from the front of the path 
            until we reach the first duplicate node.
            We also have to adapt the constraint. */
         s = G1toG[s];
         while( cycle.front() != b ) {
            t = cycle.front();
            cycle.pop_front();
            r -= VM[t];
            node_used[t]--;
            z -= cyclecosts.front();
            cyclecosts.pop_front();
         }

         /* Add the constraint to the subproblem. */
         if( z < 1. - feps ) {
            r.normalize();
            double rhs = floor((double(cycle.size()) - 1.) / 2.);
            S.add_basic_constraint(r<=rhs);
            if(S.configuration("STABLESET_Debug_Out")=="true") {
               cerr << "exactOddCycleSeparation: added constraint! rhs = ";
               cerr << 1 << 1 << rhs << endl;
            }
            numfound++;
         }
      }
   }
   if(S.configuration("STABLESET_Debug_Out")=="true") {
      cerr << "exactOddCycleSeparation: constraints added:" ;
      cerr << 1 << 0 << numfound << endl;
   }
   G = G_bak;
   G1 = G1_bak;

   return numfound;
}

template <typename Graph>
int STABLESET<Graph>::heuristicCliqueSeparation(subproblem& S, int maxnum){
   vertex_descriptor v;
   //switch the strategie for finding a start Vertex
   for(int strat = 0; strat < 2; strat++){
      vertex_descriptor startVertex = NULL;
      BGL_FORALL_VERTICES_T(v, G, Graph){
         if(startVertex == NULL){
            startVertex = v;
            continue;
         }
         if(strat == 0)
            if(fabs(S.value(VM[v])-.5) < fabs(S.value(VM[startVertex])-.5))
               startVertex = v;
            else
               if(S.value(VM[v]) < 1 && (S.value(VM[v]) > S.value(VM[startVertex])))
                  startVertex = v;
      }
      map<vertex_descriptor, bool> vertex_map;
      pair<adjacency_iterator, adjacency_iterator> vertexRange = adjacent_vertices(startVertex, G);
      adjacency_iterator adIt;
      //put all adjacent vertices in the vertex_map
      for(adIt = vertexRange.first; adIt != vertexRange.second; adIt++)
         vertex_map.insert(make_pair(*adIt,false));
      vertex_map.insert(make_pair(startVertex, true));
      bool cliqueFound = false;
      typename map<vertex_descriptor, bool> ::iterator mapIt;
      typename map<vertex_descriptor, bool> ::iterator mapIt2;
      //look for cliques
      while(! cliqueFound) {
         for(mapIt = vertex_map.begin(); mapIt != vertex_map.end(); mapIt++){
            if(mapIt->second == true)
               continue;
            vertexRange = adjacent_vertices(mapIt->first, G);
            for(mapIt2 = vertex_map.begin(); mapIt2 != vertex_map.end(); mapIt2++){
               if(mapIt2 != mapIt) {
                  for(adIt = vertexRange.first; adIt != vertexRange.second; adIt++)
                     if(*adIt == mapIt2->first)
                        break;
                  if(adIt == vertexRange.second)
                     vertex_map.erase(mapIt2);
               }
            }
            mapIt->second = true;
            break;
         }
         if(mapIt == vertex_map.end())
            cliqueFound = true;
      }
      row r;
      //create the cliqueconstraint
      for(mapIt = vertex_map.begin(); mapIt != vertex_map.end(); mapIt++)
         r += VM[mapIt->first];
      cons_obj* c = r<=1;
      if(c->violated(S)){
         S.add_basic_constraint(c);
         return 1;
      }
      delete c;
   }
   return 0;
}



template <typename Graph>
typename STABLESET<Graph>::status STABLESET<Graph>::feasible(solution& S) {
   if(S.originating_subproblem()->configuration("STABLESET_Debug_Out")=="true")
      cerr << "STABLESET::feasible\n";
   /* A solution is feasible if the edge inequalities
    * x_i + x_j <= 1. 
    * for all edges ij \in E are satisfied.
    * Since these inequalities are added in init() and the solution 
    * is integer, the test below is obsolete and could be skipped.
    */
   BGL_FORALL_EDGES_T(e, G, Graph) {
      if( S.value(VM[source(e, G)]) && S.value(VM[target(e, G)]))
         return infeasible_solution;
   }

   return feasible_solution;
}

template <typename Graph>
typename STABLESET<Graph>::status 
STABLESET<Graph>::standard_separation(subproblem& S) {
   if(S.configuration("STABLESET_Debug_Out")=="true")
      cerr << "STABLESET::standard_separation\n";
   int num = 2000;
   if( heuristicCliqueSeparation(S, num) )
      return constraint_found;
   else {
      if( exactOddCycleSeparation(S, num) )
         return constraint_found;
      else
         return no_constraint_found;
   }
}

template <typename Graph>
typename STABLESET<Graph>::status
STABLESET<Graph>::fast_separation(subproblem& S) {
   if(S.configuration("STABLESET_Debug_Out")=="true")
      cerr << "STABLESET::fast_separation\n";
   int num = 2000;
   if( heuristicCliqueSeparation(S, num) )
      return constraint_found;
   else
      return no_constraint_found;
}

template <typename Graph>
void STABLESET<Graph>::info() {
   cout<<"Configurations:\n";
   cout<<"STABLESET_Debug_Out [true|false]\n";
};
#endif

---