alpar@399: /* -*- mode: C++; indent-tabs-mode: nil; -*-
alpar@399:  *
alpar@399:  * This file is a part of LEMON, a generic C++ optimization library.
alpar@399:  *
alpar@877:  * Copyright (C) 2003-2010
alpar@399:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@399:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@399:  *
alpar@399:  * Permission to use, modify and distribute this software is granted
alpar@399:  * provided that this copyright notice appears in all copies. For
alpar@399:  * precise terms see the accompanying LICENSE file.
alpar@399:  *
alpar@399:  * This software is provided "AS IS" with no warranty of any kind,
alpar@399:  * express or implied, and with no claim as to its suitability for any
alpar@399:  * purpose.
alpar@399:  *
alpar@399:  */
alpar@399: 
alpar@399: #ifndef LEMON_CIRCULATION_H
alpar@399: #define LEMON_CIRCULATION_H
alpar@399: 
alpar@399: #include <lemon/tolerance.h>
alpar@399: #include <lemon/elevator.h>
kpeter@622: #include <limits>
alpar@399: 
alpar@399: ///\ingroup max_flow
alpar@399: ///\file
kpeter@402: ///\brief Push-relabel algorithm for finding a feasible circulation.
alpar@399: ///
alpar@399: namespace lemon {
alpar@399: 
alpar@399:   /// \brief Default traits class of Circulation class.
alpar@399:   ///
alpar@399:   /// Default traits class of Circulation class.
kpeter@610:   ///
kpeter@610:   /// \tparam GR Type of the digraph the algorithm runs on.
kpeter@610:   /// \tparam LM The type of the lower bound map.
kpeter@610:   /// \tparam UM The type of the upper bound (capacity) map.
kpeter@610:   /// \tparam SM The type of the supply map.
kpeter@492:   template <typename GR, typename LM,
kpeter@610:             typename UM, typename SM>
alpar@399:   struct CirculationDefaultTraits {
alpar@399: 
kpeter@402:     /// \brief The type of the digraph the algorithm runs on.
kpeter@492:     typedef GR Digraph;
alpar@399: 
kpeter@610:     /// \brief The type of the lower bound map.
alpar@399:     ///
kpeter@610:     /// The type of the map that stores the lower bounds on the arcs.
kpeter@610:     /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
kpeter@610:     typedef LM LowerMap;
alpar@399: 
kpeter@610:     /// \brief The type of the upper bound (capacity) map.
alpar@399:     ///
kpeter@610:     /// The type of the map that stores the upper bounds (capacities)
kpeter@610:     /// on the arcs.
kpeter@610:     /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
kpeter@610:     typedef UM UpperMap;
alpar@399: 
kpeter@610:     /// \brief The type of supply map.
alpar@399:     ///
alpar@877:     /// The type of the map that stores the signed supply values of the
alpar@877:     /// nodes.
kpeter@610:     /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
kpeter@610:     typedef SM SupplyMap;
alpar@399: 
kpeter@641:     /// \brief The type of the flow and supply values.
kpeter@641:     typedef typename SupplyMap::Value Value;
alpar@399: 
kpeter@402:     /// \brief The type of the map that stores the flow values.
alpar@399:     ///
kpeter@402:     /// The type of the map that stores the flow values.
kpeter@610:     /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
kpeter@610:     /// concept.
kpeter@713: #ifdef DOXYGEN
kpeter@713:     typedef GR::ArcMap<Value> FlowMap;
kpeter@713: #else
kpeter@641:     typedef typename Digraph::template ArcMap<Value> FlowMap;
kpeter@713: #endif
alpar@399: 
alpar@399:     /// \brief Instantiates a FlowMap.
alpar@399:     ///
alpar@399:     /// This function instantiates a \ref FlowMap.
kpeter@610:     /// \param digraph The digraph for which we would like to define
alpar@399:     /// the flow map.
alpar@399:     static FlowMap* createFlowMap(const Digraph& digraph) {
alpar@399:       return new FlowMap(digraph);
alpar@399:     }
alpar@399: 
kpeter@402:     /// \brief The elevator type used by the algorithm.
alpar@399:     ///
kpeter@402:     /// The elevator type used by the algorithm.
alpar@399:     ///
kpeter@713:     /// \sa Elevator, LinkedElevator
kpeter@713: #ifdef DOXYGEN
kpeter@713:     typedef lemon::Elevator<GR, GR::Node> Elevator;
kpeter@713: #else
alpar@399:     typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
kpeter@713: #endif
alpar@399: 
alpar@399:     /// \brief Instantiates an Elevator.
alpar@399:     ///
kpeter@402:     /// This function instantiates an \ref Elevator.
kpeter@610:     /// \param digraph The digraph for which we would like to define
alpar@399:     /// the elevator.
alpar@399:     /// \param max_level The maximum level of the elevator.
alpar@399:     static Elevator* createElevator(const Digraph& digraph, int max_level) {
alpar@399:       return new Elevator(digraph, max_level);
alpar@399:     }
alpar@399: 
alpar@399:     /// \brief The tolerance used by the algorithm
alpar@399:     ///
alpar@399:     /// The tolerance used by the algorithm to handle inexact computation.
kpeter@641:     typedef lemon::Tolerance<Value> Tolerance;
alpar@399: 
alpar@399:   };
alpar@399: 
kpeter@402:   /**
kpeter@402:      \brief Push-relabel algorithm for the network circulation problem.
alpar@399: 
alpar@399:      \ingroup max_flow
kpeter@610:      This class implements a push-relabel algorithm for the \e network
kpeter@610:      \e circulation problem.
kpeter@402:      It is to find a feasible circulation when lower and upper bounds
kpeter@610:      are given for the flow values on the arcs and lower bounds are
kpeter@610:      given for the difference between the outgoing and incoming flow
kpeter@610:      at the nodes.
kpeter@402: 
alpar@399:      The exact formulation of this problem is the following.
kpeter@622:      Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$
kpeter@622:      \f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and
kpeter@622:      upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$
kpeter@610:      holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
kpeter@610:      denotes the signed supply values of the nodes.
kpeter@610:      If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
kpeter@610:      supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
kpeter@610:      \f$-sup(u)\f$ demand.
kpeter@622:      A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$
kpeter@610:      solution of the following problem.
kpeter@402: 
kpeter@610:      \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
kpeter@610:      \geq sup(u) \quad \forall u\in V, \f]
kpeter@610:      \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]
alpar@877: 
kpeter@610:      The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
kpeter@610:      zero or negative in order to have a feasible solution (since the sum
kpeter@610:      of the expressions on the left-hand side of the inequalities is zero).
kpeter@610:      It means that the total demand must be greater or equal to the total
kpeter@610:      supply and all the supplies have to be carried out from the supply nodes,
kpeter@610:      but there could be demands that are not satisfied.
kpeter@610:      If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
kpeter@610:      constraints have to be satisfied with equality, i.e. all demands
kpeter@610:      have to be satisfied and all supplies have to be used.
alpar@877: 
kpeter@610:      If you need the opposite inequalities in the supply/demand constraints
kpeter@610:      (i.e. the total demand is less than the total supply and all the demands
kpeter@610:      have to be satisfied while there could be supplies that are not used),
kpeter@610:      then you could easily transform the problem to the above form by reversing
kpeter@610:      the direction of the arcs and taking the negative of the supply values
kpeter@610:      (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
kpeter@610: 
kpeter@622:      This algorithm either calculates a feasible circulation, or provides
kpeter@622:      a \ref barrier() "barrier", which prooves that a feasible soultion
kpeter@622:      cannot exist.
kpeter@622: 
kpeter@610:      Note that this algorithm also provides a feasible solution for the
kpeter@610:      \ref min_cost_flow "minimum cost flow problem".
kpeter@402: 
kpeter@492:      \tparam GR The type of the digraph the algorithm runs on.
kpeter@610:      \tparam LM The type of the lower bound map. The default
kpeter@492:      map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
kpeter@610:      \tparam UM The type of the upper bound (capacity) map.
kpeter@610:      The default map type is \c LM.
kpeter@610:      \tparam SM The type of the supply map. The default map type is
kpeter@492:      \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
kpeter@825:      \tparam TR The traits class that defines various types used by the
kpeter@825:      algorithm. By default, it is \ref CirculationDefaultTraits
kpeter@825:      "CirculationDefaultTraits<GR, LM, UM, SM>".
kpeter@825:      In most cases, this parameter should not be set directly,
kpeter@825:      consider to use the named template parameters instead.
alpar@399:   */
kpeter@402: #ifdef DOXYGEN
kpeter@492: template< typename GR,
kpeter@492:           typename LM,
kpeter@492:           typename UM,
kpeter@610:           typename SM,
kpeter@492:           typename TR >
kpeter@402: #else
kpeter@492: template< typename GR,
kpeter@492:           typename LM = typename GR::template ArcMap<int>,
kpeter@492:           typename UM = LM,
kpeter@610:           typename SM = typename GR::template NodeMap<typename UM::Value>,
kpeter@610:           typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
kpeter@402: #endif
alpar@399:   class Circulation {
kpeter@402:   public:
alpar@399: 
kpeter@402:     ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
kpeter@492:     typedef TR Traits;
kpeter@402:     ///The type of the digraph the algorithm runs on.
alpar@399:     typedef typename Traits::Digraph Digraph;
kpeter@641:     ///The type of the flow and supply values.
kpeter@641:     typedef typename Traits::Value Value;
alpar@399: 
kpeter@610:     ///The type of the lower bound map.
kpeter@610:     typedef typename Traits::LowerMap LowerMap;
kpeter@610:     ///The type of the upper bound (capacity) map.
kpeter@610:     typedef typename Traits::UpperMap UpperMap;
kpeter@610:     ///The type of the supply map.
kpeter@610:     typedef typename Traits::SupplyMap SupplyMap;
kpeter@402:     ///The type of the flow map.
alpar@399:     typedef typename Traits::FlowMap FlowMap;
kpeter@402: 
kpeter@402:     ///The type of the elevator.
alpar@399:     typedef typename Traits::Elevator Elevator;
kpeter@402:     ///The type of the tolerance.
alpar@399:     typedef typename Traits::Tolerance Tolerance;
alpar@399: 
kpeter@402:   private:
kpeter@402: 
kpeter@402:     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
alpar@399: 
alpar@399:     const Digraph &_g;
alpar@399:     int _node_num;
alpar@399: 
kpeter@610:     const LowerMap *_lo;
kpeter@610:     const UpperMap *_up;
kpeter@610:     const SupplyMap *_supply;
alpar@399: 
alpar@399:     FlowMap *_flow;
alpar@399:     bool _local_flow;
alpar@399: 
alpar@399:     Elevator* _level;
alpar@399:     bool _local_level;
alpar@399: 
kpeter@641:     typedef typename Digraph::template NodeMap<Value> ExcessMap;
alpar@399:     ExcessMap* _excess;
alpar@399: 
alpar@399:     Tolerance _tol;
alpar@399:     int _el;
alpar@399: 
alpar@399:   public:
alpar@399: 
alpar@399:     typedef Circulation Create;
alpar@399: 
kpeter@402:     ///\name Named Template Parameters
alpar@399: 
alpar@399:     ///@{
alpar@399: 
kpeter@559:     template <typename T>
alpar@401:     struct SetFlowMapTraits : public Traits {
kpeter@559:       typedef T FlowMap;
alpar@399:       static FlowMap *createFlowMap(const Digraph&) {
alpar@399:         LEMON_ASSERT(false, "FlowMap is not initialized");
alpar@399:         return 0; // ignore warnings
alpar@399:       }
alpar@399:     };
alpar@399: 
alpar@399:     /// \brief \ref named-templ-param "Named parameter" for setting
alpar@399:     /// FlowMap type
alpar@399:     ///
alpar@399:     /// \ref named-templ-param "Named parameter" for setting FlowMap
kpeter@402:     /// type.
kpeter@559:     template <typename T>
alpar@401:     struct SetFlowMap
kpeter@610:       : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
kpeter@559:                            SetFlowMapTraits<T> > {
kpeter@610:       typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
kpeter@559:                           SetFlowMapTraits<T> > Create;
alpar@399:     };
alpar@399: 
kpeter@559:     template <typename T>
alpar@401:     struct SetElevatorTraits : public Traits {
kpeter@559:       typedef T Elevator;
alpar@399:       static Elevator *createElevator(const Digraph&, int) {
alpar@399:         LEMON_ASSERT(false, "Elevator is not initialized");
alpar@399:         return 0; // ignore warnings
alpar@399:       }
alpar@399:     };
alpar@399: 
alpar@399:     /// \brief \ref named-templ-param "Named parameter" for setting
alpar@399:     /// Elevator type
alpar@399:     ///
alpar@399:     /// \ref named-templ-param "Named parameter" for setting Elevator
kpeter@402:     /// type. If this named parameter is used, then an external
kpeter@402:     /// elevator object must be passed to the algorithm using the
kpeter@402:     /// \ref elevator(Elevator&) "elevator()" function before calling
kpeter@402:     /// \ref run() or \ref init().
kpeter@402:     /// \sa SetStandardElevator
kpeter@559:     template <typename T>
alpar@401:     struct SetElevator
kpeter@610:       : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
kpeter@559:                            SetElevatorTraits<T> > {
kpeter@610:       typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
kpeter@559:                           SetElevatorTraits<T> > Create;
alpar@399:     };
alpar@399: 
kpeter@559:     template <typename T>
alpar@401:     struct SetStandardElevatorTraits : public Traits {
kpeter@559:       typedef T Elevator;
alpar@399:       static Elevator *createElevator(const Digraph& digraph, int max_level) {
alpar@399:         return new Elevator(digraph, max_level);
alpar@399:       }
alpar@399:     };
alpar@399: 
alpar@399:     /// \brief \ref named-templ-param "Named parameter" for setting
kpeter@402:     /// Elevator type with automatic allocation
alpar@399:     ///
alpar@399:     /// \ref named-templ-param "Named parameter" for setting Elevator
kpeter@402:     /// type with automatic allocation.
kpeter@402:     /// The Elevator should have standard constructor interface to be
kpeter@402:     /// able to automatically created by the algorithm (i.e. the
kpeter@402:     /// digraph and the maximum level should be passed to it).
kpeter@786:     /// However, an external elevator object could also be passed to the
kpeter@402:     /// algorithm with the \ref elevator(Elevator&) "elevator()" function
kpeter@402:     /// before calling \ref run() or \ref init().
kpeter@402:     /// \sa SetElevator
kpeter@559:     template <typename T>
alpar@401:     struct SetStandardElevator
kpeter@610:       : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
kpeter@559:                        SetStandardElevatorTraits<T> > {
kpeter@610:       typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
kpeter@559:                       SetStandardElevatorTraits<T> > Create;
alpar@399:     };
alpar@399: 
alpar@399:     /// @}
alpar@399: 
alpar@399:   protected:
alpar@399: 
alpar@399:     Circulation() {}
alpar@399: 
alpar@399:   public:
alpar@399: 
kpeter@610:     /// Constructor.
alpar@399: 
alpar@399:     /// The constructor of the class.
kpeter@610:     ///
kpeter@610:     /// \param graph The digraph the algorithm runs on.
kpeter@610:     /// \param lower The lower bounds for the flow values on the arcs.
alpar@877:     /// \param upper The upper bounds (capacities) for the flow values
kpeter@610:     /// on the arcs.
kpeter@610:     /// \param supply The signed supply values of the nodes.
kpeter@610:     Circulation(const Digraph &graph, const LowerMap &lower,
kpeter@610:                 const UpperMap &upper, const SupplyMap &supply)
kpeter@610:       : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
kpeter@610:         _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
kpeter@610:         _excess(NULL) {}
alpar@399: 
kpeter@402:     /// Destructor.
alpar@399:     ~Circulation() {
alpar@399:       destroyStructures();
alpar@399:     }
alpar@399: 
kpeter@402: 
alpar@399:   private:
alpar@399: 
kpeter@622:     bool checkBoundMaps() {
kpeter@622:       for (ArcIt e(_g);e!=INVALID;++e) {
kpeter@622:         if (_tol.less((*_up)[e], (*_lo)[e])) return false;
kpeter@622:       }
kpeter@622:       return true;
kpeter@622:     }
kpeter@622: 
alpar@399:     void createStructures() {
alpar@399:       _node_num = _el = countNodes(_g);
alpar@399: 
alpar@399:       if (!_flow) {
alpar@399:         _flow = Traits::createFlowMap(_g);
alpar@399:         _local_flow = true;
alpar@399:       }
alpar@399:       if (!_level) {
alpar@399:         _level = Traits::createElevator(_g, _node_num);
alpar@399:         _local_level = true;
alpar@399:       }
alpar@399:       if (!_excess) {
alpar@399:         _excess = new ExcessMap(_g);
alpar@399:       }
alpar@399:     }
alpar@399: 
alpar@399:     void destroyStructures() {
alpar@399:       if (_local_flow) {
alpar@399:         delete _flow;
alpar@399:       }
alpar@399:       if (_local_level) {
alpar@399:         delete _level;
alpar@399:       }
alpar@399:       if (_excess) {
alpar@399:         delete _excess;
alpar@399:       }
alpar@399:     }
alpar@399: 
alpar@399:   public:
alpar@399: 
kpeter@610:     /// Sets the lower bound map.
alpar@399: 
kpeter@610:     /// Sets the lower bound map.
kpeter@402:     /// \return <tt>(*this)</tt>
kpeter@610:     Circulation& lowerMap(const LowerMap& map) {
alpar@399:       _lo = &map;
alpar@399:       return *this;
alpar@399:     }
alpar@399: 
kpeter@610:     /// Sets the upper bound (capacity) map.
alpar@399: 
kpeter@610:     /// Sets the upper bound (capacity) map.
kpeter@402:     /// \return <tt>(*this)</tt>
kpeter@622:     Circulation& upperMap(const UpperMap& map) {
alpar@399:       _up = &map;
alpar@399:       return *this;
alpar@399:     }
alpar@399: 
kpeter@610:     /// Sets the supply map.
alpar@399: 
kpeter@610:     /// Sets the supply map.
kpeter@402:     /// \return <tt>(*this)</tt>
kpeter@610:     Circulation& supplyMap(const SupplyMap& map) {
kpeter@610:       _supply = &map;
alpar@399:       return *this;
alpar@399:     }
alpar@399: 
kpeter@402:     /// \brief Sets the flow map.
kpeter@402:     ///
alpar@399:     /// Sets the flow map.
kpeter@402:     /// If you don't use this function before calling \ref run() or
kpeter@402:     /// \ref init(), an instance will be allocated automatically.
kpeter@402:     /// The destructor deallocates this automatically allocated map,
kpeter@402:     /// of course.
kpeter@402:     /// \return <tt>(*this)</tt>
alpar@399:     Circulation& flowMap(FlowMap& map) {
alpar@399:       if (_local_flow) {
alpar@399:         delete _flow;
alpar@399:         _local_flow = false;
alpar@399:       }
alpar@399:       _flow = &map;
alpar@399:       return *this;
alpar@399:     }
alpar@399: 
kpeter@402:     /// \brief Sets the elevator used by algorithm.
alpar@399:     ///
kpeter@402:     /// Sets the elevator used by algorithm.
kpeter@402:     /// If you don't use this function before calling \ref run() or
kpeter@402:     /// \ref init(), an instance will be allocated automatically.
kpeter@402:     /// The destructor deallocates this automatically allocated elevator,
kpeter@402:     /// of course.
kpeter@402:     /// \return <tt>(*this)</tt>
alpar@399:     Circulation& elevator(Elevator& elevator) {
alpar@399:       if (_local_level) {
alpar@399:         delete _level;
alpar@399:         _local_level = false;
alpar@399:       }
alpar@399:       _level = &elevator;
alpar@399:       return *this;
alpar@399:     }
alpar@399: 
kpeter@402:     /// \brief Returns a const reference to the elevator.
alpar@399:     ///
kpeter@402:     /// Returns a const reference to the elevator.
kpeter@402:     ///
kpeter@402:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@402:     /// using this function.
kpeter@420:     const Elevator& elevator() const {
alpar@399:       return *_level;
alpar@399:     }
alpar@399: 
kpeter@689:     /// \brief Sets the tolerance used by the algorithm.
kpeter@402:     ///
kpeter@689:     /// Sets the tolerance object used by the algorithm.
kpeter@689:     /// \return <tt>(*this)</tt>
kpeter@688:     Circulation& tolerance(const Tolerance& tolerance) {
alpar@399:       _tol = tolerance;
alpar@399:       return *this;
alpar@399:     }
alpar@399: 
kpeter@402:     /// \brief Returns a const reference to the tolerance.
alpar@399:     ///
kpeter@689:     /// Returns a const reference to the tolerance object used by
kpeter@689:     /// the algorithm.
alpar@399:     const Tolerance& tolerance() const {
kpeter@688:       return _tol;
alpar@399:     }
alpar@399: 
kpeter@402:     /// \name Execution Control
kpeter@402:     /// The simplest way to execute the algorithm is to call \ref run().\n
kpeter@713:     /// If you need better control on the initial solution or the execution,
kpeter@713:     /// you have to call one of the \ref init() functions first, then
kpeter@402:     /// the \ref start() function.
alpar@399: 
alpar@399:     ///@{
alpar@399: 
alpar@399:     /// Initializes the internal data structures.
alpar@399: 
kpeter@402:     /// Initializes the internal data structures and sets all flow values
kpeter@402:     /// to the lower bound.
alpar@399:     void init()
alpar@399:     {
kpeter@622:       LEMON_DEBUG(checkBoundMaps(),
kpeter@622:         "Upper bounds must be greater or equal to the lower bounds");
kpeter@622: 
alpar@399:       createStructures();
alpar@399: 
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n) {
alpar@611:         (*_excess)[n] = (*_supply)[n];
alpar@399:       }
alpar@399: 
alpar@399:       for (ArcIt e(_g);e!=INVALID;++e) {
alpar@399:         _flow->set(e, (*_lo)[e]);
kpeter@581:         (*_excess)[_g.target(e)] += (*_flow)[e];
kpeter@581:         (*_excess)[_g.source(e)] -= (*_flow)[e];
alpar@399:       }
alpar@399: 
alpar@399:       // global relabeling tested, but in general case it provides
alpar@399:       // worse performance for random digraphs
alpar@399:       _level->initStart();
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@399:         _level->initAddItem(n);
alpar@399:       _level->initFinish();
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@399:         if(_tol.positive((*_excess)[n]))
alpar@399:           _level->activate(n);
alpar@399:     }
alpar@399: 
kpeter@402:     /// Initializes the internal data structures using a greedy approach.
alpar@399: 
kpeter@402:     /// Initializes the internal data structures using a greedy approach
kpeter@402:     /// to construct the initial solution.
alpar@399:     void greedyInit()
alpar@399:     {
kpeter@622:       LEMON_DEBUG(checkBoundMaps(),
kpeter@622:         "Upper bounds must be greater or equal to the lower bounds");
kpeter@622: 
alpar@399:       createStructures();
alpar@399: 
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n) {
alpar@611:         (*_excess)[n] = (*_supply)[n];
alpar@399:       }
alpar@399: 
alpar@399:       for (ArcIt e(_g);e!=INVALID;++e) {
kpeter@622:         if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
alpar@399:           _flow->set(e, (*_up)[e]);
kpeter@581:           (*_excess)[_g.target(e)] += (*_up)[e];
kpeter@581:           (*_excess)[_g.source(e)] -= (*_up)[e];
kpeter@622:         } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
alpar@399:           _flow->set(e, (*_lo)[e]);
kpeter@581:           (*_excess)[_g.target(e)] += (*_lo)[e];
kpeter@581:           (*_excess)[_g.source(e)] -= (*_lo)[e];
alpar@399:         } else {
kpeter@641:           Value fc = -(*_excess)[_g.target(e)];
alpar@399:           _flow->set(e, fc);
kpeter@581:           (*_excess)[_g.target(e)] = 0;
kpeter@581:           (*_excess)[_g.source(e)] -= fc;
alpar@399:         }
alpar@399:       }
alpar@399: 
alpar@399:       _level->initStart();
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@399:         _level->initAddItem(n);
alpar@399:       _level->initFinish();
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@399:         if(_tol.positive((*_excess)[n]))
alpar@399:           _level->activate(n);
alpar@399:     }
alpar@399: 
kpeter@402:     ///Executes the algorithm
alpar@399: 
kpeter@402:     ///This function executes the algorithm.
kpeter@402:     ///
kpeter@402:     ///\return \c true if a feasible circulation is found.
alpar@399:     ///
alpar@399:     ///\sa barrier()
kpeter@402:     ///\sa barrierMap()
alpar@399:     bool start()
alpar@399:     {
alpar@399: 
alpar@399:       Node act;
alpar@399:       Node bact=INVALID;
alpar@399:       Node last_activated=INVALID;
alpar@399:       while((act=_level->highestActive())!=INVALID) {
alpar@399:         int actlevel=(*_level)[act];
alpar@399:         int mlevel=_node_num;
kpeter@641:         Value exc=(*_excess)[act];
alpar@399: 
alpar@399:         for(OutArcIt e(_g,act);e!=INVALID; ++e) {
alpar@399:           Node v = _g.target(e);
kpeter@641:           Value fc=(*_up)[e]-(*_flow)[e];
alpar@399:           if(!_tol.positive(fc)) continue;
alpar@399:           if((*_level)[v]<actlevel) {
alpar@399:             if(!_tol.less(fc, exc)) {
alpar@399:               _flow->set(e, (*_flow)[e] + exc);
kpeter@581:               (*_excess)[v] += exc;
alpar@399:               if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@399:                 _level->activate(v);
kpeter@581:               (*_excess)[act] = 0;
alpar@399:               _level->deactivate(act);
alpar@399:               goto next_l;
alpar@399:             }
alpar@399:             else {
alpar@399:               _flow->set(e, (*_up)[e]);
kpeter@581:               (*_excess)[v] += fc;
alpar@399:               if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@399:                 _level->activate(v);
alpar@399:               exc-=fc;
alpar@399:             }
alpar@399:           }
alpar@399:           else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
alpar@399:         }
alpar@399:         for(InArcIt e(_g,act);e!=INVALID; ++e) {
alpar@399:           Node v = _g.source(e);
kpeter@641:           Value fc=(*_flow)[e]-(*_lo)[e];
alpar@399:           if(!_tol.positive(fc)) continue;
alpar@399:           if((*_level)[v]<actlevel) {
alpar@399:             if(!_tol.less(fc, exc)) {
alpar@399:               _flow->set(e, (*_flow)[e] - exc);
kpeter@581:               (*_excess)[v] += exc;
alpar@399:               if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@399:                 _level->activate(v);
kpeter@581:               (*_excess)[act] = 0;
alpar@399:               _level->deactivate(act);
alpar@399:               goto next_l;
alpar@399:             }
alpar@399:             else {
alpar@399:               _flow->set(e, (*_lo)[e]);
kpeter@581:               (*_excess)[v] += fc;
alpar@399:               if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@399:                 _level->activate(v);
alpar@399:               exc-=fc;
alpar@399:             }
alpar@399:           }
alpar@399:           else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
alpar@399:         }
alpar@399: 
kpeter@581:         (*_excess)[act] = exc;
alpar@399:         if(!_tol.positive(exc)) _level->deactivate(act);
alpar@399:         else if(mlevel==_node_num) {
alpar@399:           _level->liftHighestActiveToTop();
alpar@399:           _el = _node_num;
alpar@399:           return false;
alpar@399:         }
alpar@399:         else {
alpar@399:           _level->liftHighestActive(mlevel+1);
alpar@399:           if(_level->onLevel(actlevel)==0) {
alpar@399:             _el = actlevel;
alpar@399:             return false;
alpar@399:           }
alpar@399:         }
alpar@399:       next_l:
alpar@399:         ;
alpar@399:       }
alpar@399:       return true;
alpar@399:     }
alpar@399: 
kpeter@402:     /// Runs the algorithm.
alpar@399: 
kpeter@402:     /// This function runs the algorithm.
kpeter@402:     ///
kpeter@402:     /// \return \c true if a feasible circulation is found.
kpeter@402:     ///
kpeter@402:     /// \note Apart from the return value, c.run() is just a shortcut of
kpeter@402:     /// the following code.
alpar@399:     /// \code
kpeter@402:     ///   c.greedyInit();
kpeter@402:     ///   c.start();
alpar@399:     /// \endcode
alpar@399:     bool run() {
alpar@399:       greedyInit();
alpar@399:       return start();
alpar@399:     }
alpar@399: 
alpar@399:     /// @}
alpar@399: 
alpar@399:     /// \name Query Functions
kpeter@402:     /// The results of the circulation algorithm can be obtained using
kpeter@402:     /// these functions.\n
kpeter@402:     /// Either \ref run() or \ref start() should be called before
kpeter@402:     /// using them.
alpar@399: 
alpar@399:     ///@{
alpar@399: 
kpeter@641:     /// \brief Returns the flow value on the given arc.
kpeter@402:     ///
kpeter@641:     /// Returns the flow value on the given arc.
kpeter@402:     ///
kpeter@402:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@402:     /// using this function.
kpeter@641:     Value flow(const Arc& arc) const {
kpeter@402:       return (*_flow)[arc];
kpeter@402:     }
kpeter@402: 
kpeter@402:     /// \brief Returns a const reference to the flow map.
kpeter@402:     ///
kpeter@402:     /// Returns a const reference to the arc map storing the found flow.
kpeter@402:     ///
kpeter@402:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@402:     /// using this function.
kpeter@420:     const FlowMap& flowMap() const {
kpeter@402:       return *_flow;
kpeter@402:     }
kpeter@402: 
alpar@399:     /**
kpeter@402:        \brief Returns \c true if the given node is in a barrier.
kpeter@402: 
alpar@399:        Barrier is a set \e B of nodes for which
kpeter@402: 
kpeter@610:        \f[ \sum_{uv\in A: u\in B} upper(uv) -
kpeter@610:            \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
kpeter@402: 
kpeter@402:        holds. The existence of a set with this property prooves that a
kpeter@402:        feasible circualtion cannot exist.
kpeter@402: 
kpeter@402:        This function returns \c true if the given node is in the found
kpeter@402:        barrier. If a feasible circulation is found, the function
kpeter@402:        gives back \c false for every node.
kpeter@402: 
kpeter@402:        \pre Either \ref run() or \ref init() must be called before
kpeter@402:        using this function.
kpeter@402: 
kpeter@402:        \sa barrierMap()
alpar@399:        \sa checkBarrier()
alpar@399:     */
kpeter@420:     bool barrier(const Node& node) const
kpeter@402:     {
kpeter@402:       return (*_level)[node] >= _el;
kpeter@402:     }
kpeter@402: 
kpeter@402:     /// \brief Gives back a barrier.
kpeter@402:     ///
kpeter@402:     /// This function sets \c bar to the characteristic vector of the
kpeter@402:     /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
kpeter@402:     /// node map with \c bool (or convertible) value type.
kpeter@402:     ///
kpeter@402:     /// If a feasible circulation is found, the function gives back an
kpeter@402:     /// empty set, so \c bar[v] will be \c false for all nodes \c v.
kpeter@402:     ///
kpeter@402:     /// \note This function calls \ref barrier() for each node,
kpeter@559:     /// so it runs in O(n) time.
kpeter@402:     ///
kpeter@402:     /// \pre Either \ref run() or \ref init() must be called before
kpeter@402:     /// using this function.
kpeter@402:     ///
kpeter@402:     /// \sa barrier()
kpeter@402:     /// \sa checkBarrier()
kpeter@402:     template<class BarrierMap>
kpeter@420:     void barrierMap(BarrierMap &bar) const
alpar@399:     {
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@399:         bar.set(n, (*_level)[n] >= _el);
alpar@399:     }
alpar@399: 
alpar@399:     /// @}
alpar@399: 
alpar@399:     /// \name Checker Functions
kpeter@402:     /// The feasibility of the results can be checked using
kpeter@402:     /// these functions.\n
kpeter@402:     /// Either \ref run() or \ref start() should be called before
kpeter@402:     /// using them.
alpar@399: 
alpar@399:     ///@{
alpar@399: 
kpeter@402:     ///Check if the found flow is a feasible circulation
kpeter@402: 
kpeter@402:     ///Check if the found flow is a feasible circulation,
kpeter@402:     ///
kpeter@420:     bool checkFlow() const {
alpar@399:       for(ArcIt e(_g);e!=INVALID;++e)
alpar@399:         if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@399:         {
kpeter@641:           Value dif=-(*_supply)[n];
alpar@399:           for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
alpar@399:           for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
alpar@399:           if(_tol.negative(dif)) return false;
alpar@399:         }
alpar@399:       return true;
alpar@399:     }
alpar@399: 
alpar@399:     ///Check whether or not the last execution provides a barrier
alpar@399: 
kpeter@402:     ///Check whether or not the last execution provides a barrier.
alpar@399:     ///\sa barrier()
kpeter@402:     ///\sa barrierMap()
kpeter@420:     bool checkBarrier() const
alpar@399:     {
kpeter@641:       Value delta=0;
kpeter@641:       Value inf_cap = std::numeric_limits<Value>::has_infinity ?
kpeter@641:         std::numeric_limits<Value>::infinity() :
kpeter@641:         std::numeric_limits<Value>::max();
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@399:         if(barrier(n))
kpeter@610:           delta-=(*_supply)[n];
alpar@399:       for(ArcIt e(_g);e!=INVALID;++e)
alpar@399:         {
alpar@399:           Node s=_g.source(e);
alpar@399:           Node t=_g.target(e);
kpeter@622:           if(barrier(s)&&!barrier(t)) {
kpeter@622:             if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
kpeter@622:             delta+=(*_up)[e];
kpeter@622:           }
alpar@399:           else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
alpar@399:         }
alpar@399:       return _tol.negative(delta);
alpar@399:     }
alpar@399: 
alpar@399:     /// @}
alpar@399: 
alpar@399:   };
alpar@399: 
alpar@399: }
alpar@399: 
alpar@399: #endif