RhodeCode

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Copyright (C) 2003-2010

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_CAPACITY_SCALING_H

#define LEMON_CAPACITY_SCALING_H

/// \ingroup min_cost_flow_algs

///

/// \file

/// \brief Capacity Scaling algorithm for finding a minimum cost flow.

#include <vector>

#include <limits>

#include <lemon/core.h>

#include <lemon/bin_heap.h>

namespace lemon {

  /// \brief Default traits class of CapacityScaling algorithm.

  ///

  /// Default traits class of CapacityScaling algorithm.

  /// \tparam GR Digraph type.

  /// \tparam V The number type used for flow amounts, capacity bounds

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

  /// \tparam C The number type used for costs and potentials.

  /// By default it is the same as \c V.

  template <typename GR, typename V = int, typename C = V>

  struct CapacityScalingDefaultTraits

  {

    /// The type of the digraph

    typedef GR Digraph;

    /// The type of the flow amounts, capacity bounds and supply values

    typedef V Value;

    /// The type of the arc costs

    typedef C Cost;

    /// \brief The type of the heap used for internal Dijkstra computations.

    ///

    /// The type of the heap used for internal Dijkstra computations.

    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,

    /// its priority type must be \c Cost and its cross reference type

    /// must be \ref RangeMap "RangeMap<int>".

    typedef BinHeap<Cost, RangeMap<int> > Heap;

  };

  /// \addtogroup min_cost_flow_algs

  /// @{

  /// \brief Implementation of the Capacity Scaling algorithm for

  /// finding a \ref min_cost_flow "minimum cost flow".

  ///

  /// \ref CapacityScaling implements the capacity scaling version

  /// of the successive shortest path algorithm for finding a

  /// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,

  /// \ref edmondskarp72theoretical. It is an efficient dual

  /// solution method.

  ///

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

  /// \tparam V The number type used for flow amounts, capacity bounds

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

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

  /// algorithm. By default, it is the same as \c V.

  /// \tparam TR The traits class that defines various types used by the

  /// algorithm. By default, it is \ref CapacityScalingDefaultTraits

  /// "CapacityScalingDefaultTraits<GR, V, C>".

  /// In most cases, this parameter should not be set directly,

  /// consider to use the named template parameters instead.

  ///

  /// \warning Both \c V and \c C must be signed number types.

  /// \warning This algorithm does not support negative costs for

  /// arcs having infinite upper bound.

#ifdef DOXYGEN

  template <typename GR, typename V, typename C, typename TR>

#else

  template < typename GR, typename V = int, typename C = V,

             typename TR = CapacityScalingDefaultTraits<GR, V, C> >

#endif

  class CapacityScaling

  {

  public:

    /// The type of the digraph

    typedef typename TR::Digraph Digraph;

    /// The type of the flow amounts, capacity bounds and supply values

    typedef typename TR::Value Value;

    /// The type of the arc costs

    typedef typename TR::Cost Cost;

    /// The type of the heap used for internal Dijkstra computations

    typedef typename TR::Heap Heap;

    /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm

    typedef TR Traits;

  public:

    /// \brief Problem type constants for the \c run() function.

    ///

    /// Enum type containing the problem type constants that can be

    /// returned by the \ref run() function of the algorithm.

    enum ProblemType {

      /// The problem has no feasible solution (flow).

      INFEASIBLE,

      /// The problem has optimal solution (i.e. it is feasible and

      /// bounded), and the algorithm has found optimal flow and node

      /// potentials (primal and dual solutions).

      OPTIMAL,

      /// The digraph contains an arc of negative cost and infinite

      /// upper bound. It means that the objective function is unbounded

      /// on that arc, however, note that it could actually be bounded

      /// over the feasible flows, but this algroithm cannot handle

      /// these cases.

      UNBOUNDED

    };

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(GR);

    typedef std::vector<int> IntVector;

    typedef std::vector<Value> ValueVector;

    typedef std::vector<Cost> CostVector;

    typedef std::vector<char> BoolVector;

    // Note: vector<char> is used instead of vector<bool> for efficiency reasons

  private:

    // Data related to the underlying digraph

    const GR &_graph;

    int _node_num;

    int _arc_num;

    int _res_arc_num;

    int _root;

    // Parameters of the problem

    bool _have_lower;

    Value _sum_supply;

    // Data structures for storing the digraph

    IntNodeMap _node_id;

    IntArcMap _arc_idf;

    IntArcMap _arc_idb;

    IntVector _first_out;

    BoolVector _forward;

    IntVector _source;

    IntVector _target;

    IntVector _reverse;

    // Node and arc data

    ValueVector _lower;

    ValueVector _upper;

    CostVector _cost;

    ValueVector _supply;

    ValueVector _res_cap;

    CostVector _pi;

    ValueVector _excess;

    IntVector _excess_nodes;

    IntVector _deficit_nodes;

    Value _delta;

    int _factor;

    IntVector _pred;

  public:

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

    ///

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

    /// It is \c std::numeric_limits<Value>::infinity() if available,

    /// \c std::numeric_limits<Value>::max() otherwise.

    const Value INF;

  private:

    // Special implementation of the Dijkstra algorithm for finding

    // shortest paths in the residual network of the digraph with

    // respect to the reduced arc costs and modifying the node

    // potentials according to the found distance labels.

    class ResidualDijkstra

    {

    private:

      int _node_num;

      bool _geq;

      const IntVector &_first_out;

      const IntVector &_target;

      const CostVector &_cost;

      const ValueVector &_res_cap;

      const ValueVector &_excess;

      CostVector &_pi;

      IntVector &_pred;

      IntVector _proc_nodes;

      CostVector _dist;

    public:

      ResidualDijkstra(CapacityScaling& cs) :

        _node_num(cs._node_num), _geq(cs._sum_supply < 0),

        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),

        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),

        _pred(cs._pred), _dist(cs._node_num)

      {}

      int run(int s, Value delta = 1) {

        RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);

        Heap heap(heap_cross_ref);

        heap.push(s, 0);

        _pred[s] = -1;

        _proc_nodes.clear();

        // Process nodes

        while (!heap.empty() && _excess[heap.top()] > -delta) {

          int u = heap.top(), v;

          Cost d = heap.prio() + _pi[u], dn;

          _dist[u] = heap.prio();

          _proc_nodes.push_back(u);

          heap.pop();

          // Traverse outgoing residual arcs

          int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;

          for (int a = _first_out[u]; a != last_out; ++a) {

            if (_res_cap[a] < delta) continue;

            v = _target[a];

            switch (heap.state(v)) {

              case Heap::PRE_HEAP:

                heap.push(v, d + _cost[a] - _pi[v]);

                _pred[v] = a;

                break;

              case Heap::IN_HEAP:

                dn = d + _cost[a] - _pi[v];

                if (dn < heap[v]) {

                  heap.decrease(v, dn);

                  _pred[v] = a;

                }

                break;

              case Heap::POST_HEAP:

                break;

            }

          }

        }

        if (heap.empty()) return -1;

        // Update potentials of processed nodes

        int t = heap.top();

        Cost dt = heap.prio();

        for (int i = 0; i < int(_proc_nodes.size()); ++i) {

          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;

        }

        return t;

      }

    }; //class ResidualDijkstra

  public:

    /// \name Named Template Parameters

    /// @{

    template <typename T>

    struct SetHeapTraits : public Traits {

      typedef T Heap;

    };

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

    /// \c Heap type.

    ///

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

    /// type, which is used for internal Dijkstra computations.

    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,

    /// its priority type must be \c Cost and its cross reference type

    /// must be \ref RangeMap "RangeMap<int>".

    template <typename T>

    struct SetHeap

      : public CapacityScaling<GR, V, C, SetHeapTraits<T> > {

      typedef  CapacityScaling<GR, V, C, SetHeapTraits<T> > Create;

    };

    /// @}

  protected:

    CapacityScaling() {}

  public:

    /// \brief Constructor.

    ///

    /// The constructor of the class.

    ///

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

    CapacityScaling(const GR& graph) :

      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),

      INF(std::numeric_limits<Value>::has_infinity ?

          std::numeric_limits<Value>::infinity() :

          std::numeric_limits<Value>::max())

    {

      // Check the number types

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

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

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

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

      // Reset data structures

      reset();

    }

    /// \name Parameters

    /// The parameters of the algorithm can be specified using these

    /// functions.

    /// @{

    /// \brief Set the lower bounds on the arcs.

    ///

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

    /// Its \c Value type must be convertible to the \c Value type

    /// of the algorithm.

    ///

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

    template <typename LowerMap>

    CapacityScaling& lowerMap(const LowerMap& map) {

      _have_lower = true;

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

        _lower[_arc_idf[a]] = map[a];

        _lower[_arc_idb[a]] = map[a];

      }

      return *this;

    }

    /// \brief Set the upper bounds (capacities) on the arcs.

    ///

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

    ///

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

    /// Its \c Value type must be convertible to the \c Value type

    /// of the algorithm.

    ///

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

    template<typename UpperMap>

    CapacityScaling& upperMap(const UpperMap& map) {

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

        _upper[_arc_idf[a]] = map[a];

      }

      return *this;

    }

    /// \brief Set the costs of the arcs.

    ///

    /// This function sets the costs of the arcs.

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

    /// will be set to \c 1 on all arcs.

    ///

    /// \param map An arc map storing the costs.

    /// Its \c Value type must be convertible to the \c Cost type

    /// of the algorithm.

    ///

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

    template<typename CostMap>

    CapacityScaling& costMap(const CostMap& map) {

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

        _cost[_arc_idf[a]] =  map[a];

        _cost[_arc_idb[a]] = -map[a];

      }

      return *this;

    }

    /// \brief Set the supply values of the nodes.

    ///

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

    ///

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

    /// Its \c Value type must be convertible to the \c Value type

    /// of the algorithm.

    ///

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

    template<typename SupplyMap>

    CapacityScaling& supplyMap(const SupplyMap& map) {

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

        _supply[_node_id[n]] = map[n];

      }

      return *this;

    }

    /// \brief Set single source and target nodes and a supply value.

    ///

    /// This function sets a single source node and a single target node

    /// and the required flow value.

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

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

    ///

    /// Using this function has the same effect as using \ref supplyMap()

    /// with a map in which \c k is assigned to \c s, \c -k is

    /// assigned to \c t and all other nodes have zero supply value.

    ///

    /// \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>

    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {

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

        _supply[i] = 0;

      }

      _supply[_node_id[s]] =  k;

      _supply[_node_id[t]] = -k;

      return *this;

    }

    /// @}

    /// \name Execution control

    /// The algorithm can be executed using \ref run().

    /// @{

    /// \brief Run the algorithm.

    ///

    /// This function runs the algorithm.

    /// The paramters can be specified using functions \ref lowerMap(),

    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().

    /// For example,

    /// \code

    ///   CapacityScaling<ListDigraph> cs(graph);

    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)

    ///     .supplyMap(sup).run();

    /// \endcode

    ///

    /// This function can be called more than once. All the given parameters

    /// are kept for the next call, unless \ref resetParams() or \ref reset()

    /// is used, thus only the modified parameters have to be set again.

    /// If the underlying digraph was also modified after the construction

    /// of the class (or the last \ref reset() call), then the \ref reset()

    /// function must be called.

    ///

    /// \param factor The capacity scaling factor. It must be larger than

    /// one to use scaling. If it is less or equal to one, then scaling

    /// will be disabled.

    ///

    /// \return \c INFEASIBLE if no feasible flow exists,

    /// \n \c OPTIMAL if the problem has optimal solution

    /// (i.e. it is feasible and bounded), and the algorithm has found

    /// optimal flow and node potentials (primal and dual solutions),

    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost

    /// and infinite upper bound. It means that the objective function

    /// is unbounded on that arc, however, note that it could actually be

    /// bounded over the feasible flows, but this algroithm cannot handle

    /// these cases.

    ///

    /// \see ProblemType

    /// \see resetParams(), reset()

    ProblemType run(int factor = 4) {

      _factor = factor;

      ProblemType pt = init();

      if (pt != OPTIMAL) return pt;

      return start();

    }

    /// \brief Reset all the parameters that have been given before.

    ///

    /// This function resets all the paramaters that have been given

    /// before using functions \ref lowerMap(), \ref upperMap(),

    /// \ref costMap(), \ref supplyMap(), \ref stSupply().

    ///

    /// It is useful for multiple \ref run() calls. Basically, all the given

    /// parameters are kept for the next \ref run() call, unless

    /// \ref resetParams() or \ref reset() is used.

    /// If the underlying digraph was also modified after the construction

    /// of the class or the last \ref reset() call, then the \ref reset()

    /// function must be used, otherwise \ref resetParams() is sufficient.

    ///

    /// For example,

    /// \code

    ///   CapacityScaling<ListDigraph> cs(graph);

    ///

    ///   // First run

    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)

    ///     .supplyMap(sup).run();

    ///

    ///   // Run again with modified cost map (resetParams() is not called,

    ///   // so only the cost map have to be set again)

    ///   cost[e] += 100;

    ///   cs.costMap(cost).run();

    ///

    ///   // Run again from scratch using resetParams()

    ///   // (the lower bounds will be set to zero on all arcs)

    ///   cs.resetParams();

    ///   cs.upperMap(capacity).costMap(cost)

    ///     .supplyMap(sup).run();

    /// \endcode

    ///

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

    ///

    /// \see reset(), run()

    CapacityScaling& resetParams() {

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

        _supply[i] = 0;

      }

      for (int j = 0; j != _res_arc_num; ++j) {

        _lower[j] = 0;

        _upper[j] = INF;

        _cost[j] = _forward[j] ? 1 : -1;

      }

      _have_lower = false;

      return *this;

    }

    /// \brief Reset the internal data structures and all the parameters

    /// that have been given before.

    ///

    /// This function resets the internal data structures and all the

    /// paramaters that have been given before using functions \ref lowerMap(),

    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().

    ///

    /// It is useful for multiple \ref run() calls. Basically, all the given

    /// parameters are kept for the next \ref run() call, unless

    /// \ref resetParams() or \ref reset() is used.

    /// If the underlying digraph was also modified after the construction

    /// of the class or the last \ref reset() call, then the \ref reset()

    /// function must be used, otherwise \ref resetParams() is sufficient.

    ///

    /// See \ref resetParams() for examples.

    ///

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

    ///

    /// \see resetParams(), run()

    CapacityScaling& reset() {

      // Resize vectors

      _node_num = countNodes(_graph);

      _arc_num = countArcs(_graph);

      _res_arc_num = 2 * (_arc_num + _node_num);

      _root = _node_num;

      ++_node_num;

      _first_out.resize(_node_num + 1);

      _forward.resize(_res_arc_num);

      _source.resize(_res_arc_num);

      _target.resize(_res_arc_num);

      _reverse.resize(_res_arc_num);

      _lower.resize(_res_arc_num);

      _upper.resize(_res_arc_num);

      _cost.resize(_res_arc_num);

      _supply.resize(_node_num);

      _res_cap.resize(_res_arc_num);

      _pi.resize(_node_num);

      _excess.resize(_node_num);

      _pred.resize(_node_num);

      // Copy the graph

      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;

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

        _node_id[n] = i;

      }

      i = 0;

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

        _first_out[i] = j;

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

          _arc_idf[a] = j;

          _forward[j] = true;

          _source[j] = i;

          _target[j] = _node_id[_graph.runningNode(a)];

        }

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

          _arc_idb[a] = j;

          _forward[j] = false;

          _source[j] = i;

          _target[j] = _node_id[_graph.runningNode(a)];

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Copyright (C) 2003-2010

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_NETWORK_SIMPLEX_H

#define LEMON_NETWORK_SIMPLEX_H

/// \ingroup min_cost_flow_algs

///

/// \file

/// \brief Network Simplex algorithm for finding a minimum cost flow.

#include <vector>

#include <limits>

#include <algorithm>

#include <lemon/core.h>

#include <lemon/math.h>

namespace lemon {

  /// \addtogroup min_cost_flow_algs

  /// @{

  /// \brief Implementation of the primal Network Simplex algorithm

  /// for finding a \ref min_cost_flow "minimum cost flow".

  ///

  /// \ref NetworkSimplex implements the primal Network Simplex algorithm

  /// for finding a \ref min_cost_flow "minimum cost flow"

  /// \ref amo93networkflows, \ref dantzig63linearprog,

  /// \ref kellyoneill91netsimplex.

  /// This algorithm is a highly efficient specialized version of the

  /// linear programming simplex method directly for the minimum cost

  /// flow problem.

  ///

  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest

  /// implementations available in LEMON for this problem.

  /// Furthermore, this class supports both directions of the supply/demand

  /// inequality constraints. For more information, see \ref SupplyType.

  ///

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

  /// \tparam V The number type used for flow amounts, capacity bounds

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

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

  /// algorithm. By default, it is the same as \c V.

  ///

  /// \warning Both \c V and \c C must be signed number types.

  /// \warning All input data (capacities, supply values, and costs) 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.

  template <typename GR, typename V = int, typename C = V>

  class NetworkSimplex

  {

  public:

    /// The type of the flow amounts, capacity bounds and supply values

    typedef V Value;

    /// The type of the arc costs

    typedef C Cost;

  public:

    /// \brief Problem type constants for the \c run() function.

    ///

    /// Enum type containing the problem type constants that can be

    /// returned by the \ref run() function of the algorithm.

    enum ProblemType {

      /// The problem has no feasible solution (flow).

      INFEASIBLE,

      /// The problem has optimal solution (i.e. it is feasible and

      /// bounded), and the algorithm has found optimal flow and node

      /// potentials (primal and dual solutions).

      OPTIMAL,

      /// The objective function of the problem is unbounded, i.e.

      /// there is a directed cycle having negative total cost and

      /// infinite upper bound.

      UNBOUNDED

    };

    /// \brief Constants for selecting the type of the supply constraints.

    ///

    /// Enum type containing constants for selecting the supply type,

    /// i.e. the direction of the inequalities in the supply/demand

    /// constraints of the \ref min_cost_flow "minimum cost flow problem".

    ///

    /// The default supply type is \c GEQ, the \c LEQ type can be

    /// selected using \ref supplyType().

    /// The equality form is a special case of both supply types.

    enum SupplyType {

      /// This option means that there are <em>"greater or equal"</em>

      /// supply/demand constraints in the definition of the problem.

      GEQ,

      /// This option means that there are <em>"less or equal"</em>

      /// supply/demand constraints in the definition of the problem.

      LEQ

    };

    /// \brief Constants for selecting the pivot rule.

    ///

    /// Enum type containing constants for selecting the pivot rule for

    /// the \ref run() function.

    ///

    /// of the algorithm.

    enum PivotRule {

      /// The \e First \e Eligible pivot rule.

      /// The next eligible arc is selected in a wraparound fashion

      /// in every iteration.

      FIRST_ELIGIBLE,

      /// The \e Best \e Eligible pivot rule.

      /// The best eligible arc is selected in every iteration.

      BEST_ELIGIBLE,

      /// The \e Block \e Search pivot rule.

      /// A specified number of arcs are examined in every iteration

      /// in a wraparound fashion and the best eligible arc is selected

      /// from this block.

      BLOCK_SEARCH,

      /// The \e Candidate \e List pivot rule.

      /// In a major iteration a candidate list is built from eligible arcs

      /// in a wraparound fashion and in the following minor iterations

      /// the best eligible arc is selected from this list.

      CANDIDATE_LIST,

      /// The \e Altering \e Candidate \e List pivot rule.

      /// It is a modified version of the Candidate List method.

      /// candidate list and extends this list in every iteration.

      ALTERING_LIST

    };

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(GR);

    typedef std::vector<int> IntVector;

    typedef std::vector<Value> ValueVector;

    typedef std::vector<Cost> CostVector;

    typedef std::vector<signed char> CharVector;

    // Note: vector<signed char> is used instead of vector<ArcState> and

    // vector<ArcDirection> for efficiency reasons

    // State constants for arcs

    enum ArcState {

      STATE_UPPER = -1,

      STATE_TREE  =  0,

      STATE_LOWER =  1

    };

    // Direction constants for tree arcs

    enum ArcDirection {

      DIR_DOWN = -1,

      DIR_UP   =  1

    };

  private:

    // Data related to the underlying digraph

    const GR &_graph;

    int _node_num;

    int _arc_num;

    int _all_arc_num;

    int _search_arc_num;

    // Parameters of the problem

    bool _have_lower;

    SupplyType _stype;

    Value _sum_supply;

    // Data structures for storing the digraph

    IntNodeMap _node_id;

    IntArcMap _arc_id;

    IntVector _source;

    IntVector _target;

    bool _arc_mixing;

    // Node and arc data

    ValueVector _lower;

    ValueVector _upper;

    ValueVector _cap;

    CostVector _cost;

    ValueVector _supply;

    ValueVector _flow;

    CostVector _pi;

    // Data for storing the spanning tree structure

    IntVector _parent;

    IntVector _pred;

    IntVector _thread;

    IntVector _rev_thread;

    IntVector _succ_num;

    IntVector _last_succ;

    CharVector _pred_dir;

    CharVector _state;

    IntVector _dirty_revs;

    int _root;

    // Temporary data used in the current pivot iteration

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

    Value delta;

    const Value MAX;

  public:

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

    ///

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

    /// It is \c std::numeric_limits<Value>::infinity() if available,

    /// \c std::numeric_limits<Value>::max() otherwise.

    const Value INF;

  private:

    // Implementation of the First Eligible pivot rule

    class FirstEligiblePivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const CharVector &_state;

      const CostVector &_pi;

      int &_in_arc;

      int _search_arc_num;

      // Pivot rule data

      int _next_arc;

    public:

      // Constructor

      FirstEligiblePivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),

        _next_arc(0)

      {}

      // Find next entering arc

      bool findEnteringArc() {

        Cost c;

        for (int e = _next_arc; e != _search_arc_num; ++e) {

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < 0) {

            _in_arc = e;

            _next_arc = e + 1;

            return true;

          }

        }

        for (int e = 0; e != _next_arc; ++e) {

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < 0) {

            _in_arc = e;

            _next_arc = e + 1;

            return true;

          }

        }

        return false;

      }

    }; //class FirstEligiblePivotRule

    // Implementation of the Best Eligible pivot rule

    class BestEligiblePivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const CharVector &_state;

      const CostVector &_pi;

      int &_in_arc;

      int _search_arc_num;

    public:

      // Constructor

      BestEligiblePivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)

      {}

      // Find next entering arc

      bool findEnteringArc() {

        Cost c, min = 0;

        for (int e = 0; e != _search_arc_num; ++e) {

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < min) {

            min = c;

            _in_arc = e;

          }

        }

        return min < 0;

      }

    }; //class BestEligiblePivotRule

    // Implementation of the Block Search pivot rule

    class BlockSearchPivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const CharVector &_state;

      const CostVector &_pi;

      int &_in_arc;

      int _search_arc_num;

      // Pivot rule data

      int _block_size;

      int _next_arc;

    public:

      // Constructor

      BlockSearchPivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),

        _next_arc(0)

      {

        // The main parameters of the pivot rule

        const double BLOCK_SIZE_FACTOR = 1.0;

        const int MIN_BLOCK_SIZE = 10;

        _block_size = std::max( int(BLOCK_SIZE_FACTOR *

                                    std::sqrt(double(_search_arc_num))),

                                MIN_BLOCK_SIZE );

      }

      // Find next entering arc

      bool findEnteringArc() {

        Cost c, min = 0;

        int cnt = _block_size;

        int e;

        for (e = _next_arc; e != _search_arc_num; ++e) {

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < min) {

            min = c;

            _in_arc = e;

          }

          if (--cnt == 0) {

            if (min < 0) goto search_end;

            cnt = _block_size;

          }

        }

        for (e = 0; e != _next_arc; ++e) {

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < min) {

            min = c;

            _in_arc = e;

          }

          if (--cnt == 0) {

            if (min < 0) goto search_end;

            cnt = _block_size;

          }

        }

        if (min >= 0) return false;

      search_end:

        _next_arc = e;

        return true;

      }

    }; //class BlockSearchPivotRule

    // Implementation of the Candidate List pivot rule

    class CandidateListPivotRule

    {

    private:

      // References to the NetworkSimplex class

      const IntVector  &_source;

      const IntVector  &_target;

      const CostVector &_cost;

      const CharVector &_state;

      const CostVector &_pi;

      int &_in_arc;

      int _search_arc_num;

      // Pivot rule data

      IntVector _candidates;

      int _list_length, _minor_limit;

      int _curr_length, _minor_count;

      int _next_arc;

    public:

      /// Constructor

      CandidateListPivotRule(NetworkSimplex &ns) :

        _source(ns._source), _target(ns._target),

        _cost(ns._cost), _state(ns._state), _pi(ns._pi),

        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),

        _next_arc(0)

      {

        // The main parameters of the pivot rule

        const double LIST_LENGTH_FACTOR = 0.25;

        const int MIN_LIST_LENGTH = 10;

        const double MINOR_LIMIT_FACTOR = 0.1;

        const int MIN_MINOR_LIMIT = 3;

        _list_length = std::max( int(LIST_LENGTH_FACTOR *

                                     std::sqrt(double(_search_arc_num))),

                                 MIN_LIST_LENGTH );

        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),

                                 MIN_MINOR_LIMIT );

        _curr_length = _minor_count = 0;

        _candidates.resize(_list_length);

      }

      /// Find next entering arc

      bool findEnteringArc() {

        Cost min, c;

        int e;

        if (_curr_length > 0 && _minor_count < _minor_limit) {

          // Minor iteration: select the best eligible arc from the

          // current candidate list

          ++_minor_count;

          min = 0;

          for (int i = 0; i < _curr_length; ++i) {

            e = _candidates[i];

            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

            if (c < min) {

              min = c;

              _in_arc = e;

            }

            else if (c >= 0) {

              _candidates[i--] = _candidates[--_curr_length];

            }

          }

          if (min < 0) return true;

        }

        // Major iteration: build a new candidate list

        min = 0;

        _curr_length = 0;

        for (e = _next_arc; e != _search_arc_num; ++e) {

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < 0) {

            _candidates[_curr_length++] = e;

            if (c < min) {

              min = c;

              _in_arc = e;

            }

            if (_curr_length == _list_length) goto search_end;

          }

        }

        for (e = 0; e != _next_arc; ++e) {

          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (c < 0) {

            _candidates[_curr_length++] = e;

            if (c < min) {

              min = c;