spantree.cc

#ifndef SPAN_TREE_CC
#define SPAN_TREE_CC

#include<scil/scil.h>
#include<scil/global.h>
#include<boost/graph/iteration_macros.hpp>
#include<boost/foreach.hpp>
#include<boost/graph/boykov_kolmogorov_max_flow.hpp>
#include<boost/graph/connected_components.hpp>


using namespace SCIL;
using namespace std;
using namespace boost;

using boost::edge_capacity;
using boost::edge_flow;

#define VALONE 10000

template
<typename Graph>
class SpanTree<Graph>::st_cutset_inequality : public cons_obj {

   private: 
      typedef typename graph_traits<Graph>::edge_descriptor edge;
      typedef typename graph_traits<Graph>::vertex_descriptor vertex;

   private:
      Graph& G; 
      var_map<edge>& VM;
      std::map<vertex,bool> cut;
      double rhs;
  
   public:
  
      st_cutset_inequality(Graph& Gt, var_map<edge>& VM_, std::list<vertex>& L) 
         : cons_obj(Greater, 1), G(Gt), VM(VM_) {
   
         vertex v;
         BGL_FORALL_VERTICES_T(v,G,Graph)
            cut[v] = false;
 
         BOOST_FOREACH(v,L)
            cut[v] = true;
      }

      void non_zero_entries(row& r) { 
         BGL_FORALL_EDGES_T(e,G,Graph) {
            if(coeff(e)==1) r+=VM[e];
         }
      }

   virtual
      ~st_cutset_inequality() {};

   double coeff(edge e) { 
      if(cut[source(e,G)] != cut[target(e,G)]) return 1;
      return(0);
   };
};


template
<typename Graph>
SpanTree<Graph>::SpanTree(Graph& G_, var_map<edge_descriptor>& VM_)
   : G(G_), VM(VM_)
{ 
  // TODO: some more structure-requirements???
   /*if( is_undirected(G) == false ){
      cerr << "spantree.cc: The Graph has to be ";
      cerr << "undirected" << endl;
      exit(1);
   }*/

   edge_descriptor e;
   BGL_FORALL_EDGES_T(e,G,Graph) 
   if( VM[e].type() != Vartype_Binary ) {
      cerr << "spantree.cc: All variables in VM need to be ";
      cerr << "Vartype_Binary" << endl;
      exit(1);
   }
};

template
<typename Graph>
void SpanTree<Graph>::init(subproblem& S) {
   if (S.configuration("SpanTree_Debug_Out")=="true")
      cerr << "SpanTree::init\n";
   row r;
   H_edge e1,e2;   
   H_edge r1,r2;

   BGL_FORALL_EDGES_T(e,G,Graph) { 
      r+=VM[e]; 
   }
   S.add_basic_constraint(r == num_vertices(G)-1, Static);

   BGL_FORALL_VERTICES_T(v,G,Graph){
      GtoH[v] = add_vertex(H);
   }

   // add edges e, inv(e) to H 
   property_map<bdGraph, edge_reverse_t>::type 
      rev = get(edge_reverse, H);   
   BGL_FORALL_EDGES_T(e,G,Graph){
      GtoH1[e] = add_edge( GtoH[source(e,G)], GtoH[target(e,G)], H ).first;   
      GtoH2[e] = add_edge( GtoH[target(e,G)], GtoH[source(e,G)], H ).first;
      rev[GtoH1[e]] = GtoH2[e];
      rev[GtoH2[e]] = GtoH1[e];
   }

   son = add_vertex(H);
   tan = add_vertex(H);

   BGL_FORALL_VERTICES_T(v,G,Graph) {
      e1 = add_edge(son, GtoH[v], H).first;
      e2 = add_edge(GtoH[v], tan, H).first;

      // add reverse edges too:
      r1 = add_edge(GtoH[v], son, H).first;
      r2 = add_edge(tan, GtoH[v], H).first;
      rev[e1] = r1;
      rev[r1] = e1;
      rev[e2] = r2;
      rev[r2] = e2;
   }

}

     
template
<typename Graph>
void SpanTree<Graph>::min_cut(H_vertex u, property_map<bdGraph, edge_capacity_t>::type& C,  
   property_map<bdGraph, edge_residual_capacity_t>::type& resC,
   property_map<bdGraph, vertex_reached_t>::type& R, int& k) {

   R[u]=true;
   k++;

   BGL_FORALL_OUTEDGES_T(u,e,H,bdGraph)
      if( (resC[e]>0) && (!R[target(e,H)]) )     // resC != 0  <=>  flow != max
         min_cut(target(e,H), C, resC, R, k);

   BGL_FORALL_INEDGES_T(u,e,H,bdGraph)
      if((resC[e]<C[e]) && (!R[source(e,H)])) {  // resC != C  <=>  flow != 0
         min_cut(source(e,H), C, resC, R, k);
   };
};


template
<typename Graph>
typename SpanTree<Graph>::status SpanTree<Graph>::standard_separation(subproblem& S) {
   if (S.configuration("SpanTree_Debug_Out")=="true")
      cerr << "SpanTree::standard_separation\n";

   status sta=no_constraint_found;
   // Preparation: Calculate Capacity
   int j, k;

   std::vector<int> Con(num_vertices(H));

   property_map<bdGraph, edge_capacity_t>::type Cap
      = get(edge_capacity, H);

   BGL_FORALL_EDGES_T(e, G, Graph) { 
      j=(int) (VALONE*S.value(VM[e]));
      if(j<0) j=0;
      Cap[GtoH1[e]]=j;
      Cap[GtoH2[e]]=j;
      Con[GtoH[source(e,G)]]+=j;
      Con[GtoH[target(e,G)]]+=j;
   };
  
   BGL_FORALL_OUTEDGES_T(son,he,H,bdGraph) {
      Cap[he]=std::max(Con[target(he,H)]-2*VALONE,0);
   }

   BGL_FORALL_INEDGES_T(tan,he,H,bdGraph) {
      Cap[he]=std::max(2*VALONE-Con[source(he,H)],0);
   }
  
   if(sta==constraint_found) return sta;                   

   row q = 0;
   bool cancel;
   property_map<bdGraph, vertex_reached_t>::type        
      R = get(vertex_reached_t(),H);

   property_map<bdGraph, edge_residual_capacity_t>::type
      resCap = get(edge_residual_capacity, H);

   BGL_FORALL_OUTEDGES_T(son, he, H, bdGraph) {
      k=Cap[he];
      Cap[he]=1000*VALONE;

      j = boykov_kolmogorov_max_flow(H, son, tan);
      // Make Constraint

      BGL_FORALL_VERTICES_T(hv,H,bdGraph) 
         R[hv]=false;       
      j=0;
      min_cut(son, Cap, resCap, R, j);
      row r;
      double d=0;
      BGL_FORALL_EDGES_T(f, G,Graph) {
         if((R[GtoH[source(f,G)]]) && (R[GtoH[target(f,G)]])) {
            r+=VM[f];
            d+=S.value(VM[f]);
         };
      }
      if(d>((double) j)-1.99) {
         cancel = false;
         if(q.size()!=0){
            q-=r;
            q.normalize(true);
            if(q.size()==0)
               cancel = true;       
         }
         q = r;
         if(!cancel){
           S.add_basic_constraint(r<=j-2, Dynamic);
           sta=constraint_found;
         }
      };
      Cap[he]=k;
   };
   return sta; 
}

template
<typename Graph>
typename SpanTree<Graph>::status SpanTree<Graph>::feasible(solution& S) {
   if (S.has_os() && S.originating_subproblem()->configuration("SpanTree_Debug_Out")=="true")
      cerr << "SpanTree::init\n";
   std::map<vertex,vertex> GtoH;
   Graph H;
  
   BGL_FORALL_VERTICES_T(v,G,Graph)
      GtoH[v] = add_vertex(H);

   BGL_FORALL_EDGES_T(e,G,Graph)
      if(S.value(VM[e])>0.5) {
         add_edge(GtoH[source(e,G)], GtoH[target(e,G)], H);
         add_edge(GtoH[target(e,G)], GtoH[source(e,G)], H);     // regard H as undirected (to be more precise: bidirected graph)
      }

   std::vector<int> component(num_vertices(H));
   int num = connected_components(H, &component[0]);

   if( num==1 ) return feasible_solution;
   if (S.has_os() && S.originating_subproblem()->configuration("SpanTree_Debug_Out")=="true")
      cerr << "integer feasible, but infeasible\n";
   return infeasible_solution;
}

template
<typename Graph>
void SpanTree<Graph>::info() {
   cout<<"Configurations:\n";
   cout<<"SpanTree_Debug_Out [true|false]\n";
};


#undef VALONE

#endif

---