boost/graph/boykov_kolmogorov_max_flow.hpp

   1 //  Copyright (c) 2006, Stephan Diederich
   2 //
   3 //  This code may be used under either of the following two licences:
   4 //
   5 //    Permission is hereby granted, free of charge, to any person
   6 //    obtaining a copy of this software and associated documentation
   7 //    files (the "Software"), to deal in the Software without
   8 //    restriction, including without limitation the rights to use,
   9 //    copy, modify, merge, publish, distribute, sublicense, and/or
  10 //    sell copies of the Software, and to permit persons to whom the
  11 //    Software is furnished to do so, subject to the following
  12 //    conditions:
  13 //
  14 //    The above copyright notice and this permission notice shall be
  15 //    included in all copies or substantial portions of the Software.
  16 //
  17 //    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  18 //    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  19 //    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  20 //    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  21 //    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  22 //    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  23 //    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  24 //    OTHER DEALINGS IN THE SOFTWARE. OF SUCH DAMAGE.
  25 //
  26 //  Or:
  27 //
  28 //    Distributed under the Boost Software License, Version 1.0.
  29 //    (See accompanying file LICENSE_1_0.txt or copy at
  30 //    http://www.boost.org/LICENSE_1_0.txt)
  31
  32 #ifndef BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
  33 #define BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
  34
  35 #include <boost/config.hpp>
  36 #include <boost/assert.hpp>
  37 #include <vector>
  38 #include <list>
  39 #include <utility>
  40 #include <iosfwd>
  41 #include <algorithm> // for std::min and std::max
  42
  43 #include <boost/pending/queue.hpp>
  44 #include <boost/limits.hpp>
  45 #include <boost/property_map/property_map.hpp>
  46 #include <boost/none_t.hpp>
  47 #include <boost/graph/graph_concepts.hpp>
  48 #include <boost/graph/named_function_params.hpp>
  49 #include <boost/graph/lookup_edge.hpp>
  50 #include <boost/concept/assert.hpp>
  51
  52 // The algorithm impelemented here is described in:
  53 //
  54 // Boykov, Y., Kolmogorov, V. "An Experimental Comparison of Min-Cut/Max-Flow
  55 // Algorithms for Energy Minimization in Vision", In IEEE Transactions on
  56 // Pattern Analysis and Machine Intelligence, vol. 26, no. 9, pp. 1124-1137,
  57 // Sep 2004.
  58 //
  59 // For further reading, also see:
  60 //
  61 // Kolmogorov, V. "Graph Based Algorithms for Scene Reconstruction from Two or
  62 // More Views". PhD thesis, Cornell University, Sep 2003.
  63
  64 namespace boost {
  65
  66 namespace detail {
  67
  68 template <class Graph,
  69           class EdgeCapacityMap,
  70           class ResidualCapacityEdgeMap,
  71           class ReverseEdgeMap,
  72           class PredecessorMap,
  73           class ColorMap,
  74           class DistanceMap,
  75           class IndexMap>
  76 class bk_max_flow {
  77   typedef typename property_traits<EdgeCapacityMap>::value_type tEdgeVal;
  78   typedef graph_traits<Graph> tGraphTraits;
  79   typedef typename tGraphTraits::vertex_iterator vertex_iterator;
  80   typedef typename tGraphTraits::vertex_descriptor vertex_descriptor;
  81   typedef typename tGraphTraits::edge_descriptor edge_descriptor;
  82   typedef typename tGraphTraits::edge_iterator edge_iterator;
  83   typedef typename tGraphTraits::out_edge_iterator out_edge_iterator;
  84   typedef boost::queue<vertex_descriptor> tQueue;                               //queue of vertices, used in adoption-stage
  85   typedef typename property_traits<ColorMap>::value_type tColorValue;
  86   typedef color_traits<tColorValue> tColorTraits;
  87   typedef typename property_traits<DistanceMap>::value_type tDistanceVal;
  88
  89     public:
  90       bk_max_flow(Graph& g,
  91                   EdgeCapacityMap cap,
  92                   ResidualCapacityEdgeMap res,
  93                   ReverseEdgeMap rev,
  94                   PredecessorMap pre,
  95                   ColorMap color,
  96                   DistanceMap dist,
  97                   IndexMap idx,
  98                   vertex_descriptor src,
  99                   vertex_descriptor sink):
 100       m_g(g),
 101       m_index_map(idx),
 102       m_cap_map(cap),
 103       m_res_cap_map(res),
 104       m_rev_edge_map(rev),
 105       m_pre_map(pre),
 106       m_tree_map(color),
 107       m_dist_map(dist),
 108       m_source(src),
 109       m_sink(sink),
 110       m_active_nodes(),
 111       m_in_active_list_vec(num_vertices(g), false),
 112       m_in_active_list_map(make_iterator_property_map(m_in_active_list_vec.begin(), m_index_map)),
 113       m_has_parent_vec(num_vertices(g), false),
 114       m_has_parent_map(make_iterator_property_map(m_has_parent_vec.begin(), m_index_map)),
 115       m_time_vec(num_vertices(g), 0),
 116       m_time_map(make_iterator_property_map(m_time_vec.begin(), m_index_map)),
 117       m_flow(0),
 118       m_time(1),
 119       m_last_grow_vertex(graph_traits<Graph>::null_vertex()){
 120         // initialize the color-map with gray-values
 121         vertex_iterator vi, v_end;
 122         for(boost::tie(vi, v_end) = vertices(m_g); vi != v_end; ++vi){
 123           set_tree(*vi, tColorTraits::gray());
 124         }
 125         // Initialize flow to zero which means initializing
 126         // the residual capacity equal to the capacity
 127         edge_iterator ei, e_end;
 128         for(boost::tie(ei, e_end) = edges(m_g); ei != e_end; ++ei) {
 129           put(m_res_cap_map, *ei, get(m_cap_map, *ei));
 130           BOOST_ASSERT(get(m_rev_edge_map, get(m_rev_edge_map, *ei)) == *ei); //check if the reverse edge map is build up properly
 131         }
 132         //init the search trees with the two terminals
 133         set_tree(m_source, tColorTraits::black());
 134         set_tree(m_sink, tColorTraits::white());
 135         put(m_time_map, m_source, 1);
 136         put(m_time_map, m_sink, 1);
 137       }
 138
 139       tEdgeVal max_flow(){
 140         //augment direct paths from SOURCE->SINK and SOURCE->VERTEX->SINK
 141         augment_direct_paths();
 142         //start the main-loop
 143         while(true){
 144           bool path_found;
 145           edge_descriptor connecting_edge;
 146           boost::tie(connecting_edge, path_found) = grow(); //find a path from source to sink
 147           if(!path_found){
 148             //we're finished, no more paths were found
 149             break;
 150           }
 151           ++m_time;
 152           augment(connecting_edge); //augment that path
 153           adopt(); //rebuild search tree structure
 154         }
 155         return m_flow;
 156       }
 157
 158       // the complete class is protected, as we want access to members in
 159       // derived test-class (see test/boykov_kolmogorov_max_flow_test.cpp)
 160     protected:
 161       void augment_direct_paths(){
 162         // in a first step, we augment all direct paths from source->NODE->sink
 163         // and additionally paths from source->sink. This improves especially
 164         // graphcuts for segmentation, as most of the nodes have source/sink
 165         // connects but shouldn't have an impact on other maxflow problems
 166         // (this is done in grow() anyway)
 167         out_edge_iterator ei, e_end;
 168         for(boost::tie(ei, e_end) = out_edges(m_source, m_g); ei != e_end; ++ei){
 169           edge_descriptor from_source = *ei;
 170           vertex_descriptor current_node = target(from_source, m_g);
 171           if(current_node == m_sink){
 172             tEdgeVal cap = get(m_res_cap_map, from_source);
 173             put(m_res_cap_map, from_source, 0);
 174             m_flow += cap;
 175             continue;
 176           }
 177           edge_descriptor to_sink;
 178           bool is_there;
 179           boost::tie(to_sink, is_there) = lookup_edge(current_node, m_sink, m_g);
 180           if(is_there){
 181             tEdgeVal cap_from_source = get(m_res_cap_map, from_source);
 182             tEdgeVal cap_to_sink = get(m_res_cap_map, to_sink);
 183             if(cap_from_source > cap_to_sink){
 184               set_tree(current_node, tColorTraits::black());
 185               add_active_node(current_node);
 186               set_edge_to_parent(current_node, from_source);
 187               put(m_dist_map, current_node, 1);
 188               put(m_time_map, current_node, 1);
 189               // add stuff to flow and update residuals. we dont need to
 190               // update reverse_edges, as incoming/outgoing edges to/from
 191               // source/sink don't count for max-flow
 192               put(m_res_cap_map, from_source, get(m_res_cap_map, from_source) - cap_to_sink);
 193               put(m_res_cap_map, to_sink, 0);
 194               m_flow += cap_to_sink;
 195             } else if(cap_to_sink > 0){
 196               set_tree(current_node, tColorTraits::white());
 197               add_active_node(current_node);
 198               set_edge_to_parent(current_node, to_sink);
 199               put(m_dist_map, current_node, 1);
 200               put(m_time_map, current_node, 1);
 201               // add stuff to flow and update residuals. we dont need to update
 202               // reverse_edges, as incoming/outgoing edges to/from source/sink
 203               // don't count for max-flow
 204               put(m_res_cap_map, to_sink, get(m_res_cap_map, to_sink) - cap_from_source);
 205               put(m_res_cap_map, from_source, 0);
 206               m_flow += cap_from_source;
 207             }
 208           } else if(get(m_res_cap_map, from_source)){
 209             // there is no sink connect, so we can't augment this path, but to
 210             // avoid adding m_source to the active nodes, we just activate this
 211             // node and set the approciate things
 212             set_tree(current_node, tColorTraits::black());
 213             set_edge_to_parent(current_node, from_source);
 214             put(m_dist_map, current_node, 1);
 215             put(m_time_map, current_node, 1);
 216             add_active_node(current_node);
 217           }
 218         }
 219         for(boost::tie(ei, e_end) = out_edges(m_sink, m_g); ei != e_end; ++ei){
 220           edge_descriptor to_sink = get(m_rev_edge_map, *ei);
 221           vertex_descriptor current_node = source(to_sink, m_g);
 222           if(get(m_res_cap_map, to_sink)){
 223             set_tree(current_node, tColorTraits::white());
 224             set_edge_to_parent(current_node, to_sink);
 225             put(m_dist_map, current_node, 1);
 226             put(m_time_map, current_node, 1);
 227             add_active_node(current_node);
 228           }
 229         }
 230       }
 231
 232       /**
 233        * Returns a pair of an edge and a boolean. if the bool is true, the
 234        * edge is a connection of a found path from s->t , read "the link" and
 235        * source(returnVal, m_g) is the end of the path found in the source-tree
 236        * target(returnVal, m_g) is the beginning of the path found in the sink-tree
 237        */
 238       std::pair<edge_descriptor, bool> grow(){
 239         BOOST_ASSERT(m_orphans.empty());
 240         vertex_descriptor current_node;
 241         while((current_node = get_next_active_node()) != graph_traits<Graph>::null_vertex()){ //if there is one
 242           BOOST_ASSERT(get_tree(current_node) != tColorTraits::gray() &&
 243                        (has_parent(current_node) ||
 244                          current_node == m_source ||
 245                          current_node == m_sink));
 246
 247           if(get_tree(current_node) == tColorTraits::black()){
 248             //source tree growing
 249             out_edge_iterator ei, e_end;
 250             if(current_node != m_last_grow_vertex){
 251               m_last_grow_vertex = current_node;
 252               boost::tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g);
 253             }
 254             for(; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it) {
 255               edge_descriptor out_edge = *m_last_grow_edge_it;
 256               if(get(m_res_cap_map, out_edge) > 0){ //check if we have capacity left on this edge
 257                 vertex_descriptor other_node = target(out_edge, m_g);
 258                 if(get_tree(other_node) == tColorTraits::gray()){ //it's a free node
 259                   set_tree(other_node, tColorTraits::black()); //aquire other node to our search tree
 260                   set_edge_to_parent(other_node, out_edge);   //set us as parent
 261                   put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);  //and update the distance-heuristic
 262                   put(m_time_map, other_node, get(m_time_map, current_node));
 263                   add_active_node(other_node);
 264                 } else if(get_tree(other_node) == tColorTraits::black()) {
 265                   // we do this to get shorter paths. check if we are nearer to
 266                   // the source as its parent is
 267                   if(is_closer_to_terminal(current_node, other_node)){
 268                     set_edge_to_parent(other_node, out_edge);
 269                     put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);
 270                     put(m_time_map, other_node, get(m_time_map, current_node));
 271                   }
 272                 } else{
 273                   BOOST_ASSERT(get_tree(other_node)==tColorTraits::white());
 274                   //kewl, found a path from one to the other search tree, return
 275                   // the connecting edge in src->sink dir
 276                   return std::make_pair(out_edge, true);
 277                 }
 278               }
 279             } //for all out-edges
 280           } //source-tree-growing
 281           else{
 282             BOOST_ASSERT(get_tree(current_node) == tColorTraits::white());
 283             out_edge_iterator ei, e_end;
 284             if(current_node != m_last_grow_vertex){
 285               m_last_grow_vertex = current_node;
 286               boost::tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g);
 287             }
 288             for(; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it){
 289               edge_descriptor in_edge = get(m_rev_edge_map, *m_last_grow_edge_it);
 290               if(get(m_res_cap_map, in_edge) > 0){ //check if there is capacity left
 291                 vertex_descriptor other_node = source(in_edge, m_g);
 292                 if(get_tree(other_node) == tColorTraits::gray()){ //it's a free node
 293                   set_tree(other_node, tColorTraits::white());      //aquire that node to our search tree
 294                   set_edge_to_parent(other_node, in_edge);          //set us as parent
 295                   add_active_node(other_node);                      //activate that node
 296                   put(m_dist_map, other_node, get(m_dist_map, current_node) + 1); //set its distance
 297                   put(m_time_map, other_node, get(m_time_map, current_node));//and time
 298                 } else if(get_tree(other_node) == tColorTraits::white()){
 299                   if(is_closer_to_terminal(current_node, other_node)){
 300                     //we are closer to the sink than its parent is, so we "adopt" him
 301                     set_edge_to_parent(other_node, in_edge);
 302                     put(m_dist_map, other_node, get(m_dist_map, current_node) + 1);
 303                     put(m_time_map, other_node, get(m_time_map, current_node));
 304                   }
 305                 } else{
 306                   BOOST_ASSERT(get_tree(other_node)==tColorTraits::black());
 307                   //kewl, found a path from one to the other search tree,
 308                   // return the connecting edge in src->sink dir
 309                   return std::make_pair(in_edge, true);
 310                 }
 311               }
 312             } //for all out-edges
 313           } //sink-tree growing
 314
 315           //all edges of that node are processed, and no more paths were found.
 316           // remove if from the front of the active queue
 317           finish_node(current_node);
 318         } //while active_nodes not empty
 319
 320         //no active nodes anymore and no path found, we're done
 321         return std::make_pair(edge_descriptor(), false);
 322       }
 323
 324       /**
 325        * augments path from s->t and updates residual graph
 326        * source(e, m_g) is the end of the path found in the source-tree
 327        * target(e, m_g) is the beginning of the path found in the sink-tree
 328        * this phase generates orphans on satured edges, if the attached verts are
 329        * from different search-trees orphans are ordered in distance to
 330        * sink/source. first the farest from the source are front_inserted into
 331        * the orphans list, and after that the sink-tree-orphans are
 332        * front_inserted. when going to adoption stage the orphans are popped_front,
 333        * and so we process the nearest verts to the terminals first
 334        */
 335       void augment(edge_descriptor e) {
 336         BOOST_ASSERT(get_tree(target(e, m_g)) == tColorTraits::white());
 337         BOOST_ASSERT(get_tree(source(e, m_g)) == tColorTraits::black());
 338         BOOST_ASSERT(m_orphans.empty());
 339
 340         const tEdgeVal bottleneck = find_bottleneck(e);
 341         //now we push the found flow through the path
 342         //for each edge we saturate we have to look for the verts that belong to that edge, one of them becomes an orphans
 343         //now process the connecting edge
 344         put(m_res_cap_map, e, get(m_res_cap_map, e) - bottleneck);
 345         BOOST_ASSERT(get(m_res_cap_map, e) >= 0);
 346         put(m_res_cap_map, get(m_rev_edge_map, e), get(m_res_cap_map, get(m_rev_edge_map, e)) + bottleneck);
 347
 348         //now we follow the path back to the source
 349         vertex_descriptor current_node = source(e, m_g);
 350         while(current_node != m_source){
 351           edge_descriptor pred = get_edge_to_parent(current_node);
 352           put(m_res_cap_map, pred, get(m_res_cap_map, pred) - bottleneck);
 353           BOOST_ASSERT(get(m_res_cap_map, pred) >= 0);
 354           put(m_res_cap_map, get(m_rev_edge_map, pred), get(m_res_cap_map, get(m_rev_edge_map, pred)) + bottleneck);
 355           if(get(m_res_cap_map, pred) == 0){
 356             set_no_parent(current_node);
 357             m_orphans.push_front(current_node);
 358           }
 359           current_node = source(pred, m_g);
 360         }
 361         //then go forward in the sink-tree
 362         current_node = target(e, m_g);
 363         while(current_node != m_sink){
 364           edge_descriptor pred = get_edge_to_parent(current_node);
 365           put(m_res_cap_map, pred, get(m_res_cap_map, pred) - bottleneck);
 366           BOOST_ASSERT(get(m_res_cap_map, pred) >= 0);
 367           put(m_res_cap_map, get(m_rev_edge_map, pred), get(m_res_cap_map, get(m_rev_edge_map, pred)) + bottleneck);
 368           if(get(m_res_cap_map, pred) == 0){
 369             set_no_parent(current_node);
 370             m_orphans.push_front(current_node);
 371           }
 372           current_node = target(pred, m_g);
 373         }
 374         //and add it to the max-flow
 375         m_flow += bottleneck;
 376       }
 377
 378       /**
 379        * returns the bottleneck of a s->t path (end_of_path is last vertex in
 380        * source-tree, begin_of_path is first vertex in sink-tree)
 381        */
 382       inline tEdgeVal find_bottleneck(edge_descriptor e){
 383         BOOST_USING_STD_MIN();
 384         tEdgeVal minimum_cap = get(m_res_cap_map, e);
 385         vertex_descriptor current_node = source(e, m_g);
 386         //first go back in the source tree
 387         while(current_node != m_source){
 388           edge_descriptor pred = get_edge_to_parent(current_node);
 389           minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(minimum_cap, get(m_res_cap_map, pred));
 390           current_node = source(pred, m_g);
 391         }
 392         //then go forward in the sink-tree
 393         current_node = target(e, m_g);
 394         while(current_node != m_sink){
 395           edge_descriptor pred = get_edge_to_parent(current_node);
 396           minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(minimum_cap, get(m_res_cap_map, pred));
 397           current_node = target(pred, m_g);
 398         }
 399         return minimum_cap;
 400       }
 401
 402       /**
 403        * rebuild search trees
 404        * empty the queue of orphans, and find new parents for them or just drop
 405        * them from the search trees
 406        */
 407       void adopt(){
 408         while(!m_orphans.empty() || !m_child_orphans.empty()){
 409           vertex_descriptor current_node;
 410           if(m_child_orphans.empty()){
 411             //get the next orphan from the main-queue  and remove it
 412             current_node = m_orphans.front();
 413             m_orphans.pop_front();
 414           } else{
 415             current_node = m_child_orphans.front();
 416             m_child_orphans.pop();
 417           }
 418           if(get_tree(current_node) == tColorTraits::black()){
 419             //we're in the source-tree
 420             tDistanceVal min_distance = (std::numeric_limits<tDistanceVal>::max)();
 421             edge_descriptor new_parent_edge;
 422             out_edge_iterator ei, e_end;
 423             for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
 424               const edge_descriptor in_edge = get(m_rev_edge_map, *ei);
 425               BOOST_ASSERT(target(in_edge, m_g) == current_node); //we should be the target of this edge
 426               if(get(m_res_cap_map, in_edge) > 0){
 427                 vertex_descriptor other_node = source(in_edge, m_g);
 428                 if(get_tree(other_node) == tColorTraits::black() && has_source_connect(other_node)){
 429                   if(get(m_dist_map, other_node) < min_distance){
 430                     min_distance = get(m_dist_map, other_node);
 431                     new_parent_edge = in_edge;
 432                   }
 433                 }
 434               }
 435             }
 436             if(min_distance != (std::numeric_limits<tDistanceVal>::max)()){
 437               set_edge_to_parent(current_node, new_parent_edge);
 438               put(m_dist_map, current_node, min_distance + 1);
 439               put(m_time_map, current_node, m_time);
 440             } else{
 441               put(m_time_map, current_node, 0);
 442               for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
 443                 edge_descriptor in_edge = get(m_rev_edge_map, *ei);
 444                 vertex_descriptor other_node = source(in_edge, m_g);
 445                 if(get_tree(other_node) == tColorTraits::black() && other_node != m_source){
 446                   if(get(m_res_cap_map, in_edge) > 0){
 447                     add_active_node(other_node);
 448                   }
 449                   if(has_parent(other_node) && source(get_edge_to_parent(other_node), m_g) == current_node){
 450                     //we are the parent of that node
 451                     //it has to find a new parent, too
 452                     set_no_parent(other_node);
 453                     m_child_orphans.push(other_node);
 454                   }
 455                 }
 456               }
 457               set_tree(current_node, tColorTraits::gray());
 458             } //no parent found
 459           } //source-tree-adoption
 460           else{
 461             //now we should be in the sink-tree, check that...
 462             BOOST_ASSERT(get_tree(current_node) == tColorTraits::white());
 463             out_edge_iterator ei, e_end;
 464             edge_descriptor new_parent_edge;
 465             tDistanceVal min_distance = (std::numeric_limits<tDistanceVal>::max)();
 466             for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
 467               const edge_descriptor out_edge = *ei;
 468               if(get(m_res_cap_map, out_edge) > 0){
 469                 const vertex_descriptor other_node = target(out_edge, m_g);
 470                 if(get_tree(other_node) == tColorTraits::white() && has_sink_connect(other_node))
 471                   if(get(m_dist_map, other_node) < min_distance){
 472                     min_distance = get(m_dist_map, other_node);
 473                     new_parent_edge = out_edge;
 474                   }
 475               }
 476             }
 477             if(min_distance != (std::numeric_limits<tDistanceVal>::max)()){
 478               set_edge_to_parent(current_node, new_parent_edge);
 479               put(m_dist_map, current_node, min_distance + 1);
 480               put(m_time_map, current_node, m_time);
 481             } else{
 482               put(m_time_map, current_node, 0);
 483               for(boost::tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){
 484                 const edge_descriptor out_edge = *ei;
 485                 const vertex_descriptor other_node = target(out_edge, m_g);
 486                 if(get_tree(other_node) == tColorTraits::white() && other_node != m_sink){
 487                   if(get(m_res_cap_map, out_edge) > 0){
 488                     add_active_node(other_node);
 489                   }
 490                   if(has_parent(other_node) && target(get_edge_to_parent(other_node), m_g) == current_node){
 491                     //we were it's parent, so it has to find a new one, too
 492                     set_no_parent(other_node);
 493                     m_child_orphans.push(other_node);
 494                   }
 495                 }
 496               }
 497               set_tree(current_node, tColorTraits::gray());
 498             } //no parent found
 499           } //sink-tree adoption
 500         } //while !orphans.empty()
 501       } //adopt
 502
 503       /**
 504        * return next active vertex if there is one, otherwise a null_vertex
 505        */
 506       inline vertex_descriptor get_next_active_node(){
 507         while(true){
 508           if(m_active_nodes.empty())
 509             return graph_traits<Graph>::null_vertex();
 510           vertex_descriptor v = m_active_nodes.front();
 511
 512       //if it has no parent, this node can't be active (if its not source or sink)
 513       if(!has_parent(v) && v != m_source && v != m_sink){
 514             m_active_nodes.pop();
 515             put(m_in_active_list_map, v, false);
 516           } else{
 517             BOOST_ASSERT(get_tree(v) == tColorTraits::black() || get_tree(v) == tColorTraits::white());
 518             return v;
 519           }
 520         }
 521       }
 522
 523       /**
 524        * adds v as an active vertex, but only if its not in the list already
 525        */
 526       inline void add_active_node(vertex_descriptor v){
 527         BOOST_ASSERT(get_tree(v) != tColorTraits::gray());
 528         if(get(m_in_active_list_map, v)){
 529           if (m_last_grow_vertex == v) {
 530               m_last_grow_vertex = graph_traits<Graph>::null_vertex();
 531           }
 532           return;
 533         } else{
 534           put(m_in_active_list_map, v, true);
 535           m_active_nodes.push(v);
 536         }
 537       }
 538
 539       /**
 540        * 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)
 541        */
 542       inline void finish_node(vertex_descriptor v){
 543         BOOST_ASSERT(m_active_nodes.front() == v);
 544         m_active_nodes.pop();
 545         put(m_in_active_list_map, v, false);
 546         m_last_grow_vertex = graph_traits<Graph>::null_vertex();
 547       }
 548
 549       /**
 550        * removes a vertex from the queue of active nodes (actually this does nothing,
 551        * but checks if this node has no parent edge, as this is the criteria for
 552        * being no more active)
 553        */
 554       inline void remove_active_node(vertex_descriptor v){
 555         BOOST_ASSERT(!has_parent(v));
 556       }
 557
 558       /**
 559        * returns the search tree of v; tColorValue::black() for source tree,
 560        * white() for sink tree, gray() for no tree
 561        */
 562       inline tColorValue get_tree(vertex_descriptor v) const {
 563         return get(m_tree_map, v);
 564       }
 565
 566       /**
 567        * sets search tree of v; tColorValue::black() for source tree, white()
 568        * for sink tree, gray() for no tree
 569        */
 570       inline void set_tree(vertex_descriptor v, tColorValue t){
 571         put(m_tree_map, v, t);
 572       }
 573
 574       /**
 575        * returns edge to parent vertex of v;
 576        */
 577       inline edge_descriptor get_edge_to_parent(vertex_descriptor v) const{
 578         return get(m_pre_map, v);
 579       }
 580
 581       /**
 582        * returns true if the edge stored in m_pre_map[v] is a valid entry
 583        */
 584       inline bool has_parent(vertex_descriptor v) const{
 585         return get(m_has_parent_map, v);
 586       }
 587
 588       /**
 589        * sets edge to parent vertex of v;
 590        */
 591       inline void set_edge_to_parent(vertex_descriptor v, edge_descriptor f_edge_to_parent){
 592         BOOST_ASSERT(get(m_res_cap_map, f_edge_to_parent) > 0);
 593         put(m_pre_map, v, f_edge_to_parent);
 594         put(m_has_parent_map, v, true);
 595       }
 596
 597       /**
 598        * removes the edge to parent of v (this is done by invalidating the
 599        * entry an additional map)
 600        */
 601       inline void set_no_parent(vertex_descriptor v){
 602         put(m_has_parent_map, v, false);
 603       }
 604
 605       /**
 606        * checks if vertex v has a connect to the sink-vertex (@var m_sink)
 607        * @param v the vertex which is checked
 608        * @return true if a path to the sink was found, false if not
 609        */
 610       inline bool has_sink_connect(vertex_descriptor v){
 611         tDistanceVal current_distance = 0;
 612         vertex_descriptor current_vertex = v;
 613         while(true){
 614           if(get(m_time_map, current_vertex) == m_time){
 615             //we found a node which was already checked this round. use it for distance calculations
 616             current_distance += get(m_dist_map, current_vertex);
 617             break;
 618           }
 619           if(current_vertex == m_sink){
 620             put(m_time_map, m_sink, m_time);
 621             break;
 622           }
 623           if(has_parent(current_vertex)){
 624             //it has a parent, so get it
 625             current_vertex = target(get_edge_to_parent(current_vertex), m_g);
 626             ++current_distance;
 627           } else{
 628             //no path found
 629             return false;
 630           }
 631         }
 632         current_vertex=v;
 633         while(get(m_time_map, current_vertex) != m_time){
 634           put(m_dist_map, current_vertex, current_distance);
 635           --current_distance;
 636           put(m_time_map, current_vertex, m_time);
 637           current_vertex = target(get_edge_to_parent(current_vertex), m_g);
 638         }
 639         return true;
 640       }
 641
 642       /**
 643        * checks if vertex v has a connect to the source-vertex (@var m_source)
 644        * @param v the vertex which is checked
 645        * @return true if a path to the source was found, false if not
 646        */
 647       inline bool has_source_connect(vertex_descriptor v){
 648         tDistanceVal current_distance = 0;
 649         vertex_descriptor current_vertex = v;
 650         while(true){
 651           if(get(m_time_map, current_vertex) == m_time){
 652             //we found a node which was already checked this round. use it for distance calculations
 653             current_distance += get(m_dist_map, current_vertex);
 654             break;
 655           }
 656           if(current_vertex == m_source){
 657             put(m_time_map, m_source, m_time);
 658             break;
 659           }
 660           if(has_parent(current_vertex)){
 661             //it has a parent, so get it
 662             current_vertex = source(get_edge_to_parent(current_vertex), m_g);
 663             ++current_distance;
 664           } else{
 665             //no path found
 666             return false;
 667           }
 668         }
 669         current_vertex=v;
 670         while(get(m_time_map, current_vertex) != m_time){
 671             put(m_dist_map, current_vertex, current_distance);
 672             --current_distance;
 673             put(m_time_map, current_vertex, m_time);
 674             current_vertex = source(get_edge_to_parent(current_vertex), m_g);
 675         }
 676         return true;
 677       }
 678
 679       /**
 680        * returns true, if p is closer to a terminal than q
 681        */
 682       inline bool is_closer_to_terminal(vertex_descriptor p, vertex_descriptor q){
 683         //checks the timestamps first, to build no cycles, and after that the real distance
 684         return (get(m_time_map, q) <= get(m_time_map, p) &&
 685                 get(m_dist_map, q) > get(m_dist_map, p)+1);
 686       }
 687
 688       ////////
 689       // member vars
 690       ////////
 691       Graph& m_g;
 692       IndexMap m_index_map;
 693       EdgeCapacityMap m_cap_map;
 694       ResidualCapacityEdgeMap m_res_cap_map;
 695       ReverseEdgeMap m_rev_edge_map;
 696       PredecessorMap m_pre_map; //stores paths found in the growth stage
 697       ColorMap m_tree_map; //maps each vertex into one of the two search tree or none (gray())
 698       DistanceMap m_dist_map; //stores distance to source/sink nodes
 699       vertex_descriptor m_source;
 700       vertex_descriptor m_sink;
 701
 702       tQueue m_active_nodes;
 703       std::vector<bool> m_in_active_list_vec;
 704       iterator_property_map<std::vector<bool>::iterator, IndexMap> m_in_active_list_map;
 705
 706       std::list<vertex_descriptor> m_orphans;
 707       tQueue m_child_orphans; // we use a second queuqe for child orphans, as they are FIFO processed
 708
 709       std::vector<bool> m_has_parent_vec;
 710       iterator_property_map<std::vector<bool>::iterator, IndexMap> m_has_parent_map;
 711
 712       std::vector<long> m_time_vec; //timestamp of each node, used for sink/source-path calculations
 713       iterator_property_map<std::vector<long>::iterator, IndexMap> m_time_map;
 714       tEdgeVal m_flow;
 715       long m_time;
 716       vertex_descriptor m_last_grow_vertex;
 717       out_edge_iterator m_last_grow_edge_it;
 718       out_edge_iterator m_last_grow_edge_end;
 719 };
 720
 721 } //namespace boost::detail
 722
 723 /**
 724   * non-named-parameter version, given everything
 725   * this is the catch all version
 726   */
 727 template<class Graph,
 728          class CapacityEdgeMap,
 729          class ResidualCapacityEdgeMap,
 730          class ReverseEdgeMap, class PredecessorMap,
 731          class ColorMap,
 732          class DistanceMap,
 733          class IndexMap>
 734 typename property_traits<CapacityEdgeMap>::value_type
 735 boykov_kolmogorov_max_flow(Graph& g,
 736                            CapacityEdgeMap cap,
 737                            ResidualCapacityEdgeMap res_cap,
 738                            ReverseEdgeMap rev_map,
 739                            PredecessorMap pre_map,
 740                            ColorMap color,
 741                            DistanceMap dist,
 742                            IndexMap idx,
 743                            typename graph_traits<Graph>::vertex_descriptor src,
 744                            typename graph_traits<Graph>::vertex_descriptor sink)
 745 {
 746   typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;
 747   typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;
 748
 749   //as this method is the last one before we instantiate the solver, we do the concept checks here
 750   BOOST_CONCEPT_ASSERT(( VertexListGraphConcept<Graph> )); //to have vertices(), num_vertices(),
 751   BOOST_CONCEPT_ASSERT(( EdgeListGraphConcept<Graph> )); //to have edges()
 752   BOOST_CONCEPT_ASSERT(( IncidenceGraphConcept<Graph> )); //to have source(), target() and out_edges()
 753   BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<CapacityEdgeMap, edge_descriptor> )); //read flow-values from edges
 754   BOOST_CONCEPT_ASSERT(( ReadWritePropertyMapConcept<ResidualCapacityEdgeMap, edge_descriptor> )); //write flow-values to residuals
 755   BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<ReverseEdgeMap, edge_descriptor> )); //read out reverse edges
 756   BOOST_CONCEPT_ASSERT(( ReadWritePropertyMapConcept<PredecessorMap, vertex_descriptor> )); //store predecessor there
 757   BOOST_CONCEPT_ASSERT(( ReadWritePropertyMapConcept<ColorMap, vertex_descriptor> )); //write corresponding tree
 758   BOOST_CONCEPT_ASSERT(( ReadWritePropertyMapConcept<DistanceMap, vertex_descriptor> )); //write distance to source/sink
 759   BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<IndexMap, vertex_descriptor> )); //get index 0...|V|-1
 760   BOOST_ASSERT(num_vertices(g) >= 2 && src != sink);
 761
 762   detail::bk_max_flow<
 763     Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, ReverseEdgeMap,
 764     PredecessorMap, ColorMap, DistanceMap, IndexMap
 765   > algo(g, cap, res_cap, rev_map, pre_map, color, dist, idx, src, sink);
 766
 767   return algo.max_flow();
 768 }
 769
 770 /**
 771  * non-named-parameter version, given capacity, residucal_capacity,
 772  * reverse_edges, and an index map.
 773  */
 774 template<class Graph,
 775          class CapacityEdgeMap,
 776          class ResidualCapacityEdgeMap,
 777          class ReverseEdgeMap,
 778          class IndexMap>
 779 typename property_traits<CapacityEdgeMap>::value_type
 780 boykov_kolmogorov_max_flow(Graph& g,
 781                            CapacityEdgeMap cap,
 782                            ResidualCapacityEdgeMap res_cap,
 783                            ReverseEdgeMap rev,
 784                            IndexMap idx,
 785                            typename graph_traits<Graph>::vertex_descriptor src,
 786                            typename graph_traits<Graph>::vertex_descriptor sink)
 787 {
 788   typename graph_traits<Graph>::vertices_size_type n_verts = num_vertices(g);
 789   std::vector<typename graph_traits<Graph>::edge_descriptor> predecessor_vec(n_verts);
 790   std::vector<default_color_type> color_vec(n_verts);
 791   std::vector<typename graph_traits<Graph>::vertices_size_type> distance_vec(n_verts);
 792   return
 793     boykov_kolmogorov_max_flow(
 794       g, cap, res_cap, rev,
 795       make_iterator_property_map(predecessor_vec.begin(), idx),
 796       make_iterator_property_map(color_vec.begin(), idx),
 797       make_iterator_property_map(distance_vec.begin(), idx),
 798       idx, src, sink);
 799 }
 800
 801 /**
 802  * non-named-parameter version, some given: capacity, residual_capacity,
 803  * reverse_edges, color_map and an index map. Use this if you are interested in
 804  * the minimum cut, as the color map provides that info.
 805  */
 806 template<class Graph,
 807          class CapacityEdgeMap,
 808          class ResidualCapacityEdgeMap,
 809          class ReverseEdgeMap,
 810          class ColorMap,
 811          class IndexMap>
 812 typename property_traits<CapacityEdgeMap>::value_type
 813 boykov_kolmogorov_max_flow(Graph& g,
 814                            CapacityEdgeMap cap,
 815                            ResidualCapacityEdgeMap res_cap,
 816                            ReverseEdgeMap rev,
 817                            ColorMap color,
 818                            IndexMap idx,
 819                            typename graph_traits<Graph>::vertex_descriptor src,
 820                            typename graph_traits<Graph>::vertex_descriptor sink)
 821 {
 822   typename graph_traits<Graph>::vertices_size_type n_verts = num_vertices(g);
 823   std::vector<typename graph_traits<Graph>::edge_descriptor> predecessor_vec(n_verts);
 824   std::vector<typename graph_traits<Graph>::vertices_size_type> distance_vec(n_verts);
 825   return
 826     boykov_kolmogorov_max_flow(
 827       g, cap, res_cap, rev,
 828       make_iterator_property_map(predecessor_vec.begin(), idx),
 829       color,
 830       make_iterator_property_map(distance_vec.begin(), idx),
 831       idx, src, sink);
 832 }
 833
 834 /**
 835  * named-parameter version, some given
 836  */
 837 template<class Graph, class P, class T, class R>
 838 typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type
 839 boykov_kolmogorov_max_flow(Graph& g,
 840                            typename graph_traits<Graph>::vertex_descriptor src,
 841                            typename graph_traits<Graph>::vertex_descriptor sink,
 842                            const bgl_named_params<P, T, R>& params)
 843 {
 844   return
 845   boykov_kolmogorov_max_flow(
 846     g,
 847     choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),
 848     choose_pmap(get_param(params, edge_residual_capacity), g, edge_residual_capacity),
 849     choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),
 850     choose_pmap(get_param(params, vertex_predecessor), g, vertex_predecessor),
 851     choose_pmap(get_param(params, vertex_color), g, vertex_color),
 852     choose_pmap(get_param(params, vertex_distance), g, vertex_distance),
 853     choose_const_pmap(get_param(params, vertex_index), g, vertex_index),
 854     src, sink);
 855 }
 856
 857 /**
 858  * named-parameter version, none given
 859  */
 860 template<class Graph>
 861 typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type
 862 boykov_kolmogorov_max_flow(Graph& g,
 863                            typename graph_traits<Graph>::vertex_descriptor src,
 864                            typename graph_traits<Graph>::vertex_descriptor sink)
 865 {
 866   bgl_named_params<int, buffer_param_t> params(0); // bogus empty param
 867   return boykov_kolmogorov_max_flow(g, src, sink, params);
 868 }
 869
 870 } // namespace boost
 871
 872 #endif // BOOST_BOYKOV_KOLMOGOROV_MAX_FLOW_HPP
 873