RhodeCode

\ref Circulation is a preflow push-relabel algorithm implemented directly 

for finding feasible circulations, which is a somewhat different problem,

but it is strongly related to maximum flow.

For more information, see \ref Circulation.

*/

/**

@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms

@ingroup algs

\brief Algorithms for finding minimum cost flows and circulations.

This group contains the algorithms for finding minimum cost flows and

circulations. For more information about this problem and its dual

solution see \ref min_cost_flow "Minimum Cost Flow Problem".

LEMON contains several algorithms for this problem.

 - \ref NetworkSimplex Primal Network Simplex algorithm with various

   pivot strategies.

 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on

   cost scaling.

 - \ref CapacityScaling Successive Shortest %Path algorithm with optional

   capacity scaling.

 - \ref CancelAndTighten The Cancel and Tighten algorithm.

 - \ref CycleCanceling Cycle-Canceling algorithms.

In general NetworkSimplex is the most efficient implementation,

but in special cases other algorithms could be faster.

For example, if the total supply and/or capacities are rather small,

CapacityScaling is usually the fastest algorithm (without effective scaling).

*/

/**

@defgroup min_cut Minimum Cut Algorithms

@ingroup algs

\brief Algorithms for finding minimum cut in graphs.

This group contains the algorithms for finding minimum cut in graphs.

The \e minimum \e cut \e problem is to find a non-empty and non-complete

\f$X\f$ subset of the nodes with minimum overall capacity on

outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a

\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum

cut is the \f$X\f$ solution of the next optimization problem:

\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}

LEMON contains several algorithms related to minimum cut problems:

- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut

  in directed graphs.

- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for

  calculating minimum cut in undirected graphs.

- \ref GomoryHu "Gomory-Hu tree computation" for calculating

  all-pairs minimum cut in undirected graphs.

If you want to find minimum cut just between two distinict nodes,

see the \ref max_flow "maximum flow problem".

*/

/**

@defgroup graph_properties Connectivity and Other Graph Properties

@ingroup algs

\brief Algorithms for discovering the graph properties

This group contains the algorithms for discovering the graph properties

like connectivity, bipartiteness, euler property, simplicity etc.

*/

/**

@defgroup planar Planarity Embedding and Drawing

@ingroup algs

\brief Algorithms for planarity checking, embedding and drawing

This group contains the algorithms for planarity checking,

embedding and drawing.

\image html planar.png

\image latex planar.eps "Plane graph" width=\textwidth

*/

/**

@defgroup matching Matching Algorithms

@ingroup algs

\brief Algorithms for finding matchings in graphs and bipartite graphs.

This group contains the algorithms for calculating

matchings in graphs and bipartite graphs. The general matching problem is

finding a subset of the edges for which each node has at most one incident

edge.

There are several different algorithms for calculate matchings in

graphs.  The matching problems in bipartite graphs are generally

easier than in general graphs. The goal of the matching optimization

can be finding maximum cardinality, maximum weight or minimum cost

matching. The search can be constrained to find perfect or

maximum cardinality matching.

The matching algorithms implemented in LEMON:

- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm

  for calculating maximum cardinality matching in bipartite graphs.

- \ref PrBipartiteMatching Push-relabel algorithm

  for calculating maximum cardinality matching in bipartite graphs.

- \ref MaxWeightedBipartiteMatching

  Successive shortest path algorithm for calculating maximum weighted

  matching and maximum weighted bipartite matching in bipartite graphs.

- \ref MinCostMaxBipartiteMatching

  Successive shortest path algorithm for calculating minimum cost maximum

  matching in bipartite graphs.

- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating

  maximum cardinality matching in general graphs.

- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating

  maximum weighted matching in general graphs.

- \ref MaxWeightedPerfectMatching

  Edmond's blossom shrinking algorithm for calculating maximum weighted

        local_reached=false;

      }

      _reached = &m;

      return *this;

    }

    ///Sets the map that indicates which nodes are processed.

    ///Sets the map that indicates which nodes are processed.

    ///If you don't use this function before calling \ref run(Node) "run()"

    ///or \ref init(), an instance will be allocated automatically.

    ///The destructor deallocates this automatically allocated map,

    ///of course.

    ///\return <tt> (*this) </tt>

    Bfs &processedMap(ProcessedMap &m)

    {

      if(local_processed) {

        delete _processed;

        local_processed=false;

      }

      _processed = &m;

      return *this;

    }

    ///Sets the map that stores the distances of the nodes.

    ///Sets the map that stores the distances of the nodes calculated by

    ///the algorithm.

    ///If you don't use this function before calling \ref run(Node) "run()"

    ///or \ref init(), an instance will be allocated automatically.

    ///The destructor deallocates this automatically allocated map,

    ///of course.

    ///\return <tt> (*this) </tt>

    Bfs &distMap(DistMap &m)

    {

      if(local_dist) {

        delete _dist;

        local_dist=false;

      }

      _dist = &m;

      return *this;

    }

  public:

    ///\name Execution Control

    ///The simplest way to execute the BFS algorithm is to use one of the

    ///member functions called \ref run(Node) "run()".\n

    ///\ref addSource(). Finally the actual path computation can be

    ///performed with one of the \ref start() functions.

    ///@{

    ///\brief Initializes the internal data structures.

    ///

    ///Initializes the internal data structures.

    void init()

    {

      create_maps();

      _queue.resize(countNodes(*G));

      _queue_head=_queue_tail=0;

      _curr_dist=1;

      for ( NodeIt u(*G) ; u!=INVALID ; ++u ) {

        _pred->set(u,INVALID);

        _reached->set(u,false);

        _processed->set(u,false);

      }

    }

    ///Adds a new source node.

    ///Adds a new source node to the set of nodes to be processed.

    ///

    void addSource(Node s)

    {

      if(!(*_reached)[s])

        {

          _reached->set(s,true);

          _pred->set(s,INVALID);

          _dist->set(s,0);

          _queue[_queue_head++]=s;

          _queue_next_dist=_queue_head;

        }

    }

    ///Processes the next node.

    ///Processes the next node.

    ///

    ///\return The processed node.

    ///

    ///\pre The queue must not be empty.

    Node processNextNode()

    {

      if(_queue_tail==_queue_next_dist) {

        _curr_dist++;

    ///

    /// \ref named-templ-param "Named parameter" for setting ReachedMap type.

    template <class T>

    struct SetReachedMap : public BfsVisit< Digraph, Visitor,

                                            SetReachedMapTraits<T> > {

      typedef BfsVisit< Digraph, Visitor, SetReachedMapTraits<T> > Create;

    };

    ///@}

  public:

    /// \brief Constructor.

    ///

    /// Constructor.

    ///

    /// \param digraph The digraph the algorithm runs on.

    /// \param visitor The visitor object of the algorithm.

    BfsVisit(const Digraph& digraph, Visitor& visitor)

      : _digraph(&digraph), _visitor(&visitor),

        _reached(0), local_reached(false) {}

    /// \brief Destructor.

    ~BfsVisit() {

      if(local_reached) delete _reached;

    }

    /// \brief Sets the map that indicates which nodes are reached.

    ///

    /// Sets the map that indicates which nodes are reached.

    /// If you don't use this function before calling \ref run(Node) "run()"

    /// or \ref init(), an instance will be allocated automatically.

    /// The destructor deallocates this automatically allocated map,

    /// of course.

    /// \return <tt> (*this) </tt>

    BfsVisit &reachedMap(ReachedMap &m) {

      if(local_reached) {

        delete _reached;

        local_reached = false;

      }

      _reached = &m;

      return *this;

    }

  public:

    /// \name Execution Control

    /// The simplest way to execute the BFS algorithm is to use one of the

    /// member functions called \ref run(Node) "run()".\n

    /// \ref addSource(). Finally the actual path computation can be

    /// performed with one of the \ref start() functions.

    /// @{

    /// \brief Initializes the internal data structures.

    ///

    /// Initializes the internal data structures.

    void init() {

      create_maps();

      _list.resize(countNodes(*_digraph));

      _list_front = _list_back = -1;

      for (NodeIt u(*_digraph) ; u != INVALID ; ++u) {

        _reached->set(u, false);

      }

    }

    /// \brief Adds a new source node.

    ///

    /// Adds a new source node to the set of nodes to be processed.

    void addSource(Node s) {

      if(!(*_reached)[s]) {

          _reached->set(s,true);

          _visitor->start(s);

          _visitor->reach(s);

          _list[++_list_back] = s;

        }

    }

    /// \brief Processes the next node.

    ///

    /// Processes the next node.

    ///

    /// \return The processed node.

    ///

    /// \pre The queue must not be empty.

    Node processNextNode() {

      Node n = _list[++_list_front];

      _visitor->process(n);

      Arc e;

      for (_digraph->firstOut(e, n); e != INVALID; _digraph->nextOut(e)) {

        Node m = _digraph->target(e);

        if (!(*_reached)[m]) {

          _visitor->discover(e);

          _visitor->reach(m);

          _reached->set(m, true);

          _list[++_list_back] = m;

        } else {

///\file

///\brief Push-relabel algorithm for finding a feasible circulation.

///

namespace lemon {

  /// \brief Default traits class of Circulation class.

  ///

  /// Default traits class of Circulation class.

  ///

  /// \tparam GR Type of the digraph the algorithm runs on.

  /// \tparam LM The type of the lower bound map.

  /// \tparam UM The type of the upper bound (capacity) map.

  /// \tparam SM The type of the supply map.

  template <typename GR, typename LM,

            typename UM, typename SM>

  struct CirculationDefaultTraits {

    /// \brief The type of the digraph the algorithm runs on.

    typedef GR Digraph;

    /// \brief The type of the lower bound map.

    ///

    /// The type of the map that stores the lower bounds on the arcs.

    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.

    typedef LM LowerMap;

    /// \brief The type of the upper bound (capacity) map.

    ///

    /// The type of the map that stores the upper bounds (capacities)

    /// on the arcs.

    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.

    typedef UM UpperMap;

    /// \brief The type of supply map.

    ///

    /// The type of the map that stores the signed supply values of the 

    /// nodes. 

    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.

    typedef SM SupplyMap;

    /// \brief The type of the flow and supply values.

    typedef typename SupplyMap::Value Value;

    /// \brief The type of the map that stores the flow values.

    ///

    /// The type of the map that stores the flow values.

    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"

    /// concept.

    typedef typename Digraph::template ArcMap<Value> FlowMap;

    /// \brief Instantiates a FlowMap.

    ///

    /// This function instantiates a \ref FlowMap.

    /// \param digraph The digraph for which we would like to define

    /// the flow map.

    static FlowMap* createFlowMap(const Digraph& digraph) {

      return new FlowMap(digraph);

    }

    /// \brief The elevator type used by the algorithm.

    ///

    /// The elevator type used by the algorithm.

    ///

    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;

    /// \brief Instantiates an Elevator.

    ///

    /// This function instantiates an \ref Elevator.

    /// \param digraph The digraph for which we would like to define

    /// the elevator.

    /// \param max_level The maximum level of the elevator.

    static Elevator* createElevator(const Digraph& digraph, int max_level) {

      return new Elevator(digraph, max_level);

    }

    /// \brief The tolerance used by the algorithm

    ///

    /// The tolerance used by the algorithm to handle inexact computation.

    typedef lemon::Tolerance<Value> Tolerance;

  };

  /**

     \brief Push-relabel algorithm for the network circulation problem.

     \ingroup max_flow

     This class implements a push-relabel algorithm for the \e network

     \e circulation problem.

     It is to find a feasible circulation when lower and upper bounds

     are given for the flow values on the arcs and lower bounds are

     given for the difference between the outgoing and incoming flow

     at the nodes.

     The exact formulation of this problem is the following.

     Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$

     \f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and

     upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$

     holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$

     denotes the signed supply values of the nodes.

     If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$

     supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with

     \f$-sup(u)\f$ demand.

     A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$

     solution of the following problem.

     \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)

     \geq sup(u) \quad \forall u\in V, \f]

     \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]

     The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be

     zero or negative in order to have a feasible solution (since the sum

     of the expressions on the left-hand side of the inequalities is zero).

      _flow = ↦

      return *this;

    }

    /// \brief Sets the elevator used by algorithm.

    ///

    /// Sets the elevator used by algorithm.

    /// If you don't use this function before calling \ref run() or

    /// \ref init(), an instance will be allocated automatically.

    /// The destructor deallocates this automatically allocated elevator,

    /// of course.

    /// \return <tt>(*this)</tt>

    Circulation& elevator(Elevator& elevator) {

      if (_local_level) {

        delete _level;

        _local_level = false;

      }

      _level = &elevator;

      return *this;

    }

    /// \brief Returns a const reference to the elevator.

    ///

    /// Returns a const reference to the elevator.

    ///

    /// \pre Either \ref run() or \ref init() must be called before

    /// using this function.

    const Elevator& elevator() const {

      return *_level;

    }

    /// \brief Sets the tolerance used by algorithm.

    ///

    /// Sets the tolerance used by algorithm.

    Circulation& tolerance(const Tolerance& tolerance) const {

      _tol = tolerance;

      return *this;

    }

    /// \brief Returns a const reference to the tolerance.

    ///

    /// Returns a const reference to the tolerance.

    const Tolerance& tolerance() const {

      return tolerance;

    }

    /// \name Execution Control

    /// The simplest way to execute the algorithm is to call \ref run().\n

    /// the \ref start() function.

    ///@{

    /// Initializes the internal data structures.

    /// Initializes the internal data structures and sets all flow values

    /// to the lower bound.

    void init()

    {

      LEMON_DEBUG(checkBoundMaps(),

        "Upper bounds must be greater or equal to the lower bounds");

      createStructures();

      for(NodeIt n(_g);n!=INVALID;++n) {

        (*_excess)[n] = (*_supply)[n];

      }

      for (ArcIt e(_g);e!=INVALID;++e) {

        _flow->set(e, (*_lo)[e]);

        (*_excess)[_g.target(e)] += (*_flow)[e];

        (*_excess)[_g.source(e)] -= (*_flow)[e];

      }

      // global relabeling tested, but in general case it provides

      // worse performance for random digraphs

      _level->initStart();

      for(NodeIt n(_g);n!=INVALID;++n)

        _level->initAddItem(n);

      _level->initFinish();

      for(NodeIt n(_g);n!=INVALID;++n)

        if(_tol.positive((*_excess)[n]))

          _level->activate(n);

    }

    /// Initializes the internal data structures using a greedy approach.

    /// Initializes the internal data structures using a greedy approach

    /// to construct the initial solution.

    void greedyInit()

    {

      LEMON_DEBUG(checkBoundMaps(),

        "Upper bounds must be greater or equal to the lower bounds");

      createStructures();

      for(NodeIt n(_g);n!=INVALID;++n) {

        local_reached=false;

      }

      _reached = &m;

      return *this;

    }

    ///Sets the map that indicates which nodes are processed.

    ///Sets the map that indicates which nodes are processed.

    ///If you don't use this function before calling \ref run(Node) "run()"

    ///or \ref init(), an instance will be allocated automatically.

    ///The destructor deallocates this automatically allocated map,

    ///of course.

    ///\return <tt> (*this) </tt>

    Dfs &processedMap(ProcessedMap &m)

    {

      if(local_processed) {

        delete _processed;

        local_processed=false;

      }

      _processed = &m;

      return *this;

    }

    ///Sets the map that stores the distances of the nodes.

    ///Sets the map that stores the distances of the nodes calculated by

    ///the algorithm.

    ///If you don't use this function before calling \ref run(Node) "run()"

    ///or \ref init(), an instance will be allocated automatically.

    ///The destructor deallocates this automatically allocated map,

    ///of course.

    ///\return <tt> (*this) </tt>

    Dfs &distMap(DistMap &m)

    {

      if(local_dist) {

        delete _dist;

        local_dist=false;

      }

      _dist = &m;

      return *this;

    }

  public:

    ///\name Execution Control

    ///The simplest way to execute the DFS algorithm is to use one of the

    ///member functions called \ref run(Node) "run()".\n

    ///and perform the actual computation with \ref start().

    ///This procedure can be repeated if there are nodes that have not

    ///been reached.

    ///@{

    ///\brief Initializes the internal data structures.

    ///

    ///Initializes the internal data structures.

    void init()

    {

      create_maps();

      _stack.resize(countNodes(*G));

      _stack_head=-1;

      for ( NodeIt u(*G) ; u!=INVALID ; ++u ) {

        _pred->set(u,INVALID);

        _reached->set(u,false);

        _processed->set(u,false);

      }

    }

    ///Adds a new source node.

    ///Adds a new source node to the set of nodes to be processed.

    ///

    ///\pre The stack must be empty. Otherwise the algorithm gives

    ///wrong results. (One of the outgoing arcs of all the source nodes

    ///except for the last one will not be visited and distances will

    ///also be wrong.)

    void addSource(Node s)

    {

      LEMON_DEBUG(emptyQueue(), "The stack is not empty.");

      if(!(*_reached)[s])

        {

          _reached->set(s,true);

          _pred->set(s,INVALID);

          OutArcIt e(*G,s);

          if(e!=INVALID) {

            _stack[++_stack_head]=e;

            _dist->set(s,_stack_head);

          }

          else {

            _processed->set(s,true);

            _dist->set(s,0);

          }

        }

    }

    ///

    /// \ref named-templ-param "Named parameter" for setting ReachedMap type.

    template <class T>

    struct SetReachedMap : public DfsVisit< Digraph, Visitor,

                                            SetReachedMapTraits<T> > {

      typedef DfsVisit< Digraph, Visitor, SetReachedMapTraits<T> > Create;

    };

    ///@}

  public:

    /// \brief Constructor.

    ///

    /// Constructor.

    ///

    /// \param digraph The digraph the algorithm runs on.

    /// \param visitor The visitor object of the algorithm.

    DfsVisit(const Digraph& digraph, Visitor& visitor)

      : _digraph(&digraph), _visitor(&visitor),

        _reached(0), local_reached(false) {}

    /// \brief Destructor.

    ~DfsVisit() {

      if(local_reached) delete _reached;

    }

    /// \brief Sets the map that indicates which nodes are reached.

    ///

    /// Sets the map that indicates which nodes are reached.

    /// If you don't use this function before calling \ref run(Node) "run()"

    /// or \ref init(), an instance will be allocated automatically.

    /// The destructor deallocates this automatically allocated map,

    /// of course.

    /// \return <tt> (*this) </tt>

    DfsVisit &reachedMap(ReachedMap &m) {

      if(local_reached) {

        delete _reached;

        local_reached=false;

      }

      _reached = &m;

      return *this;

    }

  public:

    /// \name Execution Control

    /// The simplest way to execute the DFS algorithm is to use one of the

    /// member functions called \ref run(Node) "run()".\n

    /// and perform the actual computation with \ref start().

    /// This procedure can be repeated if there are nodes that have not

    /// been reached.

    /// @{

    /// \brief Initializes the internal data structures.

    ///

    /// Initializes the internal data structures.

    void init() {

      create_maps();

      _stack.resize(countNodes(*_digraph));

      _stack_head = -1;

      for (NodeIt u(*_digraph) ; u != INVALID ; ++u) {

        _reached->set(u, false);

      }

    }

    /// \brief Adds a new source node.

    ///

    /// Adds a new source node to the set of nodes to be processed.

    ///

    /// \pre The stack must be empty. Otherwise the algorithm gives

    /// wrong results. (One of the outgoing arcs of all the source nodes

    /// except for the last one will not be visited and distances will

    /// also be wrong.)

    void addSource(Node s)

    {

      LEMON_DEBUG(emptyQueue(), "The stack is not empty.");

      if(!(*_reached)[s]) {

          _reached->set(s,true);

          _visitor->start(s);

          _visitor->reach(s);

          Arc e;

          _digraph->firstOut(e, s);

          if (e != INVALID) {

            _stack[++_stack_head] = e;

          } else {

            _visitor->leave(s);

            _visitor->stop(s);

          }

        }

    }

    /// \brief Processes the next arc.

    ///

    /// Processes the next arc.

    ///

    ///\return <tt> (*this) </tt>

    Dijkstra &distMap(DistMap &m)

    {

      if(local_dist) {

        delete _dist;

        local_dist=false;

      }

      _dist = &m;

      return *this;

    }

    ///Sets the heap and the cross reference used by algorithm.

    ///Sets the heap and the cross reference used by algorithm.

    ///If you don't use this function before calling \ref run(Node) "run()"

    ///or \ref init(), heap and cross reference instances will be

    ///allocated automatically.

    ///The destructor deallocates these automatically allocated objects,

    ///of course.

    ///\return <tt> (*this) </tt>

    Dijkstra &heap(Heap& hp, HeapCrossRef &cr)

    {

      if(local_heap_cross_ref) {

        delete _heap_cross_ref;

        local_heap_cross_ref=false;

      }

      _heap_cross_ref = &cr;

      if(local_heap) {

        delete _heap;

        local_heap=false;

      }

      _heap = &hp;

      return *this;

    }

  private:

    void finalizeNodeData(Node v,Value dst)

    {

      _processed->set(v,true);

      _dist->set(v, dst);

    }

  public:

    ///\name Execution Control

    ///The simplest way to execute the %Dijkstra algorithm is to use

    ///one of the member functions called \ref run(Node) "run()".\n

    ///\ref addSource(). Finally the actual path computation can be

    ///performed with one of the \ref start() functions.

    ///@{

    ///\brief Initializes the internal data structures.

    ///

    ///Initializes the internal data structures.

    void init()

    {

      create_maps();

      _heap->clear();

      for ( NodeIt u(*G) ; u!=INVALID ; ++u ) {

        _pred->set(u,INVALID);

        _processed->set(u,false);

        _heap_cross_ref->set(u,Heap::PRE_HEAP);

      }

    }

    ///Adds a new source node.

    ///Adds a new source node to the priority heap.

    ///The optional second parameter is the initial distance of the node.

    ///

    ///The function checks if the node has already been added to the heap and

    ///it is pushed to the heap only if either it was not in the heap

    ///or the shortest path found till then is shorter than \c dst.

    void addSource(Node s,Value dst=OperationTraits::zero())

    {

      if(_heap->state(s) != Heap::IN_HEAP) {

        _heap->push(s,dst);

      } else if(OperationTraits::less((*_heap)[s], dst)) {

        _heap->set(s,dst);

        _pred->set(s,INVALID);

      }

    }

    ///Processes the next node in the priority heap

    ///Processes the next node in the priority heap.

    ///

    ///\return The processed node.

    ///

    ///\warning The priority heap must not be empty.

    Node processNextNode()

    {

      Node v=_heap->top();

      Value oldvalue=_heap->prio();

	  }

	  tn = (*_pred)[tn];

	} else {

	  if ((*_weight)[sn] <= value) {

	    rn = sn;

            s_root = true;

	    value = (*_weight)[sn];

	  }

	  sn = (*_pred)[sn];

	}

      }

      typename Graph::template NodeMap<bool> reached(_graph, false);

      reached[_root] = true;

      cutMap.set(_root, !s_root);

      reached[rn] = true;

      cutMap.set(rn, s_root);

      std::vector<Node> st;

      for (NodeIt n(_graph); n != INVALID; ++n) {

	st.clear();

        Node nn = n;

	while (!reached[nn]) {

	  st.push_back(nn);

	  nn = (*_pred)[nn];

	}

	while (!st.empty()) {

	  cutMap.set(st.back(), cutMap[nn]);

	  st.pop_back();

	}

      }

      return value;

    }

    ///@}

    friend class MinCutNodeIt;

    /// Iterate on the nodes of a minimum cut

    /// This iterator class lists the nodes of a minimum cut found by

    /// GomoryHu. Before using it, you must allocate a GomoryHu class

    /// and call its \ref GomoryHu::run() "run()" method.

    ///

    /// This example counts the nodes in the minimum cut separating \c s from

    /// \c t.

    /// \code

    /// gom.run();

    /// int cnt=0;

    /// \endcode

    class MinCutNodeIt

    {

      bool _side;

      typename Graph::NodeIt _node_it;

      typename Graph::template NodeMap<bool> _cut;

    public:

      /// Constructor

      /// Constructor.

      ///

      MinCutNodeIt(GomoryHu const &gomory,

                   ///< The GomoryHu class. You must call its

                   ///  run() method

                   ///  before initializing this iterator.

                   const Node& s, ///< The base node.

                   const Node& t,

                   ///< The node you want to separate from node \c s.

                   bool side=true

                   ///< If it is \c true (default) then the iterator lists

                   ///  the nodes of the component containing \c s,

                   ///  otherwise it lists the other component.

                   /// \note As the minimum cut is not always unique,

                   /// \code

                   /// MinCutNodeIt(gomory, s, t, true);

                   /// \endcode

                   /// and

                   /// \code

                   /// MinCutNodeIt(gomory, t, s, false);

                   /// \endcode

                   /// does not necessarily give the same set of nodes.

                   /// However it is ensured that

                   /// \code

                   /// MinCutNodeIt(gomory, s, t, true);

                   /// \endcode

                   /// and

                   /// \code

                   /// MinCutNodeIt(gomory, s, t, false);

                   /// \endcode

                   /// together list each node exactly once.

                   )

        : _side(side), _cut(gomory._graph)

      {

        gomory.minCutMap(s,t,_cut);

        for(_node_it=typename Graph::NodeIt(gomory._graph);

            _node_it!=INVALID && _cut[_node_it]!=_side;

            ++_node_it) {}

      }

      /// Conversion to \c Node

      /// Conversion to \c Node.

      ///

      operator typename Graph::Node() const

      {

        return _node_it;

      }

      bool operator==(Invalid) { return _node_it==INVALID; }

      bool operator!=(Invalid) { return _node_it!=INVALID; }

      /// Next node

      /// Next node.

      ///

      MinCutNodeIt &operator++()

      {

        for(++_node_it;_node_it!=INVALID&&_cut[_node_it]!=_side;++_node_it) {}

        return *this;

      }

      /// Postfix incrementation

      /// Postfix incrementation.

      ///

      /// \warning This incrementation

      /// returns a \c Node, not a \c MinCutNodeIt, as one may

      /// expect.

      typename Graph::Node operator++(int)

      {

        typename Graph::Node n=*this;

        ++(*this);

        return n;

      }

    };

    friend class MinCutEdgeIt;

    /// Iterate on the edges of a minimum cut

    /// This iterator class lists the edges of a minimum cut found by

    /// GomoryHu. Before using it, you must allocate a GomoryHu class

    /// and call its \ref GomoryHu::run() "run()" method.

    ///

    /// This example computes the value of the minimum cut separating \c s from

    /// \c t.

    /// \code

    /// gom.run();

    /// int value=0;

    ///   value+=capacities[e];

    /// \endcode

    /// The result will be the same as the value returned by

    /// \ref GomoryHu::minCutValue() "gom.minCutValue(s,t)".

    class MinCutEdgeIt

    {

      bool _side;

      const Graph &_graph;

      typename Graph::NodeIt _node_it;

      typename Graph::OutArcIt _arc_it;

      typename Graph::template NodeMap<bool> _cut;

      void step()

      {

        ++_arc_it;

        while(_node_it!=INVALID && _arc_it==INVALID)

          {

            for(++_node_it;_node_it!=INVALID&&!_cut[_node_it];++_node_it) {}

            if(_node_it!=INVALID)

              _arc_it=typename Graph::OutArcIt(_graph,_node_it);

          }

      }

    public:

      /// Constructor

      /// Constructor.

      ///

      MinCutEdgeIt(GomoryHu const &gomory,

                   ///< The GomoryHu class. You must call its

                   ///  run() method

                   ///  before initializing this iterator.

                   const Node& s,  ///< The base node.

                   const Node& t,

                   ///< The node you want to separate from node \c s.

                   bool side=true

                   ///< If it is \c true (default) then the listed arcs

                   ///  will be oriented from the

                   ///  nodes of the component containing \c s,

                   ///  otherwise they will be oriented in the opposite

                   ///  direction.

                   )

        : _graph(gomory._graph), _cut(_graph)

      {

        gomory.minCutMap(s,t,_cut);

        if(!side)

          for(typename Graph::NodeIt n(_graph);n!=INVALID;++n)

            _cut[n]=!_cut[n];

    };

    /// @}

    /// \brief Constructor.

    ///

    /// \param digraph The digraph the algorithm will run on.

    /// \param cost The cost map used by the algorithm.

    MinCostArborescence(const Digraph& digraph, const CostMap& cost)

      : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),

        _arborescence(0), local_arborescence(false),

        _arc_order(0), _node_order(0), _cost_arcs(0),

        _heap_cross_ref(0), _heap(0) {}

    /// \brief Destructor.

    ~MinCostArborescence() {

      destroyStructures();

    }

    /// \brief Sets the arborescence map.

    ///

    /// Sets the arborescence map.

    /// \return <tt>(*this)</tt>

    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {

      if (local_arborescence) {

        delete _arborescence;

      }

      local_arborescence = false;

      _arborescence = &m;

      return *this;

    }