File indexing completed on 2024-12-22 04:10:11

0001 //  SPDX-FileCopyrightText: 2006 Stephan Diederich
0002 //
0003 //  This code may be used under either of the following two licences:
0004 //
0005 //  SPDX-License-Identifier: MIT OR BSL-1.0
0006 
0007 #ifndef BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
0008 #define BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
0009 
0010 #include <boost/config.hpp>
0011 #include <boost/assert.hpp>
0012 #include <vector>
0013 #include <list>
0014 #include <utility>
0015 #include <iosfwd>
0016 #include <algorithm> // for std::min and std::max
0017 
0018 #include <boost/pending/queue.hpp>
0019 #include <boost/limits.hpp>
0020 #include <boost/property_map/property_map.hpp>
0021 #include <boost/none_t.hpp>
0022 #include <boost/graph/graph_concepts.hpp>
0023 #include <boost/graph/named_function_params.hpp>
0024 #include <boost/graph/lookup_edge.hpp>
0025 #include <boost/concept/assert.hpp>
0026 
0027 // The algorithm implemented here is described in:
0028 //
0029 // Boykov, Y., Kolmogorov, V. "An Experimental Comparison of Min-Cut/Max-Flow
0030 // Algorithms for Energy Minimization in Vision", In IEEE Transactions on
0031 // Pattern Analysis and Machine Intelligence, vol. 26, no. 9, pp. 1124-1137,
0032 // Sep 2004.
0033 //
0034 // For further reading, also see:
0035 //
0036 // Kolmogorov, V. "Graph Based Algorithms for Scene Reconstruction from Two or
0037 // More Views". PhD thesis, Cornell University, Sep 2003.
0038 
0039 namespace boost {
0040 
0041 namespace detail {
0042 
0043 template <class Graph,
0044           class EdgeCapacityMap,
0045           class ResidualCapacityEdgeMap,
0046           class ReverseEdgeMap,
0047           class PredecessorMap,
0048           class ColorMap,
0049           class DistanceMap,
0050           class IndexMap>
0051 class bk_max_flow {
0052   typedef typename property_traits<EdgeCapacityMap>::value_type tEdgeVal;
0053   typedef graph_traits<Graph> tGraphTraits;
0054   typedef typename tGraphTraits::vertex_iterator vertex_iterator;
0055   typedef typename tGraphTraits::vertex_descriptor vertex_descriptor;
0056   typedef typename tGraphTraits::edge_descriptor edge_descriptor;
0057   typedef typename tGraphTraits::edge_iterator edge_iterator;
0058   typedef typename tGraphTraits::out_edge_iterator out_edge_iterator;
0059   typedef boost::queue<vertex_descriptor> tQueue;                               //queue of vertices, used in adoption-stage
0060   typedef typename property_traits<ColorMap>::value_type tColorValue;
0061   typedef color_traits<tColorValue> tColorTraits;
0062   typedef typename property_traits<DistanceMap>::value_type tDistanceVal;
0063 
0064     public:
0065       bk_max_flow(Graph& g,
0066                   EdgeCapacityMap cap,
0067                   ResidualCapacityEdgeMap res,
0068                   ReverseEdgeMap rev,
0069                   PredecessorMap pre,
0070                   ColorMap color,
0071                   DistanceMap dist,
0072                   IndexMap idx,
0073                   vertex_descriptor src,
0074                   vertex_descriptor sink):
0075       m_g(g),
0076       m_index_map(idx),
0077       m_cap_map(cap),
0078       m_res_cap_map(res),
0079       m_rev_edge_map(rev),
0080       m_pre_map(pre),
0081       m_tree_map(color),
0082       m_dist_map(dist),
0083       m_source(src),
0084       m_sink(sink),
0085       m_active_nodes(),
0086       m_in_active_list_vec(num_vertices(g), false),
0087       m_in_active_list_map(make_iterator_property_map(m_in_active_list_vec.begin(), m_index_map)),
0088       m_has_parent_vec(num_vertices(g), false),
0089       m_has_parent_map(make_iterator_property_map(m_has_parent_vec.begin(), m_index_map)),
0090       m_time_vec(num_vertices(g), 0),
0091       m_time_map(make_iterator_property_map(m_time_vec.begin(), m_index_map)),
0092       m_flow(0),
0093       m_time(1),
0094       m_last_grow_vertex(graph_traits<Graph>::null_vertex()){
0095         // initialize the color-map with gray-values
0096         vertex_iterator vi, v_end;
0097         for(boost::tie(vi, v_end) = vertices(m_g); vi != v_end; ++vi){
0098           set_tree(*vi, tColorTraits::gray());
0099         }
0100         // Initialize flow to zero which means initializing
0101         // the residual capacity equal to the capacity
0102         edge_iterator ei, e_end;
0103         for(boost::tie(ei, e_end) = edges(m_g); ei != e_end; ++ei) {
0104           put(m_res_cap_map, *ei, get(m_cap_map, *ei));
0105           BOOST_ASSERT(get(m_rev_edge_map, get(m_rev_edge_map, *ei)) == *ei); //check if the reverse edge map is build up properly
0106         }
0107         //init the search trees with the two terminals
0108         set_tree(m_source, tColorTraits::black());
0109         set_tree(m_sink, tColorTraits::white());
0110         put(m_time_map, m_source, 1);
0111         put(m_time_map, m_sink, 1);
0112       }
0113 
0114       tEdgeVal max_flow(){
0115         //augment direct paths from SOURCE->SINK and SOURCE->VERTEX->SINK
0116         augment_direct_paths();
0117         //start the main-loop
0118         while(true){
0119           bool path_found;
0120           edge_descriptor connecting_edge;
0121           boost::tie(connecting_edge, path_found) = grow(); //find a path from source to sink
0122           if(!path_found){
0123             //we're finished, no more paths were found
0124             break;
0125           }
0126           ++m_time;
0127           augment(connecting_edge); //augment that path
0128           adopt(); //rebuild search tree structure
0129         }
0130         return m_flow;
0131       }
0132 
0133       // the complete class is protected, as we want access to members in
0134       // derived test-class (see test/boykov_kolmogorov_max_flow_test.cpp)
0135     protected:
0136       void augment_direct_paths(){
0137         // in a first step, we augment all direct paths from source->NODE->sink
0138         // and additionally paths from source->sink. This improves especially
0139         // graphcuts for segmentation, as most of the nodes have source/sink
0140         // connects but shouldn't have an impact on other maxflow problems
0141         // (this is done in grow() anyway)
0142         out_edge_iterator ei, e_end;
0143         for(boost::tie(ei, e_end) = out_edges(m_source, m_g); ei != e_end; ++ei){
0144           edge_descriptor from_source = *ei;
0145           vertex_descriptor current_node = target(from_source, m_g);
0146           if(current_node == m_sink){
0147             tEdgeVal cap = get(m_res_cap_map, from_source);
0148             put(m_res_cap_map, from_source, 0);
0149             m_flow += cap;
0150             continue;
0151           }
0152           edge_descriptor to_sink;
0153           bool is_there;
0154           boost::tie(to_sink, is_there) = lookup_edge(current_node, m_sink, m_g);
0155           if(is_there){
0156             tEdgeVal cap_from_source = get(m_res_cap_map, from_source);
0157             tEdgeVal cap_to_sink = get(m_res_cap_map, to_sink);
0158             if(cap_from_source > cap_to_sink){
0159               set_tree(current_node, tColorTraits::black());
0160               add_active_node(current_node);
0161               set_edge_to_parent(current_node, from_source);
0162               put(m_dist_map, current_node, 1);
0163               put(m_time_map, current_node, 1);
0164               // add stuff to flow and update residuals. we don't need to
0165               // update reverse_edges, as incoming/outgoing edges to/from
0166               // source/sink don't count for max-flow
0167               put(m_res_cap_map, from_source, get(m_res_cap_map, from_source) - cap_to_sink);
0168               put(m_res_cap_map, to_sink, 0);
0169               m_flow += cap_to_sink;
0170             } else if(cap_to_sink > 0){
0171               set_tree(current_node, tColorTraits::white());
0172               add_active_node(current_node);
0173               set_edge_to_parent(current_node, to_sink);
0174               put(m_dist_map, current_node, 1);
0175               put(m_time_map, current_node, 1);
0176               // add stuff to flow and update residuals. we don't need to update
0177               // reverse_edges, as incoming/outgoing edges to/from source/sink
0178               // don't count for max-flow
0179               put(m_res_cap_map, to_sink, get(m_res_cap_map, to_sink) - cap_from_source);
0180               put(m_res_cap_map, from_source, 0);
0181               m_flow += cap_from_source;
0182             }
0183           } else if(get(m_res_cap_map, from_source)){
0184             // there is no sink connect, so we can't augment this path, but to
0185             // avoid adding m_source to the active nodes, we just activate this
0186             // node and set the appropriate things
0187             set_tree(current_node, tColorTraits::black());
0188             set_edge_to_parent(current_node, from_source);
0189             put(m_dist_map, current_node, 1);
0190             put(m_time_map, current_node, 1);
0191             add_active_node(current_node);
0192           }
0193         }
0194         for(boost::tie(ei, e_end) = out_edges(m_sink, m_g); ei != e_end; ++ei){
0195           edge_descriptor to_sink = get(m_rev_edge_map, *ei);
0196           vertex_descriptor current_node = source(to_sink, m_g);
0197           if(get(m_res_cap_map, to_sink)){
0198             set_tree(current_node, tColorTraits::white());
0199             set_edge_to_parent(current_node, to_sink);
0200             put(m_dist_map, current_node, 1);
0201             put(m_time_map, current_node, 1);
0202             add_active_node(current_node);
0203           }
0204         }
0205       }
0206 
0207       /**
0208        * Returns a pair of an edge and a boolean. if the bool is true, the
0209        * edge is a connection of a found path from s->t , read "the link" and
0210        * source(returnVal, m_g) is the end of the path found in the source-tree
0211        * target(returnVal, m_g) is the beginning of the path found in the sink-tree
0212        */
0213       std::pair<edge_descriptor, bool> grow(){
0214         BOOST_ASSERT(m_orphans.empty());
0215         vertex_descriptor current_node;
0216         while((current_node = get_next_active_node()) != graph_traits<Graph>::null_vertex()){ //if there is one
0217           BOOST_ASSERT(get_tree(current_node) != tColorTraits::gray() &&
0218                        (has_parent(current_node) ||
0219                          current_node == m_source ||
0220                          current_node == m_sink));
0221 
0222           if(get_tree(current_node) == tColorTraits::black()){
0223             //source tree growing
0224             out_edge_iterator ei, e_end;
0225             if(current_node != m_last_grow_vertex){
0226               m_last_grow_vertex = current_node;
0227               boost::tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g);
0228             }
0229             for(; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it) {
0230               edge_descriptor out_edge = *m_last_grow_edge_it;
0231               if(get(m_res_cap_map, out_edge) > 0){ //check if we have capacity left on this edge
0232                 vertex_descriptor other_node = target(out_edge, m_g);
0233                 if(get_tree(other_node) == tColorTraits::gray()){ //it's a free node
0234                   set_tree(other_node, tColorTraits::black()); //acquire other node to our search tree
0235                   set_edge_to_parent(other_node, out_edge);   //set us as parent
0236                   put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);  //and update the distance-heuristic
0237                   put(m_time_map, other_node, get(m_time_map, current_node));
0238                   add_active_node(other_node);
0239                 } else if(get_tree(other_node) == tColorTraits::black()) {
0240                   // we do this to get shorter paths. check if we are nearer to
0241                   // the source as its parent is
0242                   if(is_closer_to_terminal(current_node, other_node)){
0243                     set_edge_to_parent(other_node, out_edge);
0244                     put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);
0245                     put(m_time_map, other_node, get(m_time_map, current_node));
0246                   }
0247                 } else{
0248                   BOOST_ASSERT(get_tree(other_node)==tColorTraits::white());
0249                   //kewl, found a path from one to the other search tree, return
0250                   // the connecting edge in src->sink dir
0251                   return std::make_pair(out_edge, true);
0252                 }
0253               }
0254             } //for all out-edges
0255           } //source-tree-growing
0256           else{
0257             BOOST_ASSERT(get_tree(current_node) == tColorTraits::white());
0258             out_edge_iterator ei, e_end;
0259             if(current_node != m_last_grow_vertex){
0260               m_last_grow_vertex = current_node;
0261               boost::tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g);
0262             }
0263             for(; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it){
0264               edge_descriptor in_edge = get(m_rev_edge_map, *m_last_grow_edge_it);
0265               if(get(m_res_cap_map, in_edge) > 0){ //check if there is capacity left
0266                 vertex_descriptor other_node = source(in_edge, m_g);
0267                 if(get_tree(other_node) == tColorTraits::gray()){ //it's a free node
0268                   set_tree(other_node, tColorTraits::white());      //acquire that node to our search tree
0269                   set_edge_to_parent(other_node, in_edge);          //set us as parent
0270                   add_active_node(other_node);                      //activate that node
0271                   put(m_dist_map, other_node, get(m_dist_map, current_node) + 1); //set its distance
0272                   put(m_time_map, other_node, get(m_time_map, current_node));//and time
0273                 } else if(get_tree(other_node) == tColorTraits::white()){
0274                   if(is_closer_to_terminal(current_node, other_node)){
0275                     //we are closer to the sink than its parent is, so we "adopt" him
0276                     set_edge_to_parent(other_node, in_edge);
0277                     put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);
0278                     put(m_time_map, other_node, get(m_time_map, current_node));
0279                   }
0280                 } else{
0281                   BOOST_ASSERT(get_tree(other_node)==tColorTraits::black());
0282                   //kewl, found a path from one to the other search tree,
0283                   // return the connecting edge in src->sink dir
0284                   return std::make_pair(in_edge, true);
0285                 }
0286               }
0287             } //for all out-edges
0288           } //sink-tree growing
0289 
0290           //all edges of that node are processed, and no more paths were found.
0291           // remove if from the front of the active queue
0292           finish_node(current_node);
0293         } //while active_nodes not empty
0294 
0295         //no active nodes anymore and no path found, we're done
0296         return std::make_pair(edge_descriptor(), false);
0297       }
0298 
0299       /**
0300        * augments path from s->t and updates residual graph
0301        * source(e, m_g) is the end of the path found in the source-tree
0302        * target(e, m_g) is the beginning of the path found in the sink-tree
0303        * this phase generates orphans on saturated edges, if the attached verts are
0304        * from different search-trees orphans are ordered in distance to
0305        * sink/source. first the farthest from the source are front_inserted into
0306        * the orphans list, and after that the sink-tree-orphans are
0307        * front_inserted. when going to adoption stage the orphans are popped_front,
0308        * and so we process the nearest verts to the terminals first
0309        */
0310       void augment(edge_descriptor e) {
0311         BOOST_ASSERT(get_tree(target(e, m_g)) == tColorTraits::white());
0312         BOOST_ASSERT(get_tree(source(e, m_g)) == tColorTraits::black());
0313         BOOST_ASSERT(m_orphans.empty());
0314 
0315         const tEdgeVal bottleneck = find_bottleneck(e);
0316         //now we push the found flow through the path
0317         //for each edge we saturate we have to look for the verts that belong to that edge, one of them becomes an orphans
0318         //now process the connecting edge
0319         put(m_res_cap_map, e, get(m_res_cap_map, e) - bottleneck);
0320         BOOST_ASSERT(get(m_res_cap_map, e) >= 0);
0321         put(m_res_cap_map, get(m_rev_edge_map, e), get(m_res_cap_map, get(m_rev_edge_map, e)) + bottleneck);
0322 
0323         //now we follow the path back to the source
0324         vertex_descriptor current_node = source(e, m_g);
0325         while(current_node != m_source){
0326           edge_descriptor pred = get_edge_to_parent(current_node);
0327           const tEdgeVal new_pred_cap = get(m_res_cap_map, pred) - bottleneck;
0328           put(m_res_cap_map, pred, new_pred_cap);
0329           BOOST_ASSERT(get(m_res_cap_map, pred) >= 0);
0330 
0331           const edge_descriptor pred_rev = get(m_rev_edge_map, pred);
0332           put(m_res_cap_map, pred_rev, get(m_res_cap_map, pred_rev) + bottleneck);
0333 
0334           if(new_pred_cap == 0){
0335             set_no_parent(current_node);
0336             m_orphans.push_front(current_node);
0337           }
0338           current_node = source(pred, m_g);
0339         }
0340         //then go forward in the sink-tree
0341         current_node = target(e, m_g);
0342         while(current_node != m_sink){
0343           edge_descriptor pred = get_edge_to_parent(current_node);
0344           const tEdgeVal new_pred_cap = get(m_res_cap_map, pred) - bottleneck;
0345           put(m_res_cap_map, pred, new_pred_cap);
0346           BOOST_ASSERT(get(m_res_cap_map, pred) >= 0);
0347 
0348           const edge_descriptor pred_rev = get(m_rev_edge_map, pred);
0349           put(m_res_cap_map, pred_rev, get(m_res_cap_map, pred_rev) + bottleneck);
0350 
0351           if(new_pred_cap == 0){
0352             set_no_parent(current_node);
0353             m_orphans.push_front(current_node);
0354           }
0355           current_node = target(pred, m_g);
0356         }
0357         //and add it to the max-flow
0358         m_flow += bottleneck;
0359       }
0360 
0361       /**
0362        * returns the bottleneck of a s->t path (end_of_path is last vertex in
0363        * source-tree, begin_of_path is first vertex in sink-tree)
0364        */
0365       inline tEdgeVal find_bottleneck(edge_descriptor e){
0366         BOOST_USING_STD_MIN();
0367         tEdgeVal minimum_cap = get(m_res_cap_map, e);
0368         vertex_descriptor current_node = source(e, m_g);
0369         //first go back in the source tree
0370         while(current_node != m_source){
0371           edge_descriptor pred = get_edge_to_parent(current_node);
0372           minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(minimum_cap, get(m_res_cap_map, pred));
0373           current_node = source(pred, m_g);
0374         }
0375         //then go forward in the sink-tree
0376         current_node = target(e, m_g);
0377         while(current_node != m_sink){
0378           edge_descriptor pred = get_edge_to_parent(current_node);
0379           minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(minimum_cap, get(m_res_cap_map, pred));
0380           current_node = target(pred, m_g);
0381         }
0382         return minimum_cap;
0383       }
0384 
0385       /**
0386        * rebuild search trees
0387        * empty the queue of orphans, and find new parents for them or just drop
0388        * them from the search trees
0389        */
0390       void adopt(){
0391         while(!m_orphans.empty() || !m_child_orphans.empty()){
0392           vertex_descriptor current_node;
0393           if(m_child_orphans.empty()){
0394             //get the next orphan from the main-queue  and remove it
0395             current_node = m_orphans.front();
0396             m_orphans.pop_front();
0397           } else{
0398             current_node = m_child_orphans.front();
0399             m_child_orphans.pop();
0400           }
0401           if(get_tree(current_node) == tColorTraits::black()){
0402             //we're in the source-tree
0403             tDistanceVal min_distance = (std::numeric_limits<tDistanceVal>::max)();
0404             edge_descriptor new_parent_edge;
0405             out_edge_iterator ei, e_end;
0406             for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
0407               const edge_descriptor in_edge = get(m_rev_edge_map, *ei);
0408               BOOST_ASSERT(target(in_edge, m_g) == current_node); //we should be the target of this edge
0409               if(get(m_res_cap_map, in_edge) > 0){
0410                 vertex_descriptor other_node = source(in_edge, m_g);
0411                 if(get_tree(other_node) == tColorTraits::black() && has_source_connect(other_node)){
0412                   if(get(m_dist_map, other_node) < min_distance){
0413                     min_distance = get(m_dist_map, other_node);
0414                     new_parent_edge = in_edge;
0415                   }
0416                 }
0417               }
0418             }
0419             if(min_distance != (std::numeric_limits<tDistanceVal>::max)()){
0420               set_edge_to_parent(current_node, new_parent_edge);
0421               put(m_dist_map, current_node, min_distance + 1);
0422               put(m_time_map, current_node, m_time);
0423             } else{
0424               put(m_time_map, current_node, 0);
0425               for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
0426                 edge_descriptor in_edge = get(m_rev_edge_map, *ei);
0427                 vertex_descriptor other_node = source(in_edge, m_g);
0428                 if(get_tree(other_node) == tColorTraits::black() && other_node != m_source){
0429                   if(get(m_res_cap_map, in_edge) > 0){
0430                     add_active_node(other_node);
0431                   }
0432                   if(has_parent(other_node) && source(get_edge_to_parent(other_node), m_g) == current_node){
0433                     //we are the parent of that node
0434                     //it has to find a new parent, too
0435                     set_no_parent(other_node);
0436                     m_child_orphans.push(other_node);
0437                   }
0438                 }
0439               }
0440               set_tree(current_node, tColorTraits::gray());
0441             } //no parent found
0442           } //source-tree-adoption
0443           else{
0444             //now we should be in the sink-tree, check that...
0445             BOOST_ASSERT(get_tree(current_node) == tColorTraits::white());
0446             out_edge_iterator ei, e_end;
0447             edge_descriptor new_parent_edge;
0448             tDistanceVal min_distance = (std::numeric_limits<tDistanceVal>::max)();
0449             for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
0450               const edge_descriptor out_edge = *ei;
0451               if(get(m_res_cap_map, out_edge) > 0){
0452                 const vertex_descriptor other_node = target(out_edge, m_g);
0453                 if(get_tree(other_node) == tColorTraits::white() && has_sink_connect(other_node))
0454                   if(get(m_dist_map, other_node) < min_distance){
0455                     min_distance = get(m_dist_map, other_node);
0456                     new_parent_edge = out_edge;
0457                   }
0458               }
0459             }
0460             if(min_distance != (std::numeric_limits<tDistanceVal>::max)()){
0461               set_edge_to_parent(current_node, new_parent_edge);
0462               put(m_dist_map, current_node, min_distance + 1);
0463               put(m_time_map, current_node, m_time);
0464             } else{
0465               put(m_time_map, current_node, 0);
0466               for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
0467                 const edge_descriptor out_edge = *ei;
0468                 const vertex_descriptor other_node = target(out_edge, m_g);
0469                 if(get_tree(other_node) == tColorTraits::white() && other_node != m_sink){
0470                   if(get(m_res_cap_map, out_edge) > 0){
0471                     add_active_node(other_node);
0472                   }
0473                   if(has_parent(other_node) && target(get_edge_to_parent(other_node), m_g) == current_node){
0474                     //we were it's parent, so it has to find a new one, too
0475                     set_no_parent(other_node);
0476                     m_child_orphans.push(other_node);
0477                   }
0478                 }
0479               }
0480               set_tree(current_node, tColorTraits::gray());
0481             } //no parent found
0482           } //sink-tree adoption
0483         } //while !orphans.empty()
0484       } //adopt
0485 
0486       /**
0487        * return next active vertex if there is one, otherwise a null_vertex
0488        */
0489       inline vertex_descriptor get_next_active_node(){
0490         while(true){
0491           if(m_active_nodes.empty())
0492             return graph_traits<Graph>::null_vertex();
0493           vertex_descriptor v = m_active_nodes.front();
0494 
0495       //if it has no parent, this node can't be active (if its not source or sink)
0496       if(!has_parent(v) && v != m_source && v != m_sink){
0497             m_active_nodes.pop();
0498             put(m_in_active_list_map, v, false);
0499           } else{
0500             BOOST_ASSERT(get_tree(v) == tColorTraits::black() || get_tree(v) == tColorTraits::white());
0501             return v;
0502           }
0503         }
0504       }
0505 
0506       /**
0507        * adds v as an active vertex, but only if its not in the list already
0508        */
0509       inline void add_active_node(vertex_descriptor v){
0510         BOOST_ASSERT(get_tree(v) != tColorTraits::gray());
0511         if(get(m_in_active_list_map, v)){
0512           if (m_last_grow_vertex == v) {
0513               m_last_grow_vertex = graph_traits<Graph>::null_vertex();
0514           }
0515           return;
0516         } else{
0517           put(m_in_active_list_map, v, true);
0518           m_active_nodes.push(v);
0519         }
0520       }
0521 
0522       /**
0523        * finish_node removes a node from the front of the active queue (its called in grow phase, if no more paths can be found using this node)
0524        */
0525       inline void finish_node(vertex_descriptor v){
0526         BOOST_ASSERT(m_active_nodes.front() == v);
0527         m_active_nodes.pop();
0528         put(m_in_active_list_map, v, false);
0529         m_last_grow_vertex = graph_traits<Graph>::null_vertex();
0530       }
0531 
0532       /**
0533        * removes a vertex from the queue of active nodes (actually this does nothing,
0534        * but checks if this node has no parent edge, as this is the criteria for
0535        * being no more active)
0536        */
0537       inline void remove_active_node(vertex_descriptor v){
0538         (void)v; // disable unused parameter warning
0539         BOOST_ASSERT(!has_parent(v));
0540       }
0541 
0542       /**
0543        * returns the search tree of v; tColorValue::black() for source tree,
0544        * white() for sink tree, gray() for no tree
0545        */
0546       inline tColorValue get_tree(vertex_descriptor v) const {
0547         return get(m_tree_map, v);
0548       }
0549 
0550       /**
0551        * sets search tree of v; tColorValue::black() for source tree, white()
0552        * for sink tree, gray() for no tree
0553        */
0554       inline void set_tree(vertex_descriptor v, tColorValue t){
0555         put(m_tree_map, v, t);
0556       }
0557 
0558       /**
0559        * returns edge to parent vertex of v;
0560        */
0561       inline edge_descriptor get_edge_to_parent(vertex_descriptor v) const{
0562         return get(m_pre_map, v);
0563       }
0564 
0565       /**
0566        * returns true if the edge stored in m_pre_map[v] is a valid entry
0567        */
0568       inline bool has_parent(vertex_descriptor v) const{
0569         return get(m_has_parent_map, v);
0570       }
0571 
0572       /**
0573        * sets edge to parent vertex of v;
0574        */
0575       inline void set_edge_to_parent(vertex_descriptor v, edge_descriptor f_edge_to_parent){
0576         BOOST_ASSERT(get(m_res_cap_map, f_edge_to_parent) > 0);
0577         put(m_pre_map, v, f_edge_to_parent);
0578         put(m_has_parent_map, v, true);
0579       }
0580 
0581       /**
0582        * removes the edge to parent of v (this is done by invalidating the
0583        * entry an additional map)
0584        */
0585       inline void set_no_parent(vertex_descriptor v){
0586         put(m_has_parent_map, v, false);
0587       }
0588 
0589       /**
0590        * checks if vertex v has a connect to the sink-vertex (@p m_sink)
0591        * @param v the vertex which is checked
0592        * @return true if a path to the sink was found, false if not
0593        */
0594       inline bool has_sink_connect(vertex_descriptor v){
0595         tDistanceVal current_distance = 0;
0596         vertex_descriptor current_vertex = v;
0597         while(true){
0598           if(get(m_time_map, current_vertex) == m_time){
0599             //we found a node which was already checked this round. use it for distance calculations
0600             current_distance += get(m_dist_map, current_vertex);
0601             break;
0602           }
0603           if(current_vertex == m_sink){
0604             put(m_time_map, m_sink, m_time);
0605             break;
0606           }
0607           if(has_parent(current_vertex)){
0608             //it has a parent, so get it
0609             current_vertex = target(get_edge_to_parent(current_vertex), m_g);
0610             ++current_distance;
0611           } else{
0612             //no path found
0613             return false;
0614           }
0615         }
0616         current_vertex=v;
0617         while(get(m_time_map, current_vertex) != m_time){
0618           put(m_dist_map, current_vertex, current_distance);
0619           --current_distance;
0620           put(m_time_map, current_vertex, m_time);
0621           current_vertex = target(get_edge_to_parent(current_vertex), m_g);
0622         }
0623         return true;
0624       }
0625 
0626       /**
0627        * checks if vertex v has a connect to the source-vertex (@p m_source)
0628        * @param v the vertex which is checked
0629        * @return true if a path to the source was found, false if not
0630        */
0631       inline bool has_source_connect(vertex_descriptor v){
0632         tDistanceVal current_distance = 0;
0633         vertex_descriptor current_vertex = v;
0634         while(true){
0635           if(get(m_time_map, current_vertex) == m_time){
0636             //we found a node which was already checked this round. use it for distance calculations
0637             current_distance += get(m_dist_map, current_vertex);
0638             break;
0639           }
0640           if(current_vertex == m_source){
0641             put(m_time_map, m_source, m_time);
0642             break;
0643           }
0644           if(has_parent(current_vertex)){
0645             //it has a parent, so get it
0646             current_vertex = source(get_edge_to_parent(current_vertex), m_g);
0647             ++current_distance;
0648           } else{
0649             //no path found
0650             return false;
0651           }
0652         }
0653         current_vertex=v;
0654         while(get(m_time_map, current_vertex) != m_time){
0655             put(m_dist_map, current_vertex, current_distance);
0656             --current_distance;
0657             put(m_time_map, current_vertex, m_time);
0658             current_vertex = source(get_edge_to_parent(current_vertex), m_g);
0659         }
0660         return true;
0661       }
0662 
0663       /**
0664        * returns true, if p is closer to a terminal than q
0665        */
0666       inline bool is_closer_to_terminal(vertex_descriptor p, vertex_descriptor q){
0667         //checks the timestamps first, to build no cycles, and after that the real distance
0668         return (get(m_time_map, q) <= get(m_time_map, p) &&
0669                 get(m_dist_map, q) > get(m_dist_map, p)+1);
0670       }
0671 
0672       ////////
0673       // member vars
0674       ////////
0675       Graph& m_g;
0676       IndexMap m_index_map;
0677       EdgeCapacityMap m_cap_map;
0678       ResidualCapacityEdgeMap m_res_cap_map;
0679       ReverseEdgeMap m_rev_edge_map;
0680       PredecessorMap m_pre_map; //stores paths found in the growth stage
0681       ColorMap m_tree_map; //maps each vertex into one of the two search tree or none (gray())
0682       DistanceMap m_dist_map; //stores distance to source/sink nodes
0683       vertex_descriptor m_source;
0684       vertex_descriptor m_sink;
0685 
0686       tQueue m_active_nodes;
0687       std::vector<bool> m_in_active_list_vec;
0688       iterator_property_map<std::vector<bool>::iterator, IndexMap> m_in_active_list_map;
0689 
0690       std::list<vertex_descriptor> m_orphans;
0691       tQueue m_child_orphans; // we use a second queue for child orphans, as they are FIFO processed
0692 
0693       std::vector<bool> m_has_parent_vec;
0694       iterator_property_map<std::vector<bool>::iterator, IndexMap> m_has_parent_map;
0695 
0696       std::vector<long> m_time_vec; //timestamp of each node, used for sink/source-path calculations
0697       iterator_property_map<std::vector<long>::iterator, IndexMap> m_time_map;
0698       tEdgeVal m_flow;
0699       long m_time;
0700       vertex_descriptor m_last_grow_vertex;
0701       out_edge_iterator m_last_grow_edge_it;
0702       out_edge_iterator m_last_grow_edge_end;
0703 };
0704 
0705 } //namespace boost::detail
0706 
0707 /**
0708   * non-named-parameter version, given everything
0709   * this is the catch all version
0710   */
0711 template<class Graph,
0712          class CapacityEdgeMap,
0713          class ResidualCapacityEdgeMap,
0714          class ReverseEdgeMap, class PredecessorMap,
0715          class ColorMap,
0716          class DistanceMap,
0717          class IndexMap>
0718 typename property_traits<CapacityEdgeMap>::value_type
0719 boykov_kolmogorov_max_flow(Graph& g,
0720                            CapacityEdgeMap cap,
0721                            ResidualCapacityEdgeMap res_cap,
0722                            ReverseEdgeMap rev_map,
0723                            PredecessorMap pre_map,
0724                            ColorMap color,
0725                            DistanceMap dist,
0726                            IndexMap idx,
0727                            typename graph_traits<Graph>::vertex_descriptor src,
0728                            typename graph_traits<Graph>::vertex_descriptor sink)
0729 {
0730   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
0731   typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;
0732 
0733   //as this method is the last one before we instantiate the solver, we do the concept checks here
0734   BOOST_CONCEPT_ASSERT(( VertexListGraphConcept<Graph> )); //to have vertices(), num_vertices(),
0735   BOOST_CONCEPT_ASSERT(( EdgeListGraphConcept<Graph> )); //to have edges()
0736   BOOST_CONCEPT_ASSERT(( IncidenceGraphConcept<Graph> )); //to have source(), target() and out_edges()
0737   BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<CapacityEdgeMap, edge_descriptor> )); //read flow-values from edges
0738   BOOST_CONCEPT_ASSERT(( ReadWritePropertyMapConcept<ResidualCapacityEdgeMap, edge_descriptor> )); //write flow-values to residuals
0739   BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<ReverseEdgeMap, edge_descriptor> )); //read out reverse edges
0740   BOOST_CONCEPT_ASSERT(( ReadWritePropertyMapConcept<PredecessorMap, vertex_descriptor> )); //store predecessor there
0741   BOOST_CONCEPT_ASSERT(( ReadWritePropertyMapConcept<ColorMap, vertex_descriptor> )); //write corresponding tree
0742   BOOST_CONCEPT_ASSERT(( ReadWritePropertyMapConcept<DistanceMap, vertex_descriptor> )); //write distance to source/sink
0743   BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<IndexMap, vertex_descriptor> )); //get index 0...|V|-1
0744   BOOST_ASSERT(num_vertices(g) >= 2 && src != sink);
0745 
0746   detail::bk_max_flow<
0747     Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, ReverseEdgeMap,
0748     PredecessorMap, ColorMap, DistanceMap, IndexMap
0749   > algo(g, cap, res_cap, rev_map, pre_map, color, dist, idx, src, sink);
0750 
0751   return algo.max_flow();
0752 }
0753 
0754 /**
0755  * non-named-parameter version, given capacity, residual_capacity,
0756  * reverse_edges, and an index map.
0757  */
0758 template<class Graph,
0759          class CapacityEdgeMap,
0760          class ResidualCapacityEdgeMap,
0761          class ReverseEdgeMap,
0762          class IndexMap>
0763 typename property_traits<CapacityEdgeMap>::value_type
0764 boykov_kolmogorov_max_flow(Graph& g,
0765                            CapacityEdgeMap cap,
0766                            ResidualCapacityEdgeMap res_cap,
0767                            ReverseEdgeMap rev,
0768                            IndexMap idx,
0769                            typename graph_traits<Graph>::vertex_descriptor src,
0770                            typename graph_traits<Graph>::vertex_descriptor sink)
0771 {
0772   typename graph_traits<Graph>::vertices_size_type n_verts = num_vertices(g);
0773   std::vector<typename graph_traits<Graph>::edge_descriptor> predecessor_vec(n_verts);
0774   std::vector<default_color_type> color_vec(n_verts);
0775   std::vector<typename graph_traits<Graph>::vertices_size_type> distance_vec(n_verts);
0776   return
0777     boykov_kolmogorov_max_flow(
0778       g, cap, res_cap, rev,
0779       make_iterator_property_map(predecessor_vec.begin(), idx),
0780       make_iterator_property_map(color_vec.begin(), idx),
0781       make_iterator_property_map(distance_vec.begin(), idx),
0782       idx, src, sink);
0783 }
0784 
0785 /**
0786  * non-named-parameter version, some given: capacity, residual_capacity,
0787  * reverse_edges, color_map and an index map. Use this if you are interested in
0788  * the minimum cut, as the color map provides that info.
0789  */
0790 template<class Graph,
0791          class CapacityEdgeMap,
0792          class ResidualCapacityEdgeMap,
0793          class ReverseEdgeMap,
0794          class ColorMap,
0795          class IndexMap>
0796 typename property_traits<CapacityEdgeMap>::value_type
0797 boykov_kolmogorov_max_flow(Graph& g,
0798                            CapacityEdgeMap cap,
0799                            ResidualCapacityEdgeMap res_cap,
0800                            ReverseEdgeMap rev,
0801                            ColorMap color,
0802                            IndexMap idx,
0803                            typename graph_traits<Graph>::vertex_descriptor src,
0804                            typename graph_traits<Graph>::vertex_descriptor sink)
0805 {
0806   typename graph_traits<Graph>::vertices_size_type n_verts = num_vertices(g);
0807   std::vector<typename graph_traits<Graph>::edge_descriptor> predecessor_vec(n_verts);
0808   std::vector<typename graph_traits<Graph>::vertices_size_type> distance_vec(n_verts);
0809   return
0810     boykov_kolmogorov_max_flow(
0811       g, cap, res_cap, rev,
0812       make_iterator_property_map(predecessor_vec.begin(), idx),
0813       color,
0814       make_iterator_property_map(distance_vec.begin(), idx),
0815       idx, src, sink);
0816 }
0817 
0818 /**
0819  * named-parameter version, some given
0820  */
0821 template<class Graph, class P, class T, class R>
0822 typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type
0823 boykov_kolmogorov_max_flow(Graph& g,
0824                            typename graph_traits<Graph>::vertex_descriptor src,
0825                            typename graph_traits<Graph>::vertex_descriptor sink,
0826                            const bgl_named_params<P, T, R>& params)
0827 {
0828   return
0829   boykov_kolmogorov_max_flow(
0830     g,
0831     choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),
0832     choose_pmap(get_param(params, edge_residual_capacity), g, edge_residual_capacity),
0833     choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),
0834     choose_pmap(get_param(params, vertex_predecessor), g, vertex_predecessor),
0835     choose_pmap(get_param(params, vertex_color), g, vertex_color),
0836     choose_pmap(get_param(params, vertex_distance), g, vertex_distance),
0837     choose_const_pmap(get_param(params, vertex_index), g, vertex_index),
0838     src, sink);
0839 }
0840 
0841 /**
0842  * named-parameter version, none given
0843  */
0844 template<class Graph>
0845 typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type
0846 boykov_kolmogorov_max_flow(Graph& g,
0847                            typename graph_traits<Graph>::vertex_descriptor src,
0848                            typename graph_traits<Graph>::vertex_descriptor sink)
0849 {
0850   bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
0851   return boykov_kolmogorov_max_flow(g, src, sink, params);
0852 }
0853 
0854 } // namespace boost
0855 
0856 #endif // BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
0857