RhodeCode

    ///

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

    /// nodes. 

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

    typedef SM SupplyMap;

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

    ///

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

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

    /// concept.

    /// \brief Instantiates a FlowMap.

    ///

    /// This function instantiates a \ref FlowMap.

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

    /// the flow map.

      return new Elevator(digraph, max_level);

    }

    /// \brief The tolerance used by the algorithm

    ///

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

  };

  /**

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

  public:

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

    typedef TR Traits;

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

    typedef typename Traits::Digraph Digraph;

    ///The type of the lower bound map.

    typedef typename Traits::LowerMap LowerMap;

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

    typedef typename Traits::UpperMap UpperMap;

    ///The type of the supply map.

    FlowMap *_flow;

    bool _local_flow;

    Elevator* _level;

    bool _local_level;

    ExcessMap* _excess;

    Tolerance _tol;

    int _el;

  public:

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

        } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {

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

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

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

        } else {

          _flow->set(e, fc);

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

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

        }

      }

      Node act;

      Node bact=INVALID;

      Node last_activated=INVALID;

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

        int actlevel=(*_level)[act];

        int mlevel=_node_num;

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

          Node v = _g.target(e);

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

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

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

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

              (*_excess)[v] += exc;

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

            }

          }

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

        }

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

          Node v = _g.source(e);

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

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

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

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

              (*_excess)[v] += exc;

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

    /// these functions.\n

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

    /// using them.

    ///@{

    ///

    ///

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

    /// using this function.

      return (*_flow)[arc];

    }

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

    ///

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

    ///

    bool checkFlow() const {

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

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

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

        {

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

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

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

        }

      return true;

    }

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

    ///\sa barrier()

    ///\sa barrierMap()

    bool checkBarrier() const

    {

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

        if(barrier(n))

          delta-=(*_supply)[n];

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

        {

          Node s=_g.source(e);

  /// Most of the parameters of the problem (except for the digraph)

  /// can be given using separate functions, and the algorithm can be

  /// executed using the \ref run() function. If some parameters are not

  /// specified, then default values will be used.

  ///

  /// \tparam GR The digraph type the algorithm runs on.

  /// and supply values in the algorithm. By default it is \c int.

  /// \tparam C The value type used for costs and potentials in the

  ///

  /// \warning Both value types must be signed and all input data must

  /// be integer.

  ///

  /// \note %NetworkSimplex provides five different pivot rule

  /// implementations, from which the most efficient one is used

  /// by default. For more information see \ref PivotRule.

  class NetworkSimplex

  {

  public:

    /// The flow type of the algorithm

    /// The cost type of the algorithm

    typedef C Cost;

#ifdef DOXYGEN

    /// The type of the flow map

    /// The type of the potential map

    typedef GR::NodeMap<Cost> PotentialMap;

#else

    /// The type of the flow map

    /// The type of the potential map

    typedef typename GR::template NodeMap<Cost> PotentialMap;

#endif

  public:

    };

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(GR);

    typedef typename GR::template ArcMap<Cost> CostArcMap;

    typedef std::vector<Arc> ArcVector;

    typedef std::vector<Node> NodeVector;

    typedef std::vector<int> IntVector;

    typedef std::vector<bool> BoolVector;

    typedef std::vector<Cost> CostVector;

    // State constants for arcs

    enum ArcStateEnum {

      STATE_UPPER = -1,

      STATE_TREE  =  0,

    // Data related to the underlying digraph

    const GR &_graph;

    int _node_num;

    int _arc_num;

    // Parameters of the problem

    CostArcMap *_pcost;

    bool _pstsup;

    Node _psource, _ptarget;

    SupplyType _stype;

    // Result maps

    FlowMap *_flow_map;

    PotentialMap *_potential_map;

    bool _local_flow;

    bool _local_potential;

    int _root;

    // Temporary data used in the current pivot iteration

    int in_arc, join, u_in, v_in, u_out, v_out;

    int first, second, right, last;

    int stem, par_stem, new_stem;

  public:

    /// \brief Constant for infinite upper bounds (capacities).

    ///

    /// Constant for infinite upper bounds (capacities).

  private:

    // Implementation of the First Eligible pivot rule

    class FirstEligiblePivotRule

    {

      _graph(graph),

      _plower(NULL), _pupper(NULL), _pcost(NULL),

      _psupply(NULL), _pstsup(false), _stype(GEQ),

      _flow_map(NULL), _potential_map(NULL),

      _local_flow(false), _local_potential(false),

      _node_id(graph),

    {

      // Check the value types

        "The flow type of NetworkSimplex must be signed");

      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,

        "The cost type of NetworkSimplex must be signed");

    }

    /// Destructor.

    ///

    /// This function sets the lower bounds on the arcs.

    /// If it is not used before calling \ref run(), the lower bounds

    /// will be set to zero on all arcs.

    ///

    /// \param map An arc map storing the lower bounds.

    /// of the algorithm.

    ///

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

    template <typename LowerMap>

    NetworkSimplex& lowerMap(const LowerMap& map) {

      delete _plower;

      for (ArcIt a(_graph); a != INVALID; ++a) {

        (*_plower)[a] = map[a];

      }

      return *this;

    }

    /// This function sets the upper bounds (capacities) on the arcs.

    /// If it is not used before calling \ref run(), the upper bounds

    /// will be set to \ref INF on all arcs (i.e. the flow value will be

    /// unbounded from above on each arc).

    ///

    /// \param map An arc map storing the upper bounds.

    /// of the algorithm.

    ///

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

    template<typename UpperMap>

    NetworkSimplex& upperMap(const UpperMap& map) {

      delete _pupper;

      for (ArcIt a(_graph); a != INVALID; ++a) {

        (*_pupper)[a] = map[a];

      }

      return *this;

    }

    /// This function sets the supply values of the nodes.

    /// If neither this function nor \ref stSupply() is used before

    /// calling \ref run(), the supply of each node will be set to zero.

    /// (It makes sense only if non-zero lower bounds are given.)

    ///

    /// \param map A node map storing the supply values.

    /// of the algorithm.

    ///

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

    template<typename SupplyMap>

    NetworkSimplex& supplyMap(const SupplyMap& map) {

      delete _psupply;

      _pstsup = false;

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

        (*_psupply)[n] = map[n];

      }

      return *this;

    }

    /// \param s The source node.

    /// \param t The target node.

    /// \param k The required amount of flow from node \c s to node \c t

    /// (i.e. the supply of \c s and the demand of \c t).

    ///

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

      delete _psupply;

      _psupply = NULL;

      _pstsup = true;

      _psource = s;

      _ptarget = t;

      _pstflow = k;

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

    ///

    /// This function returns the flow on the given arc.

    ///

    /// \pre \ref run() must be called before using this function.

      return (*_flow_map)[a];

    }

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

    ///

    /// This function returns a const reference to an arc map storing

        }

      }

      // Remove non-zero lower bounds

      if (_plower) {

        for (int i = 0; i != _arc_num; ++i) {

          if (c > 0) {

            if (_cap[i] < INF) _cap[i] -= c;

            _supply[_source[i]] -= c;

            _supply[_target[i]] += c;

          }

          else if (c < 0) {

      } else {

        first  = _target[in_arc];

        second = _source[in_arc];

      }

      delta = _cap[in_arc];

      int result = 0;

      int e;

      // Search the cycle along the path form the first node to the root

      for (int u = first; u != join; u = _parent[u]) {

        e = _pred[u];

        d = _forward[u] ?

    }

    // Change _flow and _state vectors

    void changeFlow(bool change) {

      // Augment along the cycle

      if (delta > 0) {

        _flow[in_arc] += val;

        for (int u = _source[in_arc]; u != join; u = _parent[u]) {

          _flow[_pred[u]] += _forward[u] ? -val : val;

        }

        for (int u = _target[in_arc]; u != join; u = _parent[u]) {

          _flow[_pred[u]] += _forward[u] ? val : -val;

    ///

    /// The type of the map that stores the arc capacities.

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

    typedef CAP CapacityMap;

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

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

    ///

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

    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.

    /// \brief Instantiates a FlowMap.

    ///

    /// This function instantiates a \ref FlowMap.

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

    /// the flow map.

      return new Elevator(digraph, max_level);

    }

    /// \brief The tolerance used by the algorithm

    ///

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

  };

  /// \ingroup max_flow

  ///

    typedef TR Traits;

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

    typedef typename Traits::Digraph Digraph;

    ///The type of the capacity map.

    typedef typename Traits::CapacityMap CapacityMap;

    ///The type of the flow values.

    ///The type of the flow map.

    typedef typename Traits::FlowMap FlowMap;

    ///The type of the elevator.

    typedef typename Traits::Elevator Elevator;

    ///The type of the tolerance.

    FlowMap* _flow;

    bool _local_flow;

    Elevator* _level;

    bool _local_level;

    ExcessMap* _excess;

    Tolerance _tolerance;

    bool _phase;

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

        _flow->set(e, flowMap[e]);

      }

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

        for (InArcIt e(_graph, n); e != INVALID; ++e) {

          excess += (*_flow)[e];

        }

        for (OutArcIt e(_graph, n); e != INVALID; ++e) {

          excess -= (*_flow)[e];

        }

        }

        queue.swap(nqueue);

      }

      _level->initFinish();

      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {

        if (_tolerance.positive(rem)) {

          Node u = _graph.target(e);

          if ((*_level)[u] == _level->maxLevel()) continue;

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

          (*_excess)[u] += rem;

          if (u != _target && !_level->active(u)) {

            _level->activate(u);

          }

        }

      }

      for (InArcIt e(_graph, _source); e != INVALID; ++e) {

        if (_tolerance.positive(rem)) {

          Node v = _graph.source(e);

          if ((*_level)[v] == _level->maxLevel()) continue;

          _flow->set(e, 0);

          (*_excess)[v] += rem;

          if (v != _target && !_level->active(v)) {

      Node n = _level->highestActive();

      int level = _level->highestActiveLevel();

      while (n != INVALID) {

        int num = _node_num;

        while (num > 0 && n != INVALID) {

          int new_level = _level->maxLevel();

          for (OutArcIt e(_graph, n); e != INVALID; ++e) {

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.target(e);

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

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

            } else if (new_level > (*_level)[v]) {

              new_level = (*_level)[v];

            }

          }

          for (InArcIt e(_graph, n); e != INVALID; ++e) {

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.source(e);

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

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

          level = _level->highestActiveLevel();

          --num;

        }

        num = _node_num * 20;

        while (num > 0 && n != INVALID) {

          int new_level = _level->maxLevel();

          for (OutArcIt e(_graph, n); e != INVALID; ++e) {

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.target(e);

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

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

            } else if (new_level > (*_level)[v]) {

              new_level = (*_level)[v];

            }

          }

          for (InArcIt e(_graph, n); e != INVALID; ++e) {

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.source(e);

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

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

          _level->activate(n);

        }

      }

      Node n;

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

        int level = _level->highestActiveLevel();

        int new_level = _level->maxLevel();

        for (OutArcIt e(_graph, n); e != INVALID; ++e) {

          if (!_tolerance.positive(rem)) continue;

          Node v = _graph.target(e);

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

            if (!_level->active(v) && v != _source) {

              _level->activate(v);

            }

          } else if (new_level > (*_level)[v]) {

            new_level = (*_level)[v];

          }

        }

        for (InArcIt e(_graph, n); e != INVALID; ++e) {

          if (!_tolerance.positive(rem)) continue;

          Node v = _graph.source(e);

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

            if (!_level->active(v) && v != _source) {

              _level->activate(v);

            }

    /// Returns the value of the maximum flow by returning the excess

    /// of the target node. This value equals to the value of

    /// the maximum flow already after the first phase of the algorithm.

    ///

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

    /// using this function.

      return (*_excess)[_target];

    }

    ///

    /// be called after the second phase of the algorithm.

    ///

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

    /// using this function.

      return (*_flow)[arc];

    }

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

    ///

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