RhodeCode

    /// node, then the \c reach parameter will be set to \c true.

    ///

    /// \param target The target node.

    /// \retval reach Indicates if the target node is reached.

    /// It should be initially \c false.

    ///

    /// \return The processed node.

    ///

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

    Node processNextNode(Node target, bool& reach) {

      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;

          reach = reach || (target == m);

        } else {

          _visitor->examine(e);

        }

      }

      return n;

    }

    /// \brief Processes the next node.

    ///

    /// Processes the next node and checks if at least one of reached

    /// nodes has \c true value in the \c nm node map. If one node

    /// with \c true value is reachable from the processed node, then the

    /// \c rnode parameter will be set to the first of such nodes.

    ///

    /// \param nm A \c bool (or convertible) node map that indicates the

    /// possible targets.

    /// \retval rnode The reached target node.

    /// It should be initially \c INVALID.

    ///

    /// \return The processed node.

    ///

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

    template <typename NM>

    Node processNextNode(const NM& nm, Node& rnode) {

      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;

          if (nm[m] && rnode == INVALID) rnode = m;

        } else {

          _visitor->examine(e);

        }

      }

      return n;

    }

    /// \brief The next node to be processed.

    ///

    /// Returns the next node to be processed or \c INVALID if the queue

    /// is empty.

    Node nextNode() const {

      return _list_front != _list_back ? _list[_list_front + 1] : INVALID;

    }

    /// \brief Returns \c false if there are nodes

    /// to be processed.

    ///

    /// Returns \c false if there are nodes

    /// to be processed in the queue.

    bool emptyQueue() const { return _list_front == _list_back; }

    /// \brief Returns the number of the nodes to be processed.

    ///

    /// Returns the number of the nodes to be processed in the queue.

    int queueSize() const { return _list_back - _list_front; }

    /// \brief Executes the algorithm.

    ///

    /// Executes the algorithm.

    ///

    /// This method runs the %BFS algorithm from the root node(s)

    /// in order to compute the shortest path to each node.

    ///

    /// The algorithm computes

    /// - the shortest path tree (forest),

    /// - the distance of each node from the root(s).

    ///

    /// \pre init() must be called and at least one root node should be added

    /// with addSource() before using this function.

    ///

    /// \note <tt>b.start()</tt> is just a shortcut of the following code.

    /// \code

    ///   while ( !b.emptyQueue() ) {

    ///     b.processNextNode();

    ///   }

    /// \endcode

    void start() {

      while ( !emptyQueue() ) processNextNode();

    }

    /// \brief Executes the algorithm until the given target node is reached.

    ///

    /// Executes the algorithm until the given target node is reached.

    ///

    /// This method runs the %BFS algorithm from the root node(s)

    /// in order to compute the shortest path to \c t.

    ///

    /// The algorithm computes

    /// - the shortest path to \c t,

    /// - the distance of \c t from the root(s).

    ///

    /// \pre init() must be called and at least one root node should be

    /// added with addSource() before using this function.

    ///

    /// \note <tt>b.start(t)</tt> is just a shortcut of the following code.

    /// \code

    ///   bool reach = false;

    ///   while ( !b.emptyQueue() && !reach ) {

    ///     b.processNextNode(t, reach);

    ///   }

    /// \endcode

    void start(Node t) {

      bool reach = false;

      while ( !emptyQueue() && !reach ) processNextNode(t, reach);

    }

    /// \brief Executes the algorithm until a condition is met.

    ///

    /// Executes the algorithm until a condition is met.

    ///

    /// This method runs the %BFS algorithm from the root node(s) in

    /// order to compute the shortest path to a node \c v with

    /// <tt>nm[v]</tt> true, if such a node can be found.

    ///

    /// \param nm must be a bool (or convertible) node map. The

    /// algorithm will stop when it reaches a node \c v with

    /// <tt>nm[v]</tt> true.

    ///

    /// \return The reached node \c v with <tt>nm[v]</tt> true or

    /// \c INVALID if no such node was found.

    ///

    /// \pre init() must be called and at least one root node should be

    /// added with addSource() before using this function.

    ///

    /// \note <tt>b.start(nm)</tt> is just a shortcut of the following code.

    /// \code

    ///   Node rnode = INVALID;

    ///   while ( !b.emptyQueue() && rnode == INVALID ) {

    ///     b.processNextNode(nm, rnode);

    ///   }

    ///   return rnode;

    /// \endcode

    template <typename NM>

    Node start(const NM &nm) {

      Node rnode = INVALID;

      while ( !emptyQueue() && rnode == INVALID ) {

        processNextNode(nm, rnode);

      }

      return rnode;

    }

    /// \brief Runs the algorithm from the given source node.

    ///

    /// This method runs the %BFS algorithm from node \c s

    /// in order to compute the shortest path to each node.

    ///

    /// The algorithm computes

    /// - the shortest path tree,

    /// - the distance of each node from the root.

    ///

    /// \note <tt>b.run(s)</tt> is just a shortcut of the following code.

    ///\code

    ///   b.init();

    ///   b.addSource(s);

    ///   b.start();

    ///\endcode

    void run(Node s) {

      init();

      addSource(s);

      start();

    }

    /// \brief Finds the shortest path between \c s and \c t.

    ///

    /// This method runs the %BFS algorithm from node \c s

    /// in order to compute the shortest path to node \c t

    /// (it stops searching when \c t is processed).

    ///

    /// \return \c true if \c t is reachable form \c s.

    ///

    /// \note Apart from the return value, <tt>b.run(s,t)</tt> is just a

    /// shortcut of the following code.

    ///\code

    ///   b.init();

    ///   b.addSource(s);

    ///   b.start(t);

    ///\endcode

    bool run(Node s,Node t) {

      init();

      addSource(s);

      start(t);

      return reached(t);

    }

    /// \brief Runs the algorithm to visit all nodes in the digraph.

    ///

    /// This method runs the %BFS algorithm in order to

    /// compute the shortest path to each node.

    ///

    /// The algorithm computes

    /// - the shortest path tree (forest),

    /// - the distance of each node from the root(s).

    ///

    /// \note <tt>b.run(s)</tt> is just a shortcut of the following code.

    ///\code

    ///  b.init();

    ///  for (NodeIt n(gr); n != INVALID; ++n) {

    ///    if (!b.reached(n)) {

    ///      b.addSource(n);

    ///      b.start();

    ///    }

    ///  }

    ///\endcode

    void run() {

      init();

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

        if (!reached(it)) {

          addSource(it);

          start();

        }

      }

    }

    ///@}

    /// \name Query Functions

    /// The results of the BFS algorithm can be obtained using these

    /// functions.\n

    /// Either \ref run(Node) "run()" or \ref start() should be called

    /// before using them.

    ///@{

    /// \brief Checks if a node is reached from the root(s).

    ///

    /// Returns \c true if \c v is reached from the root(s).

    ///

    /// \pre Either \ref run(Node) "run()" or \ref init()

    /// must be called before using this function.

    ///@}

  };

} //END OF NAMESPACE LEMON

#endif

    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.

    typedef _Traits Traits;

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

    typedef typename Traits::Digraph Digraph;

    ///The type of the flow values.

    typedef typename Traits::Value Value;

    /// The type of the lower bound capacity map.

    typedef typename Traits::LCapMap LCapMap;

    /// The type of the upper bound capacity map.

    typedef typename Traits::UCapMap UCapMap;

    /// \brief The type of the map that stores the lower bound for

    /// the supply of the nodes.

    typedef typename Traits::DeltaMap DeltaMap;

    ///The type of the flow map.

    typedef typename Traits::FlowMap FlowMap;

    ///The type of the elevator.

    typedef typename Traits::Elevator Elevator;

    ///The type of the tolerance.

    typedef typename Traits::Tolerance Tolerance;

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

    const Digraph &_g;

    int _node_num;

    const LCapMap *_lo;

    const UCapMap *_up;

    const DeltaMap *_delta;

    FlowMap *_flow;

    bool _local_flow;

    Elevator* _level;

    bool _local_level;

    typedef typename Digraph::template NodeMap<Value> ExcessMap;

    ExcessMap* _excess;

    Tolerance _tol;

    int _el;

  public:

    typedef Circulation Create;

    ///\name Named Template Parameters

    ///@{

    template <typename _FlowMap>

    struct SetFlowMapTraits : public Traits {

      typedef _FlowMap FlowMap;

      static FlowMap *createFlowMap(const Digraph&) {

        LEMON_ASSERT(false, "FlowMap is not initialized");

        return 0; // ignore warnings

      }

    };

    /// \brief \ref named-templ-param "Named parameter" for setting

    /// FlowMap type

    ///

    /// \ref named-templ-param "Named parameter" for setting FlowMap

    /// type.

    template <typename _FlowMap>

    struct SetFlowMap

      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,

                           SetFlowMapTraits<_FlowMap> > {

      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,

                          SetFlowMapTraits<_FlowMap> > Create;

    };

    template <typename _Elevator>

    struct SetElevatorTraits : public Traits {

      typedef _Elevator Elevator;

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

        LEMON_ASSERT(false, "Elevator is not initialized");

        return 0; // ignore warnings

      }

    };

    /// \brief \ref named-templ-param "Named parameter" for setting

    /// Elevator type

    ///

    /// \ref named-templ-param "Named parameter" for setting Elevator

    /// type. If this named parameter is used, then an external

    /// elevator object must be passed to the algorithm using the

    /// \ref elevator(Elevator&) "elevator()" function before calling

    /// \ref run() or \ref init().

    /// \sa SetStandardElevator

    template <typename _Elevator>

    struct SetElevator

      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,

                           SetElevatorTraits<_Elevator> > {

      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,

                          SetElevatorTraits<_Elevator> > Create;

    };

    template <typename _Elevator>

    struct SetStandardElevatorTraits : public Traits {

      typedef _Elevator Elevator;

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

        return new Elevator(digraph, max_level);

      }

    };

    /// \brief \ref named-templ-param "Named parameter" for setting

    /// Elevator type with automatic allocation

    ///

    /// \ref named-templ-param "Named parameter" for setting Elevator

    /// type with automatic allocation.

    /// The Elevator should have standard constructor interface to be

    /// able to automatically created by the algorithm (i.e. the

    /// digraph and the maximum level should be passed to it).

    /// However an external elevator object could also be passed to the

    /// algorithm with the \ref elevator(Elevator&) "elevator()" function

    /// before calling \ref run() or \ref init().

    /// \sa SetElevator

    template <typename _Elevator>

    struct SetStandardElevator

      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,

                       SetStandardElevatorTraits<_Elevator> > {

      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,

                      SetStandardElevatorTraits<_Elevator> > Create;

    };

    /// @}

  protected:

    Circulation() {}

  public:

    /// The constructor of the class.

    /// The constructor of the class.

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

    /// \param lo The lower bound capacity of the arcs.

    /// \param up The upper bound capacity of the arcs.

    /// \param delta The lower bound for the supply of the nodes.

    Circulation(const Digraph &g,const LCapMap &lo,

                const UCapMap &up,const DeltaMap &delta)

      : _g(g), _node_num(),

        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),

        _level(0), _local_level(false), _excess(0), _el() {}

    /// Destructor.

    ~Circulation() {

      destroyStructures();

    }

  private:

    void createStructures() {

      _node_num = _el = countNodes(_g);

      if (!_flow) {

        _flow = Traits::createFlowMap(_g);

        _local_flow = true;

      }

      if (!_level) {

        _level = Traits::createElevator(_g, _node_num);

        _local_level = true;

      }

      if (!_excess) {

        _excess = new ExcessMap(_g);

      }

    }

    void destroyStructures() {

      if (_local_flow) {

        delete _flow;

      }

      if (_local_level) {

        delete _level;

      }

      if (_excess) {

        delete _excess;

      }

    }

  public:

    /// Sets the lower bound capacity map.

    /// Sets the lower bound capacity map.

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

    Circulation& lowerCapMap(const LCapMap& map) {

      _lo = ↦

      return *this;

    }

    /// Sets the upper bound capacity map.

    /// Sets the upper bound capacity map.

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

    Circulation& upperCapMap(const LCapMap& map) {

      _up = ↦

      return *this;

    }

    /// Sets the lower bound map for the supply of the nodes.

    /// Sets the lower bound map for the supply of the nodes.

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

    Circulation& deltaMap(const DeltaMap& map) {

      _delta = ↦

      return *this;

    }

    /// \brief Sets the flow map.

    ///

    /// Sets the flow map.

    /// 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 map,

    /// of course.

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

    Circulation& flowMap(FlowMap& map) {

      if (_local_flow) {

        delete _flow;

        _local_flow = false;

      }

      _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.

      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

    /// If you need more control on the initial solution or the execution,

    /// first you have to call one of the \ref init() functions, then

    /// 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()

    {

      createStructures();

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

        _excess->set(n, (*_delta)[n]);

      }

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

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

        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);

        _excess->set(_g.source(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()

    {

      createStructures();

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

        _excess->set(n, (*_delta)[n]);

      }

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

        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {

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

          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);

          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);

        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {

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

          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);

          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);

        } else {

          Value fc = -(*_excess)[_g.target(e)];

          _flow->set(e, fc);

          _excess->set(_g.target(e), 0);

          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);

        }

      }

      _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);

    }

    ///Executes the algorithm

    ///This function executes the algorithm.

    ///

    ///\return \c true if a feasible circulation is found.

    ///

    ///\sa barrier()

    ///\sa barrierMap()

    bool start()

    {

      Node act;

      Node bact=INVALID;

      Node last_activated=INVALID;

      while((act=_level->highestActive())!=INVALID) {

        int actlevel=(*_level)[act];

        int mlevel=_node_num;

        Value exc=(*_excess)[act];

        for(OutArcIt e(_g,act);e!=INVALID; ++e) {

          Node v = _g.target(e);

          Value fc=(*_up)[e]-(*_flow)[e];

          if(!_tol.positive(fc)) continue;

          if((*_level)[v]<actlevel) {

            if(!_tol.less(fc, exc)) {

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

              _excess->set(v, (*_excess)[v] + exc);

              if(!_level->active(v) && _tol.positive((*_excess)[v]))

                _level->activate(v);

              _excess->set(act,0);

              _level->deactivate(act);

              goto next_l;

            }

            else {

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

              _excess->set(v, (*_excess)[v] + fc);

              if(!_level->active(v) && _tol.positive((*_excess)[v]))

                _level->activate(v);

              exc-=fc;

            }

          }

          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];

        }

        for(InArcIt e(_g,act);e!=INVALID; ++e) {

          Node v = _g.source(e);

          Value fc=(*_flow)[e]-(*_lo)[e];

          if(!_tol.positive(fc)) continue;

          if((*_level)[v]<actlevel) {

            if(!_tol.less(fc, exc)) {

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

              _excess->set(v, (*_excess)[v] + exc);

              if(!_level->active(v) && _tol.positive((*_excess)[v]))

                _level->activate(v);

              _excess->set(act,0);

              _level->deactivate(act);

              goto next_l;

            }

            else {

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

              _excess->set(v, (*_excess)[v] + fc);

              if(!_level->active(v) && _tol.positive((*_excess)[v]))

                _level->activate(v);

              exc-=fc;

            }

          }

          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];

        }

        _excess->set(act, exc);

        if(!_tol.positive(exc)) _level->deactivate(act);

        else if(mlevel==_node_num) {

          _level->liftHighestActiveToTop();

          _el = _node_num;

          return false;

        }

        else {

          _level->liftHighestActive(mlevel+1);

          if(_level->onLevel(actlevel)==0) {

            _el = actlevel;

            return false;

          }

        }

      next_l:

        ;

      }

      return true;

    }

    /// Runs the algorithm.

    /// This function runs the algorithm.

    ///

    /// \return \c true if a feasible circulation is found.

    ///

    /// \note Apart from the return value, c.run() is just a shortcut of

    /// the following code.

    /// \code

    ///   c.greedyInit();

    ///   c.start();

    /// \endcode

    bool run() {

      greedyInit();

      return start();

    }

    /// @}

    /// \name Query Functions

    /// The results of the circulation algorithm can be obtained using

    /// these functions.\n

    /// Either \ref run() or \ref start() should be called before

    /// using them.

    ///@{

    /// \brief Returns the flow on the given arc.

    ///

    /// Returns the flow on the given arc.

    ///

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

    /// using this function.

    Value flow(const Arc& arc) const {

      return (*_flow)[arc];

    }

    /// \brief Returns a const reference to the flow map.

    ///

    /// Returns a const reference to the arc map storing the found flow.

    ///

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

    /// using this function.

      return *_flow;

    }

    /**

       \brief Returns \c true if the given node is in a barrier.

       Barrier is a set \e B of nodes for which

       \f[ \sum_{a\in\delta_{out}(B)} upper(a) -

           \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]

       holds. The existence of a set with this property prooves that a

       feasible circualtion cannot exist.

       This function returns \c true if the given node is in the found

       barrier. If a feasible circulation is found, the function

       gives back \c false for every node.

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

       using this function.

       \sa barrierMap()

       \sa checkBarrier()

    */

    {

      return (*_level)[node] >= _el;

    }

    /// \brief Gives back a barrier.

    ///

    /// This function sets \c bar to the characteristic vector of the

    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"

    /// node map with \c bool (or convertible) value type.

    ///

    /// If a feasible circulation is found, the function gives back an

    /// empty set, so \c bar[v] will be \c false for all nodes \c v.

    ///

    /// \note This function calls \ref barrier() for each node,

    /// so it runs in \f$O(n)\f$ time.

    ///

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

    /// using this function.

    ///

    /// \sa barrier()

    /// \sa checkBarrier()

    template<class BarrierMap>

    {

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

        bar.set(n, (*_level)[n] >= _el);

    }

    /// @}

    /// \name Checker Functions

    /// The feasibility of the results can be checked using

    /// these functions.\n

    /// Either \ref run() or \ref start() should be called before

    /// using them.

    ///@{

    ///Check if the found flow is a feasible circulation

    ///Check if the found flow is a feasible circulation,

    ///

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

        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;

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

        {

          Value dif=-(*_delta)[n];

          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];

          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];

          if(_tol.negative(dif)) return false;

        }

      return true;

    }

    ///Check whether or not the last execution provides a barrier

    ///Check whether or not the last execution provides a barrier.

    ///\sa barrier()

    ///\sa barrierMap()

    {

      Value delta=0;

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

        if(barrier(n))

          delta-=(*_delta)[n];

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

        {

          Node s=_g.source(e);

          Node t=_g.target(e);

          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];

          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];

        }

      return _tol.negative(delta);

    }

    /// @}

  };

}

#endif

    /// If you need more control on the execution, first you have to call

    /// \ref init(), then you can add a source node with \ref addSource()

    /// 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);

          }

        }

    }

    /// \brief Processes the next arc.

    ///

    /// Processes the next arc.

    ///

    /// \return The processed arc.

    ///

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

    Arc processNextArc() {

      Arc e = _stack[_stack_head];

      Node m = _digraph->target(e);

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

        _visitor->discover(e);

        _visitor->reach(m);

        _reached->set(m, true);

        _digraph->firstOut(_stack[++_stack_head], m);

      } else {

        _visitor->examine(e);

        m = _digraph->source(e);

        _digraph->nextOut(_stack[_stack_head]);

      }

      while (_stack_head>=0 && _stack[_stack_head] == INVALID) {

        _visitor->leave(m);

        --_stack_head;

        if (_stack_head >= 0) {

          _visitor->backtrack(_stack[_stack_head]);

          m = _digraph->source(_stack[_stack_head]);

          _digraph->nextOut(_stack[_stack_head]);

        } else {

          _visitor->stop(m);

        }

      }

      return e;

    }

    /// \brief Next arc to be processed.

    ///

    /// Next arc to be processed.

    ///

    /// \return The next arc to be processed or INVALID if the stack is

    /// empty.

    Arc nextArc() const {

      return _stack_head >= 0 ? _stack[_stack_head] : INVALID;

    }

    /// \brief Returns \c false if there are nodes

    /// to be processed.

    ///

    /// Returns \c false if there are nodes

    /// to be processed in the queue (stack).

    bool emptyQueue() const { return _stack_head < 0; }

    /// \brief Returns the number of the nodes to be processed.

    ///

    /// Returns the number of the nodes to be processed in the queue (stack).

    int queueSize() const { return _stack_head + 1; }

    /// \brief Executes the algorithm.

    ///

    /// Executes the algorithm.

    ///

    /// This method runs the %DFS algorithm from the root node

    /// in order to compute the %DFS path to each node.

    ///

    /// The algorithm computes

    /// - the %DFS tree,

    /// - the distance of each node from the root in the %DFS tree.

    ///

    /// \pre init() must be called and a root node should be

    /// added with addSource() before using this function.

    ///

    /// \note <tt>d.start()</tt> is just a shortcut of the following code.

    /// \code

    ///   while ( !d.emptyQueue() ) {

    ///     d.processNextArc();

    ///   }

    /// \endcode

    void start() {

      while ( !emptyQueue() ) processNextArc();

    }

    /// \brief Executes the algorithm until the given target node is reached.

    ///

    /// Executes the algorithm until the given target node is reached.

    ///

    /// This method runs the %DFS algorithm from the root node

    /// in order to compute the DFS path to \c t.

    ///

    /// The algorithm computes

    /// - the %DFS path to \c t,

    /// - the distance of \c t from the root in the %DFS tree.

    ///

    /// \pre init() must be called and a root node should be added

    /// with addSource() before using this function.

    void start(Node t) {

      while ( !emptyQueue() && _digraph->target(_stack[_stack_head]) != t )

        processNextArc();

    }

    /// \brief Executes the algorithm until a condition is met.

    ///

    /// Executes the algorithm until a condition is met.

    ///

    /// This method runs the %DFS algorithm from the root node

    /// until an arc \c a with <tt>am[a]</tt> true is found.

    ///

    /// \param am A \c bool (or convertible) arc map. The algorithm

    /// will stop when it reaches an arc \c a with <tt>am[a]</tt> true.

    ///

    /// \return The reached arc \c a with <tt>am[a]</tt> true or

    /// \c INVALID if no such arc was found.

    ///

    /// \pre init() must be called and a root node should be added

    /// with addSource() before using this function.

    ///

    /// \warning Contrary to \ref Bfs and \ref Dijkstra, \c am is an arc map,

    /// not a node map.

    template <typename AM>

    Arc start(const AM &am) {

      while ( !emptyQueue() && !am[_stack[_stack_head]] )

        processNextArc();

      return emptyQueue() ? INVALID : _stack[_stack_head];

    }

    /// \brief Runs the algorithm from the given source node.

    ///

    /// This method runs the %DFS algorithm from node \c s.

    /// in order to compute the DFS path to each node.

    ///

    /// The algorithm computes

    /// - the %DFS tree,

    /// - the distance of each node from the root in the %DFS tree.

    ///

    /// \note <tt>d.run(s)</tt> is just a shortcut of the following code.

    ///\code

    ///   d.init();

    ///   d.addSource(s);

    ///   d.start();

    ///\endcode

    void run(Node s) {

      init();

      addSource(s);

      start();

    }

    /// \brief Finds the %DFS path between \c s and \c t.

    /// This method runs the %DFS algorithm from node \c s

    /// in order to compute the DFS path to node \c t

    /// (it stops searching when \c t is processed).

    ///

    /// \return \c true if \c t is reachable form \c s.

    ///

    /// \note Apart from the return value, <tt>d.run(s,t)</tt> is

    /// just a shortcut of the following code.

    ///\code

    ///   d.init();

    ///   d.addSource(s);

    ///   d.start(t);

    ///\endcode

    bool run(Node s,Node t) {

      init();

      addSource(s);

      start(t);

      return reached(t);

    }

    /// \brief Runs the algorithm to visit all nodes in the digraph.

    /// This method runs the %DFS algorithm in order to

    /// compute the %DFS path to each node.

    ///

    /// The algorithm computes

    /// - the %DFS tree (forest),

    /// - the distance of each node from the root(s) in the %DFS tree.

    ///

    /// \note <tt>d.run()</tt> is just a shortcut of the following code.

    ///\code

    ///   d.init();

    ///   for (NodeIt n(digraph); n != INVALID; ++n) {

    ///     if (!d.reached(n)) {

    ///       d.addSource(n);

    ///       d.start();

    ///     }

    ///   }

    ///\endcode

    void run() {

      init();

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

        if (!reached(it)) {

          addSource(it);

          start();

        }

      }

    }

    ///@}

    /// \name Query Functions