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@440:  * Copyright (C) 2003-2009
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>
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@492:   /// \tparam GR Digraph type.
kpeter@492:   /// \tparam LM Lower bound capacity map type.
kpeter@492:   /// \tparam UM Upper bound capacity map type.
kpeter@492:   /// \tparam DM Delta map type.
kpeter@492:   template <typename GR, typename LM,
kpeter@492:             typename UM, typename DM>
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: 
alpar@399:     /// \brief The type of the map that stores the circulation lower
alpar@399:     /// bound.
alpar@399:     ///
alpar@399:     /// The type of the map that stores the circulation lower bound.
alpar@399:     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
kpeter@492:     typedef LM LCapMap;
alpar@399: 
alpar@399:     /// \brief The type of the map that stores the circulation upper
alpar@399:     /// bound.
alpar@399:     ///
alpar@399:     /// The type of the map that stores the circulation upper bound.
alpar@399:     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
kpeter@492:     typedef UM UCapMap;
alpar@399: 
kpeter@402:     /// \brief The type of the map that stores the lower bound for
kpeter@402:     /// the supply of the nodes.
alpar@399:     ///
kpeter@402:     /// The type of the map that stores the lower bound for the supply
kpeter@402:     /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
alpar@399:     /// concept.
kpeter@492:     typedef DM DeltaMap;
alpar@399: 
kpeter@402:     /// \brief The type of the flow values.
alpar@399:     typedef typename DeltaMap::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.
alpar@399:     /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
alpar@399:     typedef typename Digraph::template ArcMap<Value> FlowMap;
alpar@399: 
alpar@399:     /// \brief Instantiates a FlowMap.
alpar@399:     ///
alpar@399:     /// This function instantiates a \ref FlowMap.
alpar@399:     /// \param digraph The digraph, to 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:     ///
alpar@399:     /// \sa Elevator
alpar@399:     /// \sa LinkedElevator
alpar@399:     typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
alpar@399: 
alpar@399:     /// \brief Instantiates an Elevator.
alpar@399:     ///
kpeter@402:     /// This function instantiates an \ref Elevator.
alpar@399:     /// \param digraph The digraph, to 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.
alpar@399:     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@402:      This class implements a push-relabel algorithm for the network
kpeter@402:      circulation problem.
kpeter@402:      It is to find a feasible circulation when lower and upper bounds
kpeter@402:      are given for the flow values on the arcs and lower bounds
kpeter@402:      are given for the supply values of the nodes.
kpeter@402: 
alpar@399:      The exact formulation of this problem is the following.
kpeter@402:      Let \f$G=(V,A)\f$ be a digraph,
kpeter@402:      \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
kpeter@402:      \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
kpeter@402:      \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
kpeter@402:      \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
kpeter@402:      \geq delta(v) \quad \forall v\in V, \f]
kpeter@402:      \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
kpeter@402:      \note \f$delta(v)\f$ specifies a lower bound for the supply of node
kpeter@402:      \f$v\f$. It can be either positive or negative, however note that
kpeter@402:      \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
kpeter@402:      have a feasible solution.
kpeter@402: 
kpeter@402:      \note A special case of this problem is when
kpeter@402:      \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
kpeter@402:      will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
kpeter@402:      Thus a feasible solution for the
kpeter@402:      \ref min_cost_flow "minimum cost flow" problem can be calculated
kpeter@402:      in this way.
kpeter@402: 
kpeter@492:      \tparam GR The type of the digraph the algorithm runs on.
kpeter@492:      \tparam LM The type of the lower bound capacity map. The default
kpeter@492:      map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
kpeter@492:      \tparam UM The type of the upper bound capacity map. The default
kpeter@492:      map type is \c LM.
kpeter@492:      \tparam DM The type of the map that stores the lower bound
kpeter@402:      for the supply of the nodes. The default map type is
kpeter@492:      \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
alpar@399:   */
kpeter@402: #ifdef DOXYGEN
kpeter@492: template< typename GR,
kpeter@492:           typename LM,
kpeter@492:           typename UM,
kpeter@492:           typename DM,
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@492:           typename DM = typename GR::template NodeMap<typename UM::Value>,
kpeter@492:           typename TR = CirculationDefaultTraits<GR, LM, UM, DM> >
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@402:     ///The type of the flow values.
alpar@399:     typedef typename Traits::Value Value;
alpar@399: 
kpeter@402:     /// The type of the lower bound capacity map.
alpar@399:     typedef typename Traits::LCapMap LCapMap;
kpeter@402:     /// The type of the upper bound capacity map.
alpar@399:     typedef typename Traits::UCapMap UCapMap;
kpeter@402:     /// \brief The type of the map that stores the lower bound for
kpeter@402:     /// the supply of the nodes.
alpar@399:     typedef typename Traits::DeltaMap DeltaMap;
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: 
alpar@399:     const LCapMap *_lo;
alpar@399:     const UCapMap *_up;
alpar@399:     const DeltaMap *_delta;
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@402:     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
alpar@399:       : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
kpeter@559:                            SetFlowMapTraits<T> > {
alpar@399:       typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
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
alpar@399:       : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
kpeter@559:                            SetElevatorTraits<T> > {
alpar@399:       typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
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@402:     /// 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
alpar@399:       : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
kpeter@559:                        SetStandardElevatorTraits<T> > {
alpar@399:       typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
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: 
alpar@399:     /// The constructor of the class.
alpar@399: 
alpar@399:     /// The constructor of the class.
alpar@399:     /// \param g The digraph the algorithm runs on.
alpar@399:     /// \param lo The lower bound capacity of the arcs.
alpar@399:     /// \param up The upper bound capacity of the arcs.
kpeter@402:     /// \param delta The lower bound for the supply of the nodes.
alpar@399:     Circulation(const Digraph &g,const LCapMap &lo,
alpar@399:                 const UCapMap &up,const DeltaMap &delta)
alpar@399:       : _g(g), _node_num(),
alpar@399:         _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
alpar@399:         _level(0), _local_level(false), _excess(0), _el() {}
alpar@399: 
kpeter@402:     /// Destructor.
alpar@399:     ~Circulation() {
alpar@399:       destroyStructures();
alpar@399:     }
alpar@399: 
kpeter@402: 
alpar@399:   private:
alpar@399: 
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: 
alpar@399:     /// Sets the lower bound capacity map.
alpar@399: 
alpar@399:     /// Sets the lower bound capacity map.
kpeter@402:     /// \return <tt>(*this)</tt>
alpar@399:     Circulation& lowerCapMap(const LCapMap& map) {
alpar@399:       _lo = &map;
alpar@399:       return *this;
alpar@399:     }
alpar@399: 
alpar@399:     /// Sets the upper bound capacity map.
alpar@399: 
alpar@399:     /// Sets the upper bound capacity map.
kpeter@402:     /// \return <tt>(*this)</tt>
alpar@399:     Circulation& upperCapMap(const LCapMap& map) {
alpar@399:       _up = &map;
alpar@399:       return *this;
alpar@399:     }
alpar@399: 
kpeter@402:     /// Sets the lower bound map for the supply of the nodes.
alpar@399: 
kpeter@402:     /// Sets the lower bound map for the supply of the nodes.
kpeter@402:     /// \return <tt>(*this)</tt>
alpar@399:     Circulation& deltaMap(const DeltaMap& map) {
alpar@399:       _delta = &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@402:     /// \brief Sets the tolerance used by algorithm.
kpeter@402:     ///
alpar@399:     /// Sets the tolerance used by algorithm.
alpar@399:     Circulation& tolerance(const Tolerance& tolerance) const {
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@402:     /// Returns a const reference to the tolerance.
alpar@399:     const Tolerance& tolerance() const {
alpar@399:       return tolerance;
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@402:     /// If you need more control on the initial solution or the execution,
kpeter@402:     /// first you have to call one of the \ref init() functions, 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:     {
alpar@399:       createStructures();
alpar@399: 
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n) {
kpeter@581:         (*_excess)[n] = (*_delta)[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:     {
alpar@399:       createStructures();
alpar@399: 
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n) {
kpeter@581:         (*_excess)[n] = (*_delta)[n];
alpar@399:       }
alpar@399: 
alpar@399:       for (ArcIt e(_g);e!=INVALID;++e) {
alpar@399:         if (!_tol.positive((*_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];
alpar@399:         } else if (_tol.positive((*_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 {
alpar@399:           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;
alpar@399:         Value exc=(*_excess)[act];
alpar@399: 
alpar@399:         for(OutArcIt e(_g,act);e!=INVALID; ++e) {
alpar@399:           Node v = _g.target(e);
alpar@399:           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);
alpar@399:           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@402:     /// \brief Returns the flow on the given arc.
kpeter@402:     ///
kpeter@402:     /// Returns the flow 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@402:     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@402:        \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
kpeter@402:            \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(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:         {
alpar@399:           Value dif=-(*_delta)[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:     {
alpar@399:       Value delta=0;
alpar@399:       for(NodeIt n(_g);n!=INVALID;++n)
alpar@399:         if(barrier(n))
alpar@399:           delta-=(*_delta)[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);
alpar@399:           if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
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