cut.cc

#ifndef SCIL_CUT_CC
#define SCIL_CUT_CC

#define HIDE_THRES 5

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

/*  This visitor records the distances of the vertices and stores the order
    in which the vertices are discovered in the list bfs_order
 */
template <typename Graph>
class CUT<Graph>::bfs_dist_visitor : public default_bfs_visitor {
   public:
      bfs_dist_visitor(map<vertex_descriptor, int>* _distmap,
                       list<vertex_descriptor>* _bfs_order,
                       int start_vertex)
         : distmap(_distmap), bfs_order(_bfs_order) {
         (*distmap)[start_vertex] = 0;
         bfs_order->clear();
      }
      /* The source of e is the predecessor of the target of e, which is the
         newly discovered vertex.
       */
      void tree_edge(edge_descriptor e, const Graph &G) const {
         vertex_descriptor s = source(e, G), t = target(e, G);
         (*distmap)[t] = (*distmap)[s]+1;
      }
      void discover_vertex(vertex_descriptor u, const Graph &G) const {
         bfs_order->push_back(u);
      }
      map<vertex_descriptor, int> *distmap;
      list<vertex_descriptor> *bfs_order; 
};

template <typename Graph>
void CUT<Graph>::undirected_bfs(Graph &H, vertex_descriptor r,
                                map<vertex_descriptor, int> &dst,
                                map<vertex_descriptor, edge_descriptor> &pred) {
   vertex_descriptor v;
   BGL_FORALL_VERTICES_T(v, H, Graph) {
      dst[v] = -1;
   }
   std::queue<vertex_descriptor> q;
   vertex_descriptor s, t;
   edge_descriptor e;
   dst[r] = 0;
   q.push(r);
   while( !q.empty() ) {
      s = q.front();
      q.pop();
      typename GraphTraits::out_edge_iterator oe_it, oe_it_end;
      for(tie(oe_it, oe_it_end)=out_edges(s,H); oe_it != oe_it_end; oe_it++) {
         t = target(*oe_it, H);
         if( dst[t] == -1 ) {
            dst[t] = dst[s] + 1;
            pred[t] = *oe_it;
            q.push(t);
         }
      }
   }
}

template <typename Graph, typename vertex_descriptor, typename edge_descriptor>
bool color_graph(Graph& H, vertex_descriptor& r, 
                 map<edge_descriptor,double>& cut, 
                 map<vertex_descriptor, bool>& color) {

   typedef graph_traits<Graph> GraphTraits;

   map<vertex_descriptor, int> dst;
   vertex_descriptor v;
   BGL_FORALL_VERTICES_T(v, H, Graph) {
      dst[v] = -1;
      color[v] = false;
   }

   std::queue<vertex_descriptor> q;
   vertex_descriptor s, t;
   edge_descriptor e;
   dst[r] = 0;
   q.push(r);
   while( !q.empty() ) {
      s = q.pop();
      typename GraphTraits::out_edge_iterator oe_it, oe_it_end;
      for(tie(oe_it, oe_it_end)=out_edges(s, H); oe_it != oe_it_end; oe_it++) {
         if( cut[*oe_it] < -.5 )
            continue;
         (s == source(*oe_it,H)) ? t = target(*oe_it,H) : t = source(*oe_it,H);
         if( dst[t] == -1 ) {
            dst[t] = dst[s] + 1;
            if(cut[*oe_it] > .5)
               color[t] = !color[s];
            else
              color[t] = color[s];
            q.append(t);
         } else {
            if( (cut[*oe_it] > .5 && color[t] == color[s]) 
                 || cut[*oe_it] <= .5 && color[t] != color[s] )
               return false;
         }
      }
   }
   return true;
}

template <typename Graph>
CUT<Graph>::CUT(Graph &G_, row_map<edge_descriptor> &VM_, int nc_)
   : G(G_), VM(VM_), updateConfiguration(false), nc(nc_), initialized(false) {
   edge_descriptor e;
   feps = .01;

   //Check whether Graph is undirected
   typedef typename GraphTraits::directed_category directed_category;
   bool graph_is_undirected = is_same<directed_category,undirected_tag>::value
      || is_base_and_derived<undirected_tag, directed_category>::value;
   BOOST_ASSERT(graph_is_undirected);

   BGL_FORALL_EDGES_T(e, G, Graph) {
      row_entry::list_iterator it;
      foreach( it, VM[e] )
         if( it->Var != NULL )
         if( it->Var.type() != Vartype_Binary ) {
             cerr << "cut.cc: All variables in VM need to be Vartype_Binary";
             cerr << endl;
             exit(1);
         }
   }
   //Graph has to provide vertices() and edges()
   BOOST_CONCEPT_ASSERT((VertexAndEdgeListGraphConcept<Graph>));
   //Graph has to provide out_edges()
   BOOST_CONCEPT_ASSERT((IncidenceGraphConcept<Graph>));
   //Graph has to provide clear_vertex()
   BOOST_CONCEPT_ASSERT((MutableGraphConcept<Graph>));

   //Compute connected components of G
   vector<int> component(num_vertices(G));
   int num_components = connected_components(G, &component[0]);

   //cout<<"Connected Components: "<<num_components<<endl;
   for(int i = 0; i < num_components; i++){
      row_map<edge_descriptor> tmpVM;
      Graph tmp(num_vertices(G));
      graph_components.push_back(tmp);
      Graph& U = graph_components.back();
      BGL_FORALL_EDGES_T(e, G, Graph)
         tmpVM[add_edge(source(e,G), target(e,G), U).first]=VM[e];
      for(int j = num_vertices(G) - 1; j >= 0 ; j--){
         if(component[j] != i){
            clear_vertex(vertices(U).first[j], U);
            remove_vertex(vertices(U).first[j], U);
         }
      }
      VML.push_back(tmpVM);
   }
}

template <typename Graph>
int CUT<Graph>::generate_triangle_inequalities(subproblem& S) {
   if( S.configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::generate_triangle_inequalities\n";

   unordered_map<std::string, tuple<edge_descriptor, edge_descriptor, edge_descriptor> > tris;

   //Check if the Graph has at least 3 edges
   if(num_edges(G)<3)
      return 0;

   map<vertex_descriptor,int> vertex_index_map;
   int i = 0;
   BGL_FORALL_VERTICES_T(v, G, Graph)
      vertex_index_map[v] = i++;
   stringstream ss;
   int index1,index2,index3;

   typename GraphTraits::edge_iterator e_it, f_it, g_it, e_it_end;

   //Search triangles by checking all combinations of three adjacent edges
   for(tie(e_it, e_it_end) = edges(G); e_it != e_it_end; e_it++){
      for (typename GraphTraits::out_edge_iterator f_it = out_edges(target(*e_it, G), G).first; f_it != out_edges(target(*e_it, G), G).second; f_it++){
         for (typename GraphTraits::out_edge_iterator g_it = out_edges(target(*f_it, G), G).first; g_it != out_edges(target(*f_it, G), G).second; g_it++){
            if (target(*g_it, G) == source(*e_it, G)){
               index1 = min(vertex_index_map[target(*e_it, G)], min(vertex_index_map[target(*f_it, G)], vertex_index_map[target(*g_it, G)]));
               index3 = max(vertex_index_map[target(*e_it, G)], max(vertex_index_map[target(*f_it, G)], vertex_index_map[target(*g_it, G)]));
               index2 = vertex_index_map[target(*e_it, G)] + vertex_index_map[target(*f_it, G)] + vertex_index_map[target(*g_it, G)] - index1 - index3;
               ss.str("");
               ss.clear();
               ss<<index1<<"#"<<index2<<"#"<<index3;
               tris.insert(make_pair(ss.str(), tuple <edge_descriptor, edge_descriptor, edge_descriptor> (*e_it, *f_it, *g_it)));
            }
         }
      }
   }

   //create the inequalities and add them
   typename unordered_map<std::string, tuple<edge_descriptor, 
            edge_descriptor, 
            edge_descriptor> >::const_iterator triangle;
   int gen = 0;
   foreach(triangle, tris) {
      row r1, r2, r3, r4;
      r1 += VM[get<0>(triangle->second)];
      r1 -= VM[get<1>(triangle->second)];
      r1 -= VM[get<2>(triangle->second)];

      r2 += VM[get<0>(triangle->second)];
      r2 += VM[get<1>(triangle->second)];
      r2 += VM[get<2>(triangle->second)];

      r3 -= VM[get<0>(triangle->second)];
      r3 += VM[get<1>(triangle->second)];
      r3 -= VM[get<2>(triangle->second)];

      r4 -= VM[get<0>(triangle->second)];
      r4 -= VM[get<1>(triangle->second)];
      r4 += VM[get<2>(triangle->second)];

      S.add_basic_constraint(r1<=0., Dynamic);
      S.add_basic_constraint(r2<=2., Dynamic);
      S.add_basic_constraint(r3<=0., Dynamic);
      S.add_basic_constraint(r4<=0., Dynamic);
      gen += 4;
   }
   return 0;
}

template <typename Graph>
void CUT<Graph>::init(subproblem& S) {
   if( S.configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::init\n";
   //Check whether Graph is undirected
   typedef typename GraphTraits::directed_category directed_category;
   bool graph_is_undirected = is_same<directed_category,undirected_tag>::value
      || is_base_and_derived<undirected_tag, directed_category>::value;
   BOOST_ASSERT(graph_is_undirected);
   generate_triangle_inequalities(S);
   return;
};

template <typename Graph>
void CUT<Graph>::initAuxGraph(solution S) {
   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::initAuxGraph\n";

   G1.clear();
   back.clear();
   GtoG1a.clear();
   GtoG1b.clear();
   G1toG.clear();
   isCrossEdge.clear();
   vertex_descriptor v;
   edge_descriptor e;

   /* G1 consists of two copies G' and G" of the original graph G. */
   BGL_FORALL_VERTICES_T (v, G, Graph) {
      GtoG1a[v] = add_vertex(G1);
      GtoG1b[v] = add_vertex(G1);
      G1toG[GtoG1a[v]] = v;
      G1toG[GtoG1b[v]] = v;
   }

   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "initAuxGraph: nodes added\n";

   vertex_descriptor s, t;

   /* G' and G" each contain copies of the original edges
      and 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 f;
      s = source(e,G);
      t = target(e,G);
      f = (add_edge(GtoG1a[s], GtoG1a[t], G1)).first;
      back[f] = e;
      isCrossEdge[f] = false;
      f = (add_edge(GtoG1b[s], GtoG1b[t], G1)).first;
      back[f] = e;
      isCrossEdge[f] = false;
      f = (add_edge(GtoG1a[s], GtoG1b[t], G1)).first;
      back[f] = e;
      isCrossEdge[f] = true;
      f = (add_edge(GtoG1b[s], GtoG1a[t], G1)).first;
      back[f] = e;
      isCrossEdge[f] = true;
   }

   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "initAuxGraph: G1 initialized\n";

   initialized = true;
   return;
}

template <typename Graph>
int CUT<Graph>::exactCycleSeparation(solution& S, 
                                     int maxnum, 
                                     std::list<cons_obj*>& newCons) {
   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::exactCycleSeparation\n";

   initAuxGraph(S);

   //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, double> dist_t;
   typedef map<vertex_descriptor, vertex_descriptor> pred_t;

   vertex_descriptor v;
   map<vertex_descriptor, int> node_used;
   BGL_FORALL_VERTICES_T(v, G, Graph)
      node_used[v] = 0;

   double leps = .001;
   row lastCon = nil;
 
   edge_descriptor e;
   costs_t costs;
   BGL_FORALL_EDGES_T(e, G1, Graph) {
      if( isCrossEdge[e] )
         costs[e] = 1. - S.value(VM[back[e]]);
      else
         costs[e] = S.value(VM[back[e]]);
      if( costs[e] < 0. ) costs[e] = 0.;
      //assert(costs[e]>=0.);
      //assert(costs[e]<=1.);
   }
   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "exactCycleSeparation: costs initialized\n";

   /* Shortest paths between two nodes v' and v" in G1 with weight < 1 
      correspond to violated odd cycle inequalities. */
   vector<vertex_descriptor> ln;
   BGL_FORALL_VERTICES_T(v, G, Graph)
      ln.push_back(v);
   int nn = num_vertices(G);
   if( nc > 0 ) {
      random_shuffle(ln.begin(), ln.end());
      nn = num_vertices(G)/nc;
   }

   //Vertex_index_map is needed if vertices are stored in a std::list
   vertex_index_t vertex_index_map;
   int i = 0;
   BGL_FORALL_VERTICES_T(v, G1, Graph)
      vertex_index_map[v] = i++;
   int numfound = 0;
   //Stores predecessor of vertices
   pred_t pred;
   //Stores distance of vertices to source vertex
   dist_t dist;
   vertex_descriptor s,t;

   //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);

   int ntc = 0;

   //Create copy of Graph G
   Graph G2;
   map<vertex_descriptor, vertex_descriptor> GtoG2;
   BGL_FORALL_VERTICES_T(v, G, Graph)
      GtoG2[v] = add_vertex(G2);
   BGL_FORALL_EDGES_T(e, G, Graph)
      add_edge(GtoG2[source(e,G)], GtoG2[target(e,G)], G2);

   BOOST_FOREACH(v, ln) {
      if( numfound >= maxnum || ntc >= nn )
         break;
      /* Hide nodes that have already been used in several constraints. Since
         the adjacency_list does not provide functionality to hide vertices
         the edges adjacent to the vertices are removed, so that the nodes will
         not be used again by Dijkstra and the vertex_iterators remain intact.
       */
      BGL_FORALL_VERTICES_T(s, G, Graph) {
         if( node_used[s] > HIDE_THRES && out_degree(GtoG2[s], G2) != 0 ) {
            clear_vertex(GtoG2[s], G2);
            clear_vertex(GtoG1a[s], G1);
            clear_vertex(GtoG1b[s], G1);
         }
      }

      // Degree of 0 = hidden vertex
      if( out_degree(GtoG2[v], G2) == 0 ) 
         continue;

      ntc++;

      s = GtoG1a[v];
      t = GtoG1b[v];

      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.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
         cerr << "exactCycleSeparation: DIJKSTRA w:" << w << endl;
      if( (w < 1. - leps) && (pred[t] != t) ) {
         // Follow the path starting with pred[t] until s is reached
         bool duplicates = false;
         list<edge_descriptor> cycle;
         map<edge_descriptor, bool> oddset;
         BGL_FORALL_EDGES_T(e, G, Graph)
            oddset[e] = false;
         map<vertex_descriptor, bool> seen;
         BGL_FORALL_VERTICES_T(v, G, Graph)
            seen[v] = false;
         row r;
         vertex_descriptor a, b;
         double rhs = -1.;

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

         edge_descriptor f;
         while( !duplicates ) {
            e = (edge(t, pred[t], G1)).first;
            f = back[e];
            /* Append the edge of the original graph to the cycle list. */
            cycle.push_back(f);

            /* Mark the original edge as belonging to the odd set or not. */
            if ( isCrossEdge[e] ) {
               r += VM[f];
               rhs += 1.;
               oddset[f] = true;
            }
            else
               r -= VM[f];

            /* Determine the nodes in the original graph. */
            a = G1toG[t];
            b = G1toG[pred[t]];
            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 edges from the front of the path 
            until we reach the first duplicate node.
            We also have to adapt the constraint. */
           s = G1toG[s];
         if( b != s ) {
            do {
               e = cycle.front();
               cycle.pop_front();            
               if( oddset[e] ) {
                  r -= VM[e];
                  rhs -= 1.;
               }
               else 
                  r += VM[e];

               node_used[source(e, G)]--;
               node_used[target(e, G)]--;
            } while((source(e, G) != b) && (target(e, G) != b));
            node_used[b]++;
         }
         /* Add the constraint to the subproblem. */
         r.normalize(true);
         cons_obj* c = (r<=rhs);
         //FIXME compare
         if( c->violation(S) > feps && compare(r, lastCon) != 0 ) {
            lastCon = r;
            newCons.push_back(r<=rhs);
            if(S.originating_subproblem()->configuration("CUT_Debug_Out")=="true")
               if( rhs <= -1.)
                  cerr << r << " <= " << rhs << endl;
            numfound++;
            if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
               cerr << "exactCycleSeparation: added constraint! rhs = " << 1 << 1 << rhs << endl;
         }
         delete c;
      }
   }
   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "exactCycleSeparation: constraints added:" << 1 << 0 << numfound << endl;

   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "numfound: " << numfound << endl;
   return numfound;
}

template <typename Graph>
int CUT<Graph>::forestCycleSeparation(solution& S, int maxnum, 
      list<cons_obj*>& newCons) {
   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::forestCycleSeparation\n";

   int numfound = 0;
   
   typename list<row_map<edge_descriptor> >::iterator it2 = VML.begin();

   for(typename list<Graph>::iterator it=graph_components.begin(); it != graph_components.end(); it++){

      Graph& G = *it;
      row_map<edge_descriptor> VM = *it2;

      //Typedefs for MST INPUT
      typedef map<edge_descriptor, double> weight_map_t;
      typedef map<vertex_descriptor, size_t> vertex_index_t;   

      //compute a spanning tree of maximum fractionality
      weight_map_t tc;
      vertex_descriptor v, s, t;
      edge_descriptor e;
      BGL_FORALL_EDGES_T(e, G, Graph){
         //negative weights because we need a maximum spanning tree
         tc[e] = - fabs(S.value(VM[e]) - .5);
      }
      vertex_index_t vertex_index_map;
      int i = 0;
      BGL_FORALL_VERTICES_T(v, G, Graph)
         vertex_index_map[v] = i++;

      //Property-maps required by MST
      associative_property_map<weight_map_t> pm_weight(tc);
      associative_property_map<vertex_index_t> pm_vertex_index(vertex_index_map);

      list<edge_descriptor> tree;


      vector < vertex_descriptor> p(num_vertices(G));

      kruskal_minimum_spanning_tree(G, back_inserter(tree), boost::weight_map
            (pm_weight).vertex_index_map(pm_vertex_index));

      if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
         cerr << "CUT::forestCycleSeparation: MIN_SPANNING_TREE computed\n";

      Graph T;
      map<vertex_descriptor, vertex_descriptor> TtoG, GtoT;
      BGL_FORALL_VERTICES_T(v, G, Graph) {
         t = add_vertex(T);
         TtoG[t] = v;
         GtoT[v] = t;
      }
      map<edge_descriptor, edge_descriptor> edgeTtoG;
      BOOST_FOREACH(e, tree)
         edgeTtoG[(add_edge(GtoT[source(e, G)], GtoT[target(e, G)], T)).first] = e;

      if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
         cerr << "CUT::forestCycleSeparation: Graph Copy T created\n";

      map<vertex_descriptor, int> dst;
      map<vertex_descriptor, edge_descriptor> pred;
      //FIXME choose start node
      undirected_bfs(T, *(vertices(T).first), dst, pred);

      if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
         cerr << "CUT::forestCycleSeparation: undirected_bfs done\n";

      // construct the fundamental cycles induced by the non-tree edges
      // and check whether a violated cycle inequality has been found
      BGL_FORALL_EDGES_T(e, G, Graph) {
         bool edgeInTree = false;
         edge_descriptor f;
         BOOST_FOREACH(f, tree)
            if (f == e) edgeInTree = true;
         if( edgeInTree )
            continue;
         if( numfound >= maxnum )
            break;
         s = GtoT[source(e, G)];
         t = GtoT[target(e, G)];
         int offset = dst[s] - dst[t];
         v = s;
         if( dst[s] - dst[t] < 0 ) {
            offset *= -1;
            v = t;
         }
         if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
            cerr << "CUT::forestCycleSeparation: offset computed\n";
         if( offset > 5 )
            continue;
         row r;
         double rhs = -1;
         if( S.value(VM[e]) > .5 ) {
            r += VM[e];
            rhs += 1.;
         }
         else
            r -= VM[e];
         if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
            cerr << "CUT::forestCycleSeparation: initial edge added\n";

         for( int i = 0; i < offset; i++ ) {
            if( S.value(VM[edgeTtoG[pred[v]]]) > .5 ) {
               r += VM[edgeTtoG[pred[v]]];
               rhs += 1.;
            }
            else 
               r -= VM[edgeTtoG[pred[v]]];
            (v == source(pred[v], T)) ? v = target(pred[v], T) 
               : v = source(pred[v], T);
         }
         if( dst[s] - dst[t] < 0 )
            t = v;
         else
            s = v;

         while( s != t ) {
            if( S.value(VM[edgeTtoG[pred[s]]]) > .5 ) {
               r += VM[edgeTtoG[pred[s]]];
               rhs += 1.;
            }
            else 
               r -= VM[edgeTtoG[pred[s]]];

            if( S.value(VM[edgeTtoG[pred[t]]]) > .5 ) {
               r += VM[edgeTtoG[pred[t]]];
               rhs += 1.;
            }
            else 
               r -= VM[edgeTtoG[pred[t]]];
            (s == source(pred[s], T)) ? s = target(pred[s], T) 
               : s = source(pred[s], T);
            (t == source(pred[t], T)) ? t = target(pred[t], T) 
               : t = source(pred[t], T);
         }

         r.normalize();
         //FIXME length
         if( !(int(rhs) % 2) ) { //&& r.size() < 20 ) 
            cons_obj* c = (r<=rhs);
            if( c->violation(S) > feps ) {
               newCons.push_back(r<=rhs);
               numfound++;
            }
            delete c;
         }
      }
      it2++;//connected components var_map iterator
   }
   if( S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "fcs: " << numfound << endl;
   return numfound;
}

template <typename Graph>
typename CUT<Graph>::status CUT<Graph>::feasible(solution& S) {
   if( S.has_os() && S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::feasible\n";

   double leps = .0000001;
   map<vertex_descriptor, int> dist, col;
   vertex_descriptor s;

   BGL_FORALL_VERTICES_T (s, G, Graph) {
      dist[s] = -1;
      col[s] = 0;
   }

   BGL_FORALL_VERTICES_T (s, G, Graph)
      if( dist[s] == -1 ) {
         list<vertex_descriptor> bfs_order;
         bfs_dist_visitor vis(&dist, &bfs_order, s);
         breadth_first_search(G, s, visitor(vis));
         vertex_descriptor v, w;
         col[bfs_order.front()] = 1;
         BOOST_FOREACH(v, bfs_order) {
            typename GraphTraits::out_edge_iterator oe_it, oe_it_end;
            /* Since G is undirected out_edges() returns iterator range to
               all adjacent edges.
             */
            for (tie(oe_it, oe_it_end) = out_edges(v, G);
                 oe_it != oe_it_end; oe_it++) {
               w = (source(*oe_it, G) == v) ? target(*oe_it, G) 
                                            : source(*oe_it, G);
               if(col[w] == 0) { //node w is not colored
                  col[w] = col[v];
                  if(S.value(VM[*oe_it]) > leps)
                     col[w] *= -1;
               } else { //node has already been colored
                  if(!(((col[w] == col[v]) && (S.value(VM[*oe_it]) <= leps))
                       || ((col[w] != col[v]) && (S.value(VM[*oe_it]) >= 1.-leps)))) {
                     if(S.originating_subproblem()->configuration("CUT_Debug_Out") == "true")
                        cerr << "feasible: Solution is not a feasible cut.\n";
                     return infeasible_solution;
                  }
               }
            }
         }
      }

   if( S.has_os() && S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "feasible: Solution is a feasible cut.\n";
   if( updateConfiguration ) {
      conf.clear();
      tconf.clear();
      BGL_FORALL_VERTICES_T(s, G, Graph) 
         (col[s] == 1) ? conf.push_back(s) : tconf.push_back(s);
      conf.sort();
      tconf.sort();
   }
   return feasible_solution;
}

template <typename Graph>
typename CUT<Graph>::status CUT<Graph>::standard_separation(subproblem& S) {
   if( S.configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::standard_separation\n";
   int numf = 1000;
   int nume = 100;

   list<cons_obj*> newCons;
   list<cons_obj*>::const_iterator ci;
   solution sol;
   cons_obj* c;
   cons nc;

   sol.save_solution(S);
   if( forestCycleSeparation(sol, numf, newCons) ) {
      foreach(ci, newCons)
         nc = S.add_basic_constraint(*ci);
      return constraint_found;
   }
   else if( exactCycleSeparation(sol, nume, newCons) ) {
      foreach(ci, newCons)
         nc = S.add_basic_constraint(*ci);
      return constraint_found;
   }
   else
      return no_constraint_found;
}

template <typename Graph>
typename CUT<Graph>::status CUT<Graph>::fast_separation(subproblem& S) {
   if( S.configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::fast_separation\n";
   int numf = 1000;

   list<cons_obj*> newCons;
   list<cons_obj*>::const_iterator ci;
   solution sol;
   cons_obj* c;
   cons nc;

   sol.save_solution(S);
   if( forestCycleSeparation(sol, numf, newCons) ) {
      foreach(ci, newCons)
         nc = S.add_basic_constraint(*ci);
      return constraint_found;
   } else
      return no_constraint_found;
}


template <typename Graph>
list<typename CUT<Graph>::vertex_descriptor> CUT<Graph>::getConfiguration(solution& S) {
   if( S.has_os() && S.originating_subproblem()->configuration("CUT_Debug_Out")=="true" )
      cerr << "CUT::getConfiguration\n";
   updateConfiguration = true;
   feasible(S);
   updateConfiguration = false;
   return (conf.size() < tconf.size()) ? conf : tconf;
}

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

---