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