SteinerArborescence.cc

#ifndef SCIL_STEINER_ARBORESCENCE_CC
#define SCIL_STEINER_ARBORESCENCE_CC

#define LEPS 1e-4

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

template <typename Graph>
SteinerArborescence<Graph>::SteinerArborescence(Graph& G_,
   vertex_descriptor root_, std::list<vertex_descriptor>& terminals_, 
   var_map<edge_descriptor>& VM_)
   : G(G_), root(root_), terminals(terminals_), VM(VM_) {
   feps = 1.0e-4;
   BGL_FORALL_EDGES_T(e, G, Graph)
      if(VM[e].type() != Vartype_Binary) {
         cerr << "SteinerArborescence.cc: All variables in VM";
         cerr << " need to be Vartype_Binary" << endl;
         exit(1);
      }
   //Graph has to be bidirectional
   BOOST_CONCEPT_ASSERT((boost::BidirectionalGraphConcept<Graph>));
   BOOST_CONCEPT_ASSERT((boost::VertexAndEdgeListGraphConcept<Graph>));
}

template <typename Graph>
void SteinerArborescence<Graph>::init(subproblem& S) {
   if (S.configuration("SteinerArborescence_Debug_Out")=="true")
      cerr << "SteinerArborescence::init\n";
   //the root node must not be a terminal
   typename std::list<vertex_descriptor>::iterator it;
   it = search_n(terminals.begin(), terminals.end(), 1, root);
   if( it != terminals.end() ) {
      cerr << "SteinerArborescence: The root node must not be a terminal!" << endl;
      exit(1);
   }

   //all terminals have indegree 1
   //all Steiner nodes have indegree <=1
   BGL_FORALL_VERTICES_T(v, G, Graph) {
      row r;
      typename GraphTraits::in_edge_iterator ie_it, ie_it_end;
      for(tie(ie_it, ie_it_end) = in_edges(v,G); ie_it != ie_it_end; ie_it++)
         r += VM[*ie_it];
      it = search_n(terminals.begin(), terminals.end(), 1, v);
      if( it != terminals.end() )
         S.add_basic_constraint(r==1., Static);
      else if( v != root )
         S.add_basic_constraint(r<=1., Static);
   }
   row r;
   typename GraphTraits::out_edge_iterator oe_it, oe_it_end;
   for(tie(oe_it, oe_it_end) = out_edges(root, G); oe_it != oe_it_end; oe_it++)
      r += VM[*oe_it];
   S.add_basic_constraint(r>=1.,Static);

   int i = 0;
   BGL_FORALL_VERTICES_T(v, G, Graph) {
      vertex_descriptor new_vertex_Gmf = add_vertex(Gmf);
      vertex_descriptor new_vertex_Fmf = add_vertex(Fmf);
      vertices_GtoGmf[v] = new_vertex_Gmf;
      vertices_GtoFmf[v] = new_vertex_Fmf;
      vertex_index_Gmf[new_vertex_Gmf] = i;
      vertex_index_Fmf[new_vertex_Fmf] = i;
      i++;
   }
   //Maps each original edge (u,v) in Gmf to its reverse edge (v,u)
   BGL_FORALL_EDGES_T(e, G, Graph) {
      edge_descriptor original_edge_Gmf, original_edge_Fmf, reverse_edge_Gmf, reverse_edge_Fmf;
      vertex_descriptor s = vertices_GtoGmf[source(e, G)];
      vertex_descriptor t = vertices_GtoGmf[target(e, G)];
      vertex_descriptor u = vertices_GtoFmf[source(e, G)];
      vertex_descriptor v = vertices_GtoFmf[target(e, G)];
      original_edge_Gmf = add_edge(s, t, Gmf).first;
      original_edge_Fmf = add_edge(v, u, Fmf).first;
      is_original_edge_Gmf[original_edge_Gmf] = true;
      is_original_edge_Fmf[original_edge_Fmf] = true;
      if(!edge(target(e, G), source(e, G), G).second){
         reverse_edge_Gmf = (add_edge(t, s, Gmf)).first;
         is_original_edge_Gmf[reverse_edge_Gmf] = false;
         reverse_edge_Fmf = (add_edge(u, v, Fmf)).first;
         is_original_edge_Fmf[reverse_edge_Fmf] = false;
      }
      else if(edge(t, s, Gmf).second){
         reverse_edge_Gmf = edge(t, s, Gmf).first;
         rev_Gmf[original_edge_Gmf] = reverse_edge_Gmf;
         rev_Gmf[reverse_edge_Gmf] = original_edge_Gmf;
         reverse_edge_Fmf = edge(u, v, Fmf).first;
         rev_Fmf[original_edge_Fmf] = reverse_edge_Fmf;
         rev_Fmf[reverse_edge_Fmf] = original_edge_Fmf;
      } 
      edges_GmfToG[original_edge_Gmf] = e;
      edges_FmfToG[original_edge_Fmf] = e;
   }

   return;
};

/* This visitors sets the boolean values of all vertices in map reached to true
   if a vertex is visited
 */
template <typename Graph>
class SteinerArborescence<Graph>::dfs_visitor : public default_dfs_visitor {
   public:
      dfs_visitor(map<vertex_descriptor, bool>* _reached) : reached(_reached) {}
      void discover_vertex(vertex_descriptor u, const Graph &G) const {
         (*reached)[u] = true;
      }
      map<vertex_descriptor, bool> *reached;
};

template <typename Graph>
typename SteinerArborescence<Graph>::status 
SteinerArborescence<Graph>::feasible(solution& S) {
   if (S.originating_subproblem()->configuration("SteinerArborescence_Debug_Out")=="true")
      cerr << "SteinerArborescence::feasible\n";
   map<vertex_descriptor, vertex_descriptor> GtoH;
   Graph H;
   BGL_FORALL_VERTICES_T(u1, G, Graph)
      GtoH[u1] = 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);

   //Input for depth_first_visit
   map<vertex_descriptor, bool> reached;
   typedef map<vertex_descriptor, default_color_type> color_map_t;
   color_map_t color_map;
   BGL_FORALL_VERTICES_T(v, H, Graph) {
      reached[v] = false;
      color_map[v] = white_color;
   }
   associative_property_map<color_map_t> pm_color(color_map);
   dfs_visitor vis(&reached);
   //depth_first_visit visits just the connected component containing GtoH[root]
   depth_first_visit(H, GtoH[root], vis, pm_color);

   std::list<vertex_descriptor> reached_vertices;
   BGL_FORALL_VERTICES_T(v, H, Graph)
      if(reached[v]) reached_vertices.push_back(v);
   //FIXME

   vertex_descriptor u1;
   BOOST_FOREACH(u1, terminals) {
   //  if( reached.search(GtoH[u1]) == nil )
   //   return infeasible_solution;
      if ( !reached[GtoH[u1]] ) return infeasible_solution;
   }
   return feasible_solution;
}

template <typename Graph>
int SteinerArborescence<Graph>::separate_flipped(subproblem& S) {
   if (S.configuration("SteinerArborescence_Debug_Out")=="true")
      cerr << "SteinerArborescence::seperate_flipped\n";
   int numfound = 0;
   capacity_map_t c;
   capacity_map_t d;
   capacity_map_t residual_capacity;
   BGL_FORALL_EDGES_T(e, Fmf, Graph)
      residual_capacity[e] = 0;
   BGL_FORALL_EDGES_T(e, Fmf, Graph) {
      if(is_original_edge_Fmf[e]) {
         //edge is an original edge
         c[e] = (int) (1000. * S.value(VM[edges_FmfToG[e]]));
         if( c[e] == 0 )
            c[e] = 1;
      } else {
         //edge is a reverse edge
         c[e] = 0;
      }
   }
   //ColorMap for boykov_kolmogorov_max_flow
   color_map_t color_map;
   BGL_FORALL_VERTICES_T(v, Fmf, Graph)
      color_map[v] = white_color;
   //Predecessor map for boykov_kolmogorov_max_flow
   predecessor_map_t pred_map;
   //Distance map for boykov_kolmogorov_max_flow
   vdistance_map_t vdist_map;


   //Property maps for boykov_kolmogorov_max_flow
   associative_property_map<capacity_map_t> pm_capacity(d);
   associative_property_map<capacity_map_t> pm_res_cap(residual_capacity);
   associative_property_map<rev_edge_map_t> pm_rev(rev_Fmf);
   associative_property_map<vertex_index_map_t> pm_vi(vertex_index_Fmf);
   associative_property_map<color_map_t> pm_color(color_map);
   associative_property_map<predecessor_map_t> pm_pred(pred_map);
   associative_property_map<vdistance_map_t> pm_vdist(vdist_map);

   vertex_descriptor v;
   BOOST_FOREACH(v, terminals) {
      std::list<vertex_descriptor> cut;
      BGL_FORALL_EDGES_T(e, Fmf, Graph) 
         d[e] = c[e];
      map<edge_descriptor, bool> fixed;
      BGL_FORALL_EDGES_T(e, Fmf, Graph)
         fixed[e] = false;
      int cv = 0;
      int go = 0;
      bool added = true;
      while( added ) {
         BGL_FORALL_VERTICES_T(u, Fmf, Graph)
            color_map[u] = white_color;
         cv = boykov_kolmogorov_max_flow(Fmf, pm_capacity, pm_res_cap, pm_rev,
                 pm_pred, pm_color, pm_vdist, pm_vi, vertices_GtoFmf[v],
                 vertices_GtoFmf[root]);
         go++;
         added = false;
         cut.clear();
         typename GraphTraits::vertex_iterator v_it, v_it_end;
         /* Vertices that are colored black by boykov_kolmogorov_max_flow are in
            the source tree and form a minimum s-t-cut
          */
         for(tie(v_it,v_it_end) = vertices(Fmf); v_it != v_it_end; v_it++)
            if (color_map[*v_it] == black_color) cut.push_back(*v_it);
         if (S.configuration("SteinerArborescence_Debug_Out")=="true")
            cerr << "cv: " << cv << endl;
         if( cv < 1000 && cut.size() >= 1 && cut.size() < num_vertices(Fmf) - 1 ) {
            row r;
            vertex_descriptor w;
            BOOST_FOREACH(w, cut) {
               BGL_FORALL_OUTEDGES_T(w, e, Fmf, Graph) {
                  if (search_n(cut.begin(),cut.end(),1,target(e,Fmf)) == cut.end()
                      && is_original_edge_Fmf[e] /* && !fixed[e] */ ) {
                     r += VM[edges_FmfToG[e]];
                     d[e] = 1000;
                     fixed[e] = true;
                  }
               }
            }
            cons_obj* con = (r>=1.);
            if( con->violation(S) > feps ) {
               S.add_basic_constraint(r>=1., Dynamic);
               if (S.configuration("SteinerArborescence_Debug_Out")=="true")
                  if( go > 1 )
                     cerr << go << "  added: " << r << " >= 1" << endl;
               numfound++;
               added = true;
            }
            delete con;
         }
      }
   }
 


}

template <typename Graph>
typename SteinerArborescence<Graph>::status
SteinerArborescence<Graph>::standard_separation(subproblem& S){
   return separation(S, false);
}

template <typename Graph>
typename SteinerArborescence<Graph>::status 
SteinerArborescence<Graph>::separation(subproblem& S,
   bool separate_fast) {
   if (S.configuration("SteinerArborescence_Debug_Out")=="true")
      cerr << "SteinerArborescence::standard_separation\n";

   int numfound = 0;
   capacity_map_t c;
   capacity_map_t d;
   capacity_map_t residual_capacity;
   BGL_FORALL_EDGES_T(e, Gmf, Graph)
      residual_capacity[e] = 0;
   BGL_FORALL_EDGES_T(e, Gmf, Graph) {
      if(is_original_edge_Gmf[e]) {
         //edge is an original edge
         c[e] = (int) (1000. * S.value(VM[edges_GmfToG[e]]));
         if( c[e] == 0 )
            c[e] = 1;
      } else {
         //edge is a reverse edge
         c[e] = 0;
      }
   }
   //ColorMap for boykov_kolmogorov_max_flow
   color_map_t color_map;
   BGL_FORALL_VERTICES_T(v, Gmf, Graph)
      color_map[v] = white_color;
   //Predecessor map for boykov_kolmogorov_max_flow
   predecessor_map_t pred_map;
   //Distance map for boykov_kolmogorov_max_flow
   vdistance_map_t vdist_map;


   //Property maps for boykov_kolmogorov_max_flow
   associative_property_map<capacity_map_t> pm_capacity(d);
   associative_property_map<capacity_map_t> pm_res_cap(residual_capacity);
   associative_property_map<rev_edge_map_t> pm_rev(rev_Gmf);
   associative_property_map<vertex_index_map_t> pm_vi(vertex_index_Gmf);
   associative_property_map<color_map_t> pm_color(color_map);
   associative_property_map<predecessor_map_t> pm_pred(pred_map);
   associative_property_map<vdistance_map_t> pm_vdist(vdist_map);

   vertex_descriptor v;
   BOOST_FOREACH(v, terminals) {
      std::list<vertex_descriptor> cut;
      BGL_FORALL_EDGES_T(e, Gmf, Graph) 
         d[e] = c[e];
      map<edge_descriptor, bool> fixed;
      BGL_FORALL_EDGES_T(e, Gmf, Graph)
         fixed[e] = false;
      int cv = 0;
      int go = 0;
      bool added = true;
      while( added ) {
         BGL_FORALL_VERTICES_T(u, Gmf, Graph)
            color_map[u] = white_color;
         cv = boykov_kolmogorov_max_flow(Gmf, pm_capacity, pm_res_cap, pm_rev,
                 pm_pred, pm_color, pm_vdist, pm_vi, vertices_GtoGmf[root],
                 vertices_GtoGmf[v]);
         go++;
         added = false;
         cut.clear();
         typename GraphTraits::vertex_iterator v_it, v_it_end;
         /* Vertices that are colored black by boykov_kolmogorov_max_flow are in
            the source tree and form a minimum s-t-cut
          */
         for(tie(v_it,v_it_end) = vertices(Gmf); v_it != v_it_end; v_it++)
            if (color_map[*v_it] == black_color) cut.push_back(*v_it);
         if (S.configuration("SteinerArborescence_Debug_Out")=="true")
            cerr << "cv: " << cv << endl;
         if( cv < 1000 && cut.size() >= 1 && cut.size() < num_vertices(Gmf) - 1 ) {
            row r;
            vertex_descriptor w;
            BOOST_FOREACH(w, cut) {
               BGL_FORALL_OUTEDGES_T(w, e, Gmf, Graph) {
                  if (search_n(cut.begin(),cut.end(),1,target(e,Gmf)) == cut.end()
                      && is_original_edge_Gmf[e] /* && !fixed[e] */ ) {
                     r += VM[edges_GmfToG[e]];
                     d[e] = 1000;
                     fixed[e] = true;
                  }
               }
            }
            cons_obj* con = (r>=1.);
            if( con->violation(S) > feps ) {
               S.add_basic_constraint(r>=1., Dynamic);
               if (S.configuration("SteinerArborescence_Debug_Out")=="true")
                  if( go > 1 )
                     cerr << go << "  added: " << r << " >= 1" << endl;
               numfound++;
               added = true;
            }
            delete con;
         }
      }
   }
   //Stop separating if fast separation was chosen
   if (separate_fast) {
      if ( numfound )
         return constraint_found;
      return no_constraint_found;
   }
   numfound += separate_flipped(S);
   if (S.configuration("SteinerArborescence_Debug_Out")=="true")
      cerr << "numfound: " << numfound << endl;

   //cerr << "separating type 4 ineqs" << endl;
   std::map<vertex_descriptor, double> nodecut;
   std::map<vertex_descriptor, row> noderow;
   BGL_FORALL_VERTICES_T(v, Gmf, Graph) {
      if( search_n(terminals.begin(), terminals.end(), 1, vertices_GtoGmf[v])
          != terminals.end() || vertices_GtoGmf[root] == v )
         continue;
      double t = 0;
      row r;
      BGL_FORALL_INEDGES_T(v, e, Gmf, Graph){
         if(is_original_edge_Gmf[e]){
            t += c[e];
            r += VM[edges_GmfToG[e]];
         }
      }
      nodecut[v] = t;
      noderow[v] = r;
      //cerr << nodecut[v] << " ";
   }
   //cerr << endl;

   BGL_FORALL_VERTICES_T(v, Gmf, Graph) {
      if(search_n(terminals.begin(), terminals.end(), 1, v) != terminals.end()
         || v == vertices_GtoGmf[root] )
         continue;
      std::list<vertex_descriptor> cut, cut2;
      bool stop = false;

      //Initialize data for boykov_kolmogorov_max_flow that was changed
      associative_property_map<capacity_map_t> pm_cap(c);
      BGL_FORALL_EDGES_T(e, Gmf, Graph)
         residual_capacity[e] = 0;
      BGL_FORALL_VERTICES_T(v, Gmf, Graph)
         color_map[v] = white_color;

      int cv = 0;
      cv = boykov_kolmogorov_max_flow(Gmf, pm_cap, pm_res_cap, pm_rev, pm_pred,
              pm_color, pm_vdist, pm_vi, vertices_GtoGmf[root], v);

      cut.clear();
      typename GraphTraits::vertex_iterator v_it, v_it_end;
      /* Vertices that are colored black by boykov_kolmogorov_max_flow are in
         the source tree and form a minimum s-t-cut
       */
      for(tie(v_it,v_it_end) = vertices(Gmf); v_it != v_it_end; v_it++){
         if (color_map[*v_it] == black_color)
            cut.push_back(*v_it);
         else
            cut2.push_back(*v_it);
      }
      //cerr << "cv: " << cv << endl;

      if( cut.size() >= 1 && cut.size() < num_vertices(Gmf) - 1 ) {
         vertex_descriptor w;
         BOOST_FOREACH(w, terminals) {
            if( search_n(cut.begin(),cut.end(),1,vertices_GtoGmf[w])==cut.end() ) {
               stop = true;
               //cerr << "terminal not in cut" << endl;
               break;
            }
         }
         if( stop )
            break;

         row r;
         BOOST_FOREACH(w, cut) {
            BGL_FORALL_OUTEDGES_T(w, e, Gmf, Graph) {
               if(search_n(cut.begin(),cut.end(),1,target(e, Gmf)) == cut.end()
                  && is_original_edge_Gmf[e])
                  r += VM[edges_GmfToG[e]];
            }
         }
         BOOST_FOREACH(w, cut2) {
            if( cv - nodecut[w] < 20. ) {
               cons_obj* con = (r-noderow[w]>=0.);
               if( con->violation(S) > feps ) {
                  S.add_basic_constraint(r-noderow[w]>=0., Dynamic);
                  //cerr << "added: " << r-noderow[w] << " >= 0" << endl;
                  numfound++;
               }
               delete con;
            }
         }
      }
   }

   if( numfound )
      return constraint_found;

   return no_constraint_found;
}

template <typename Graph>
typename SteinerArborescence<Graph>::status 
SteinerArborescence<Graph>::fast_separation(subproblem& S) {
   return separation(S, true);
}

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

---