//  Copyright (c) 2006, Stephan Diederich

//

//  This code may be used under either of the following two licences:

//

//    Permission is hereby granted, free of charge, to any person

//    obtaining a copy of this software and associated documentation

//    files (the "Software"), to deal in the Software without

//    restriction, including without limitation the rights to use,

//    copy, modify, merge, publish, distribute, sublicense, and/or

//    sell copies of the Software, and to permit persons to whom the

//    Software is furnished to do so, subject to the following

//    conditions:

//

//    The above copyright notice and this permission notice shall be

//    included in all copies or substantial portions of the Software.

//

//    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,

//    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES

//    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND

//    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT

//    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,

//    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING

//    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

//    OTHER DEALINGS IN THE SOFTWARE. OF SUCH DAMAGE.

//

//  Or:

//

//    Distributed under the Boost Software License, Version 1.0.

//    (See accompanying file LICENSE_1_0.txt or copy at

//    http://www.boost.org/LICENSE_1_0.txt)



#ifndef BOOST_KOLMOGOROV_MAX_FLOW_HPP

#define BOOST_KOLMOGOROV_MAX_FLOW_HPP



#include <boost/config.hpp>

#include <cassert>

#include <vector>

#include <list>

#include <utility>

#include <iosfwd>

#include <algorithm> // for std::min and std::max



#include <boost/pending/queue.hpp>

#include <boost/limits.hpp>

#include <boost/property_map.hpp>

#include <boost/none_t.hpp>

#include <boost/graph/graph_concepts.hpp>

#include <boost/graph/named_function_params.hpp>



namespace boost {

  namespace detail {



    template <class Graph,

              class EdgeCapacityMap,

              class ResidualCapacityEdgeMap,

              class ReverseEdgeMap,

              class PredecessorMap,

              class ColorMap,

              class DistanceMap,

              class IndexMap>

    class kolmogorov{

      typedef typename property_traits<EdgeCapacityMap>::value_type tEdgeVal;

      typedef graph_traits<Graph> tGraphTraits;

      typedef typename tGraphTraits::vertex_iterator vertex_iterator;

      typedef typename tGraphTraits::vertex_descriptor vertex_descriptor;

      typedef typename tGraphTraits::edge_descriptor edge_descriptor;

      typedef typename tGraphTraits::edge_iterator edge_iterator;

      typedef typename tGraphTraits::out_edge_iterator out_edge_iterator;

      typedef boost::queue<vertex_descriptor> tQueue;                               //queue of vertices, used in adoption-stage

      typedef typename property_traits<ColorMap>::value_type tColorValue;

      typedef color_traits<tColorValue> tColorTraits;

      typedef typename property_traits<DistanceMap>::value_type tDistanceVal;



        public:

          kolmogorov(Graph& g,

                     EdgeCapacityMap cap,

                     ResidualCapacityEdgeMap res,

                     ReverseEdgeMap rev,

                     PredecessorMap pre,

                     ColorMap color,

                     DistanceMap dist,

                     IndexMap idx,

                     vertex_descriptor src,

                     vertex_descriptor sink):

          m_g(g),

          m_index_map(idx),

          m_cap_map(cap),

          m_res_cap_map(res),

          m_rev_edge_map(rev),

          m_pre_map(pre),

          m_tree_map(color),

          m_dist_map(dist),

          m_source(src),

          m_sink(sink),

          m_active_nodes(),

          m_in_active_list_vec(num_vertices(g), false),

          m_in_active_list_map(make_iterator_property_map(m_in_active_list_vec.begin(), m_index_map)),

          m_has_parent_vec(num_vertices(g), false),

          m_has_parent_map(make_iterator_property_map(m_has_parent_vec.begin(), m_index_map)),

          m_time_vec(num_vertices(g), 0),

          m_time_map(make_iterator_property_map(m_time_vec.begin(), m_index_map)),

          m_flow(0),

          m_time(1),

          m_last_grow_vertex(graph_traits<Graph>::null_vertex()){

            // initialize the color-map with gray-values

            vertex_iterator vi, v_end;

            for(tie(vi, v_end) = vertices(m_g); vi != v_end; ++vi){

              set_tree(*vi, tColorTraits::gray());

            }

            // Initialize flow to zero which means initializing

            // the residual capacity equal to the capacity

            edge_iterator ei, e_end;

            for(tie(ei, e_end) = edges(m_g); ei != e_end; ++ei) {

              m_res_cap_map[*ei] = m_cap_map[*ei];

              assert(m_rev_edge_map[m_rev_edge_map[*ei]] == *ei); //check if the reverse edge map is build up properly

            }

            //init the search trees with the two terminals

            set_tree(m_source, tColorTraits::black());

            set_tree(m_sink, tColorTraits::white());

            m_time_map[m_source] = 1;

            m_time_map[m_sink] = 1;

          }

          

          ~kolmogorov(){}



          tEdgeVal max_flow(){

            //augment direct paths from SOURCE->SINK and SOURCE->VERTEX->SINK

            augment_direct_paths();

            //start the main-loop

            while(true){

              bool path_found;

              edge_descriptor connecting_edge;

              tie(connecting_edge, path_found) = grow(); //find a path from source to sink

              if(!path_found){

                //we're finished, no more paths were found

                break; 

              }

              ++m_time;

              augment(connecting_edge); //augment that path

              adopt(); //rebuild search tree structure

            }

            return m_flow;

          }



          //the complete class is protected, as we want access to members in derived test-class (see $(BOOST_ROOT)/libs/graph/test/kolmogorov_max_flow_test.cpp)

        protected:

          void augment_direct_paths(){

            //in a first step, we augment all direct paths from source->NODE->sink

            //and additionally paths from source->sink

            //this improves especially graphcuts for segmentation, as most of the nodes have source/sink connects

            //but shouldn't have an impact on other maxflow problems (this is done in grow() anyway)

            out_edge_iterator ei, e_end;

            for(tie(ei, e_end) = out_edges(m_source, m_g); ei != e_end; ++ei){

              edge_descriptor from_source = *ei;

              vertex_descriptor current_node = target(from_source, m_g);

              if(current_node == m_sink){

                tEdgeVal cap = m_res_cap_map[from_source];

                m_res_cap_map[from_source] = 0;

                m_flow += cap;

                continue;

              }

              edge_descriptor to_sink;

              bool is_there;

              tie(to_sink, is_there) = edge(current_node, m_sink, m_g);

              if(is_there){

                tEdgeVal cap_from_source = m_res_cap_map[from_source];

                tEdgeVal cap_to_sink = m_res_cap_map[to_sink];

                if(cap_from_source > cap_to_sink){

                  set_tree(current_node, tColorTraits::black());

                  add_active_node(current_node);

                  set_edge_to_parent(current_node, from_source);

                  m_dist_map[current_node] = 1;

                  m_time_map[current_node] = 1;

                  //add stuff to flow and update residuals

                  //we dont need to update reverse_edges, as incoming/outgoing edges to/from source/sink don't count for max-flow

                  m_res_cap_map[from_source] -= cap_to_sink;

                  m_res_cap_map[to_sink] = 0;

                  m_flow += cap_to_sink;

                } else if(cap_to_sink > 0){

                  set_tree(current_node, tColorTraits::white());

                  add_active_node(current_node);

                  set_edge_to_parent(current_node, to_sink);

                  m_dist_map[current_node] = 1;

                  m_time_map[current_node] = 1;

                  //add stuff to flow and update residuals

                  //we dont need to update reverse_edges, as incoming/outgoing edges to/from source/sink don't count for max-flow

                  m_res_cap_map[to_sink] -= cap_from_source;

                  m_res_cap_map[from_source] = 0;

                  m_flow += cap_from_source;

                }

              } else if(m_res_cap_map[from_source]){

                //there is no sink connect, so we can't augment this path

                //but to avoid adding m_source to the active nodes, we just activate this node and set the approciate things

                set_tree(current_node, tColorTraits::black());

                set_edge_to_parent(current_node, from_source);

                m_dist_map[current_node] = 1;

                m_time_map[current_node] = 1;

                add_active_node(current_node);

              }

            }

            for(tie(ei, e_end) = out_edges(m_sink, m_g); ei != e_end; ++ei){

              edge_descriptor to_sink = m_rev_edge_map[*ei];

              vertex_descriptor current_node = source(to_sink, m_g);

              if(m_res_cap_map[to_sink]){

                set_tree(current_node, tColorTraits::white());

                set_edge_to_parent(current_node, to_sink);

                m_dist_map[current_node] = 1;

                m_time_map[current_node] = 1;

                add_active_node(current_node);

              }

            }

          }



          /**

          * returns a pair of an edge and a boolean. if the bool is true, the edge is a connection of a found path from s->t , read "the link" and 

          *   source(returnVal, m_g) is the end of the path found in the source-tree

          *   target(returnVal, m_g) is the beginning of the path found in the sink-tree

          */

          std::pair<edge_descriptor, bool> grow(){

            assert(m_orphans.empty());

            vertex_descriptor current_node;

            while((current_node = get_next_active_node()) != graph_traits<Graph>::null_vertex()){ //if there is one

              assert(get_tree(current_node) != tColorTraits::gray()  && (has_parent(current_node) || current_node==m_source || current_node==m_sink)); 

              if(get_tree(current_node) == tColorTraits::black()){ 

                //source tree growing

                out_edge_iterator ei, e_end;

                if(current_node != m_last_grow_vertex){

                  m_last_grow_vertex = current_node;

                  tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g);

                } 

                for(; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it){

                  edge_descriptor out_edge = *m_last_grow_edge_it;

                  if(m_res_cap_map[out_edge] > 0){ //check if we have capacity left on this edge

                    vertex_descriptor other_node = target(out_edge, m_g);

                    if(get_tree(other_node) == tColorTraits::gray()){ //it's a free node

                      set_tree(other_node, tColorTraits::black()); //aquire other node to our search tree

                      set_edge_to_parent(other_node, out_edge);   //set us as parent

                      m_dist_map[other_node] = m_dist_map[current_node] + 1;  //and update the distance-heuristic

                      m_time_map[other_node] = m_time_map[current_node];

                      add_active_node(other_node);

                    } else if(get_tree(other_node) == tColorTraits::black()){

                      if(is_closer_to_terminal(current_node, other_node)){ //we do this to get shorter paths. check if we are nearer to the source as its parent is

                        set_edge_to_parent(other_node, out_edge);

                        m_dist_map[other_node] = m_dist_map[current_node] + 1;

                        m_time_map[other_node] = m_time_map[current_node];

                      }

                    } else{ 

                      assert(get_tree(other_node)==tColorTraits::white());

                      //kewl, found a path from one to the other search tree, return the connecting edge in src->sink dir

                      return std::make_pair(out_edge, true);

                    }

                  }

                } //for all out-edges

              } //source-tree-growing

              else{

                assert(get_tree(current_node) == tColorTraits::white());

                out_edge_iterator ei, e_end;

                if(current_node != m_last_grow_vertex){

                  m_last_grow_vertex = current_node;

                  tie(m_last_grow_edge_it, m_last_grow_edge_end) = out_edges(current_node, m_g);

                }

                for(; m_last_grow_edge_it != m_last_grow_edge_end; ++m_last_grow_edge_it){

                  edge_descriptor in_edge = m_rev_edge_map[*m_last_grow_edge_it];

                  if(m_res_cap_map[in_edge] > 0){ //check if there is capacity left

                    vertex_descriptor other_node = source(in_edge, m_g);

                    if(get_tree(other_node) == tColorTraits::gray()){ //it's a free node

                      set_tree(other_node, tColorTraits::white());      //aquire that node to our search tree

                      set_edge_to_parent(other_node, in_edge);          //set us as parent

                      add_active_node(other_node);                      //activate that node

                      m_dist_map[other_node] = m_dist_map[current_node] + 1; //set its distance

                      m_time_map[other_node] = m_time_map[current_node];     //and time

                    } else if(get_tree(other_node) == tColorTraits::white()){

                      if(is_closer_to_terminal(current_node, other_node)){

                        //we are closer to the sink than its parent is, so we "adopt" him

                        set_edge_to_parent(other_node, in_edge);

                        m_dist_map[other_node] = m_dist_map[current_node] + 1;

                        m_time_map[other_node] = m_time_map[current_node];

                      }

                    } else{

                      assert(get_tree(other_node)==tColorTraits::black());

                      //kewl, found a path from one to the other search tree, return the connecting edge in src->sink dir

                      return std::make_pair(in_edge, true);

                    }

                  }

                } //for all out-edges

              } //sink-tree growing

              //all edges of that node are processed, and no more paths were found. so remove if from the front of the active queue

              finish_node(current_node);

            } //while active_nodes not empty

            return std::make_pair(edge_descriptor(), false); //no active nodes anymore and no path found, we're done

          }



          /**

          * augments path from s->t and updates residual graph

          * source(e, m_g) is the end of the path found in the source-tree

          * target(e, m_g) is the beginning of the path found in the sink-tree

          * this phase generates orphans on satured edges, if the attached verts are from different search-trees

          * orphans are ordered in distance to sink/source. first the farest from the source are front_inserted into the orphans list,

          * and after that the sink-tree-orphans are front_inserted. when going to adoption stage the orphans are popped_front, and so we process the nearest 

          * verts to the terminals first

          */

          void augment(edge_descriptor e){

            assert(get_tree(target(e, m_g)) == tColorTraits::white());

            assert(get_tree(source(e, m_g)) == tColorTraits::black());

            assert(m_orphans.empty());



            const tEdgeVal bottleneck = find_bottleneck(e);

            //now we push the found flow through the path

            //for each edge we saturate we have to look for the verts that belong to that edge, one of them becomes an orphans

            //now process the connecting edge

            m_res_cap_map[e] -= bottleneck;

            assert(m_res_cap_map[e] >= 0);

            m_res_cap_map[m_rev_edge_map[e]] += bottleneck;



            //now we follow the path back to the source

            vertex_descriptor current_node = source(e, m_g);

            while(current_node != m_source){

              edge_descriptor pred = get_edge_to_parent(current_node);

              m_res_cap_map[pred] -= bottleneck;

              assert(m_res_cap_map[pred] >= 0);

              m_res_cap_map[m_rev_edge_map[pred]] += bottleneck;

              if(m_res_cap_map[pred] == 0){

                set_no_parent(current_node);

                m_orphans.push_front(current_node);

              }

              current_node = source(pred, m_g);

            }

            //then go forward in the sink-tree

            current_node = target(e, m_g);

            while(current_node != m_sink){

              edge_descriptor pred = get_edge_to_parent(current_node);

              m_res_cap_map[pred] -= bottleneck;

              assert(m_res_cap_map[pred] >= 0);

              m_res_cap_map[m_rev_edge_map[pred]] += bottleneck;

              if(m_res_cap_map[pred] == 0){

                set_no_parent(current_node);

                m_orphans.push_front(current_node);

              }

              current_node = target(pred, m_g);

            }

            //and add it to the max-flow

            m_flow += bottleneck;

          }



          /**

           * returns the bottleneck of a s->t path (end_of_path is last vertex in source-tree, begin_of_path is first vertex in sink-tree)

           */		

          inline tEdgeVal find_bottleneck(edge_descriptor e){

            BOOST_USING_STD_MIN();

            tEdgeVal minimum_cap = m_res_cap_map[e];

            vertex_descriptor current_node = source(e, m_g);

            //first go back in the source tree

            while(current_node != m_source){

              edge_descriptor pred = get_edge_to_parent(current_node);

              minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(minimum_cap, m_res_cap_map[pred]);

              current_node = source(pred, m_g);

            }

            //then go forward in the sink-tree

            current_node = target(e, m_g);

            while(current_node != m_sink){

              edge_descriptor pred = get_edge_to_parent(current_node);

              minimum_cap = min BOOST_PREVENT_MACRO_SUBSTITUTION(minimum_cap, m_res_cap_map[pred]);

              current_node = target(pred, m_g);

            }

            return minimum_cap;

          }



          /**

          * rebuild search trees

          * empty the queue of orphans, and find new parents for them or just drop them from the search trees

          */

          void adopt(){

            while(!m_orphans.empty() || !m_child_orphans.empty()){

              vertex_descriptor current_node;

              if(m_child_orphans.empty()){

                //get the next orphan from the main-queue  and remove it

                current_node = m_orphans.front();

                m_orphans.pop_front();

              } else{

                current_node = m_child_orphans.front();

                m_child_orphans.pop();

              }

              if(get_tree(current_node) == tColorTraits::black()){ 

                //we're in the source-tree

                tDistanceVal min_distance = (std::numeric_limits<tDistanceVal>::max)();

                edge_descriptor new_parent_edge;

                out_edge_iterator ei, e_end;

                for(tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){

                  const edge_descriptor in_edge = m_rev_edge_map[*ei];

                  assert(target(in_edge, m_g) == current_node); //we should be the target of this edge

                  if(m_res_cap_map[in_edge] > 0){

                    vertex_descriptor other_node = source(in_edge, m_g);

                    if(get_tree(other_node) == tColorTraits::black() && has_source_connect(other_node)){

                      if(m_dist_map[other_node] < min_distance){

                        min_distance = m_dist_map[other_node];

                        new_parent_edge = in_edge;

                      }

                    }

                  }

                }

                if(min_distance != (std::numeric_limits<tDistanceVal>::max)()){

                  set_edge_to_parent(current_node, new_parent_edge);

                  m_dist_map[current_node] = min_distance + 1;

                  m_time_map[current_node] = m_time;

                } else{

                  m_time_map[current_node] = 0;

                  for(tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){

                    edge_descriptor in_edge = m_rev_edge_map[*ei];

                    vertex_descriptor other_node = source(in_edge, m_g);

                    if(get_tree(other_node) == tColorTraits::black() && has_parent(other_node)){

                      if(m_res_cap_map[in_edge] > 0){

                        add_active_node(other_node);

                      }

                      if(source(get_edge_to_parent(other_node), m_g) == current_node){

                        //we are the parent of that node

                        //it has to find a new parent, too

                        set_no_parent(other_node);

                        m_child_orphans.push(other_node);

                      }

                    }

                  }

                  set_tree(current_node, tColorTraits::gray());

                } //no parent found

              } //source-tree-adoption

              else{

                //now we should be in the sink-tree, check that...

                assert(get_tree(current_node) == tColorTraits::white());

                out_edge_iterator ei, e_end;

                edge_descriptor new_parent_edge;

                tDistanceVal min_distance = (std::numeric_limits<tDistanceVal>::max)();

                for(tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){

                  const edge_descriptor out_edge = *ei;

                  if(m_res_cap_map[out_edge] > 0){

                    const vertex_descriptor other_node = target(out_edge, m_g);

                    if(get_tree(other_node) == tColorTraits::white() && has_sink_connect(other_node))

                      if(m_dist_map[other_node] < min_distance){

                        min_distance = m_dist_map[other_node];

                        new_parent_edge = out_edge;

                      }

                  }

                }

                if(min_distance != (std::numeric_limits<tDistanceVal>::max)()){

                  set_edge_to_parent(current_node, new_parent_edge);

                  m_dist_map[current_node] = min_distance + 1;

                  m_time_map[current_node] = m_time;

                } else{

                  m_time_map[current_node] = 0;

                  for(tie(ei, e_end) = out_edges(current_node, m_g); ei != e_end; ++ei){

                    const edge_descriptor out_edge = *ei;

                    const vertex_descriptor other_node = target(out_edge, m_g);

                    if(get_tree(other_node) == tColorTraits::white() && has_parent(other_node)){ 

                      if(m_res_cap_map[out_edge] > 0){

                        add_active_node(other_node);

                      }

                      if(target(get_edge_to_parent(other_node), m_g) == current_node){

                        //we were it's parent, so it has to find a new one, too

                        set_no_parent(other_node);

                        m_child_orphans.push(other_node);

                      }

                    }

                  }

                  set_tree(current_node, tColorTraits::gray());

                } //no parent found

              } //sink-tree adoption

            } //while !orphans.empty()

          } //adopt



          /**

          * return next active vertex if there is one, otherwise a null_vertex

          */	

          inline vertex_descriptor get_next_active_node(){

            while(true){

              if(m_active_nodes.empty())

                return graph_traits<Graph>::null_vertex();

              vertex_descriptor v = m_active_nodes.front();



              if(!has_parent(v) && v != m_source && v != m_sink){ //if it has no parent, this node can't be active(if its not source or sink)

                m_active_nodes.pop();

                m_in_active_list_map[v] = false;

              } else{

                assert(get_tree(v) == tColorTraits::black() || get_tree(v) == tColorTraits::white());

                return v;

              }

            }

          }



          /**

          * adds v as an active vertex, but only if its not in the list already

          */		

          inline void add_active_node(vertex_descriptor v){

            assert(get_tree(v) != tColorTraits::gray());

            if(m_in_active_list_map[v]){

              return;

            } else{

              m_in_active_list_map[v] = true;

              m_active_nodes.push(v);

            }

          }



          /**

           * 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)

           */

          inline void finish_node(vertex_descriptor v){

            assert(m_active_nodes.front() == v);

            m_active_nodes.pop();

            m_in_active_list_map[v] = false;

            m_last_grow_vertex = graph_traits<Graph>::null_vertex();

          }



          /**

          * removes a vertex from the queue of active nodes (actually this does nothing, 

          * but checks if this node has no parent edge, as this is the criteria for beeing no more active) 

          */

          inline void remove_active_node(vertex_descriptor v){

            assert(!has_parent(v));

          }



          /**

          * returns the search tree of v; tColorValue::black() for source tree, white() for sink tree, gray() for no tree

          */

          inline tColorValue get_tree(vertex_descriptor v) const {

            return m_tree_map[v];

          }



          /**

          * sets search tree of v; tColorValue::black() for source tree, white() for sink tree, gray() for no tree

          */

          inline void set_tree(vertex_descriptor v, tColorValue t){

            m_tree_map[v] = t;

          }



          /**

           * returns edge to parent vertex of v;

           */

          inline edge_descriptor get_edge_to_parent(vertex_descriptor v) const{

            return m_pre_map[v];

          }



          /**

           * returns true if the edge stored in m_pre_map[v] is a valid entry

           */

          inline bool has_parent(vertex_descriptor v) const{

            return m_has_parent_map[v];  

          }



          /**

           * sets edge to parent vertex of v; 

          */		

          inline void set_edge_to_parent(vertex_descriptor v, edge_descriptor f_edge_to_parent){

            assert(m_res_cap_map[f_edge_to_parent] > 0);

            m_pre_map[v] = f_edge_to_parent;

            m_has_parent_map[v] = true;

          }



          /**

           * removes the edge to parent of v (this is done by invalidating the entry an additional map)

           */

          inline void set_no_parent(vertex_descriptor v){

            m_has_parent_map[v] = false;

          }



          /**

           * checks if vertex v has a connect to the sink-vertex (@var m_sink)

           * @param v the vertex which is checked

           * @return true if a path to the sink was found, false if not

           */

          inline bool has_sink_connect(vertex_descriptor v){

            tDistanceVal current_distance = 0;

            vertex_descriptor current_vertex = v;

            while(true){

              if(m_time_map[current_vertex] == m_time){

                //we found a node which was already checked this round. use it for distance calculations

                current_distance += m_dist_map[current_vertex];

                break;

              }

              if(current_vertex == m_sink){

                m_time_map[m_sink] = m_time;

                break; 

              }

              if(has_parent(current_vertex)){

                //it has a parent, so get it

                current_vertex = target(get_edge_to_parent(current_vertex), m_g);

                ++current_distance;

              } else{

                //no path found

                return false;

              }

            }

            current_vertex=v;

            while(m_time_map[current_vertex] != m_time){

              m_dist_map[current_vertex] = current_distance--;

              m_time_map[current_vertex] = m_time;

              current_vertex = target(get_edge_to_parent(current_vertex), m_g);

            }

            return true;

          }



          /**

           * checks if vertex v has a connect to the source-vertex (@var m_source)

           * @param v the vertex which is checked

           * @return true if a path to the source was found, false if not

           */

          inline bool has_source_connect(vertex_descriptor v){

            tDistanceVal current_distance = 0;

            vertex_descriptor current_vertex = v;

            while(true){

              if(m_time_map[current_vertex] == m_time){

                //we found a node which was already checked this round. use it for distance calculations

                current_distance += m_dist_map[current_vertex];

                break;

              }

              if(current_vertex == m_source){

                m_time_map[m_source] = m_time;

                break;

              }

              if(has_parent(current_vertex)){

                //it has a parent, so get it

                current_vertex = source(get_edge_to_parent(current_vertex), m_g);

                ++current_distance;

              } else{

                //no path found

                return false;

              }

            }

            current_vertex=v;

            while(m_time_map[current_vertex] != m_time){

                m_dist_map[current_vertex] = current_distance-- ;

                m_time_map[current_vertex] = m_time;

                current_vertex = source(get_edge_to_parent(current_vertex), m_g);

            }

            return true;

          }



          /**

          * returns true, if p is closer to a terminal than q

          */

          inline bool is_closer_to_terminal(vertex_descriptor p, vertex_descriptor q){

            //checks the timestamps first, to build no cycles, and after that the real distance

            return (m_time_map[q] <= m_time_map[p] && m_dist_map[q] > m_dist_map[p]+1);

          }



          ////////

          // member vars

          ////////

          Graph& m_g;

          IndexMap m_index_map;

          EdgeCapacityMap m_cap_map;

          ResidualCapacityEdgeMap m_res_cap_map;

          ReverseEdgeMap m_rev_edge_map;

          PredecessorMap m_pre_map; //stores paths found in the growth stage

          ColorMap m_tree_map; //maps each vertex into one of the two search tree or none (gray())

          DistanceMap m_dist_map; //stores distance to source/sink nodes

          vertex_descriptor m_source;

          vertex_descriptor m_sink;



          tQueue m_active_nodes;

          std::vector<bool> m_in_active_list_vec;

          iterator_property_map<std::vector<bool>::iterator, IndexMap> m_in_active_list_map;



          std::list<vertex_descriptor> m_orphans;

          tQueue m_child_orphans; // we use a second queuqe for child orphans, as they are FIFO processed



          std::vector<bool> m_has_parent_vec;

          iterator_property_map<std::vector<bool>::iterator, IndexMap> m_has_parent_map;



          std::vector<long> m_time_vec; //timestamp of each node, used for sink/source-path calculations

          iterator_property_map<std::vector<long>::iterator, IndexMap> m_time_map;

          tEdgeVal m_flow;

          long m_time;

          vertex_descriptor m_last_grow_vertex;

          out_edge_iterator m_last_grow_edge_it;

          out_edge_iterator m_last_grow_edge_end;

    };

  } //namespace detail

  

  /**

   * non-named-parameter version, given everything

   * this is the catch all version

   */			

  template <class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap, class ReverseEdgeMap, 

    class PredecessorMap, class ColorMap, class DistanceMap, class IndexMap>

  typename property_traits<CapacityEdgeMap>::value_type

  kolmogorov_max_flow

      (Graph& g,

       CapacityEdgeMap cap,

       ResidualCapacityEdgeMap res_cap,

       ReverseEdgeMap rev_map,

       PredecessorMap pre_map,

       ColorMap color,

       DistanceMap dist,

       IndexMap idx,

       typename graph_traits<Graph>::vertex_descriptor src,

       typename graph_traits<Graph>::vertex_descriptor sink

       )

  {

    typedef typename graph_traits<Graph>::vertex_descriptor vertex_descriptor;

    typedef typename graph_traits<Graph>::edge_descriptor edge_descriptor;

    //as this method is the last one before we instantiate the solver, we do the concept checks here

    function_requires<VertexListGraphConcept<Graph> >(); //to have vertices(), num_vertices(),

    function_requires<EdgeListGraphConcept<Graph> >(); //to have edges()

    function_requires<IncidenceGraphConcept<Graph> >(); //to have source(), target() and out_edges()

    function_requires<LvaluePropertyMapConcept<CapacityEdgeMap, edge_descriptor> >(); //read flow-values from edges

    function_requires<Mutable_LvaluePropertyMapConcept<ResidualCapacityEdgeMap, edge_descriptor> >(); //write flow-values to residuals

    function_requires<LvaluePropertyMapConcept<ReverseEdgeMap, edge_descriptor> >(); //read out reverse edges

    function_requires<Mutable_LvaluePropertyMapConcept<PredecessorMap, vertex_descriptor> >(); //store predecessor there

    function_requires<Mutable_LvaluePropertyMapConcept<ColorMap, vertex_descriptor> >(); //write corresponding tree

    function_requires<Mutable_LvaluePropertyMapConcept<DistanceMap, vertex_descriptor> >(); //write distance to source/sink

    function_requires<ReadablePropertyMapConcept<IndexMap, vertex_descriptor> >(); //get index 0...|V|-1

    assert(num_vertices(g) >= 2 && src != sink);

    detail::kolmogorov<Graph, CapacityEdgeMap, ResidualCapacityEdgeMap, ReverseEdgeMap, PredecessorMap, ColorMap, DistanceMap, IndexMap>

        algo(g, cap, res_cap, rev_map, pre_map, color, dist, idx, src, sink);

        return algo.max_flow();

  }



  /**

   * non-named-parameter version, given: capacity, residucal_capacity, reverse_edges, and an index map.

   */

  template <class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap, class ReverseEdgeMap, class IndexMap>

  typename property_traits<CapacityEdgeMap>::value_type

   kolmogorov_max_flow

       (Graph& g,

       CapacityEdgeMap cap,

       ResidualCapacityEdgeMap res_cap,

       ReverseEdgeMap rev,

       IndexMap idx,

       typename graph_traits<Graph>::vertex_descriptor src,

       typename graph_traits<Graph>::vertex_descriptor sink)

   {

     typename graph_traits<Graph>::vertices_size_type n_verts = num_vertices(g);

     std::vector<typename graph_traits<Graph>::edge_descriptor> predecessor_vec(n_verts);

     std::vector<default_color_type> color_vec(n_verts);

     std::vector<typename graph_traits<Graph>::vertices_size_type> distance_vec(n_verts);

     return kolmogorov_max_flow

         (g, cap, res_cap, rev,

          make_iterator_property_map(predecessor_vec.begin(), idx),

          make_iterator_property_map(color_vec.begin(), idx),

          make_iterator_property_map(distance_vec.begin(), idx),

          idx, src, sink);

   }

   

  /**

   * non-named-parameter version, some given: capacity, residual_capacity, reverse_edges, color_map and an index map.

   * Use this if you are interested in the minimum cut, as the color map provides that info

   */

   template <class Graph, class CapacityEdgeMap, class ResidualCapacityEdgeMap, class ReverseEdgeMap, class ColorMap, class IndexMap>

   typename property_traits<CapacityEdgeMap>::value_type

   kolmogorov_max_flow

       (Graph& g,

        CapacityEdgeMap cap,

        ResidualCapacityEdgeMap res_cap,

        ReverseEdgeMap rev,

        ColorMap color,

        IndexMap idx,

        typename graph_traits<Graph>::vertex_descriptor src,

        typename graph_traits<Graph>::vertex_descriptor sink)

   {

     typename graph_traits<Graph>::vertices_size_type n_verts = num_vertices(g);

     std::vector<typename graph_traits<Graph>::edge_descriptor> predecessor_vec(n_verts);

     std::vector<typename graph_traits<Graph>::vertices_size_type> distance_vec(n_verts);



     return kolmogorov_max_flow

         (g, cap, res_cap, rev,

          make_iterator_property_map(predecessor_vec.begin(), idx),

          color,

          make_iterator_property_map(distance_vec.begin(), idx),

          idx, src, sink);

   }

   

  /**

   * named-parameter version, some given

   */

   template <class Graph, class P, class T, class R>

   typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type

   kolmogorov_max_flow

       (Graph& g,

        typename graph_traits<Graph>::vertex_descriptor src,

        typename graph_traits<Graph>::vertex_descriptor sink,

        const bgl_named_params<P, T, R>& params)

   {

     return kolmogorov_max_flow(g,

                                choose_const_pmap(get_param(params, edge_capacity), g, edge_capacity),

                                choose_pmap(get_param(params, edge_residual_capacity), g, edge_residual_capacity),

                                choose_const_pmap(get_param(params, edge_reverse), g, edge_reverse),

                                choose_pmap(get_param(params, vertex_predecessor), g, vertex_predecessor),

                                choose_pmap(get_param(params, vertex_color), g, vertex_color),

                                choose_pmap(get_param(params, vertex_distance), g, vertex_distance),

                                choose_const_pmap(get_param(params, vertex_index), g, vertex_index),

                                src, sink);

   }

   

  /**

   * named-parameter version, none given

   */

   template <class Graph>

   typename property_traits<typename property_map<Graph, edge_capacity_t>::const_type>::value_type

   kolmogorov_max_flow

       (Graph& g,

        typename graph_traits<Graph>::vertex_descriptor src,

        typename graph_traits<Graph>::vertex_descriptor sink)

   {

     bgl_named_params<int, buffer_param_t> params(0); // bogus empty param

     return kolmogorov_max_flow(g, src, sink, params);

   }

} // namespace boost



#endif // BOOST_KOLMOGOROV_MAX_FLOW_HPP



