RhodeCode

      SortFunc _sort_func;

    public:

      // Constructor

      AlteringListPivotRule(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), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)

      {

        // The main parameters of the pivot rule

        const double BLOCK_SIZE_FACTOR = 1.0;

        const int MIN_BLOCK_SIZE = 10;

        const double HEAD_LENGTH_FACTOR = 0.1;

        const int MIN_HEAD_LENGTH = 3;

        _block_size = std::max( int(BLOCK_SIZE_FACTOR *

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

                                MIN_BLOCK_SIZE );

        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),

                                 MIN_HEAD_LENGTH );

        _candidates.resize(_head_length + _block_size);

        _curr_length = 0;

      }

      // Find next entering arc

      bool findEnteringArc() {

        // Check the current candidate list

        int e;

        Cost c;

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

          e = _candidates[i];

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

          if (c < 0) {

            _cand_cost[e] = c;

          } else {

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

          }

        }

        // Extend the list

        int cnt = _block_size;

        int limit = _head_length;

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

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

          if (c < 0) {

            _cand_cost[e] = c;

            _candidates[_curr_length++] = e;

          }

          if (--cnt == 0) {

            if (_curr_length > limit) goto search_end;

            limit = 0;

            cnt = _block_size;

          }

        }

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

          _cand_cost[e] = _state[e] *

            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);

          if (_cand_cost[e] < 0) {

            _candidates[_curr_length++] = e;

          }

          if (--cnt == 0) {

            if (_curr_length > limit) goto search_end;

            limit = 0;

            cnt = _block_size;

          }

        }

        if (_curr_length == 0) return false;

      search_end:

        // Make heap of the candidate list (approximating a partial sort)

        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,

                   _sort_func );

        // Pop the first element of the heap

        _in_arc = _candidates[0];

        _next_arc = e;

        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,

                  _sort_func );

        _curr_length = std::min(_head_length, _curr_length - 1);

        return true;

      }

    }; //class AlteringListPivotRule

  public:

    /// \brief Constructor.

    ///

    /// The constructor of the class.

    ///

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

    /// mixed order in the internal data structure.

      _graph(graph), _node_id(graph), _arc_id(graph),

      _arc_mixing(arc_mixing),

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

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

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

    {

      // Check the number types

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

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

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

        "The cost type of NetworkSimplex 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>

    NetworkSimplex& lowerMap(const LowerMap& map) {

      _have_lower = true;

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

        _lower[_arc_id[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>

    NetworkSimplex& upperMap(const UpperMap& map) {

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

        _upper[_arc_id[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>

    NetworkSimplex& costMap(const CostMap& map) {

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

        _cost[_arc_id[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>

    NetworkSimplex& supplyMap(const SupplyMap& map) {

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

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

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

    ///   NetworkSimplex<ListDigraph> ns(graph);

    ///

    ///   // First run

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

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

    ///

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

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

    ///   ns.resetParams();

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

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

    /// \endcode

    ///

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

    ///

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

    NetworkSimplex& resetParams() {

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

        _supply[i] = 0;

      }

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

        _lower[i] = 0;

        _upper[i] = INF;

        _cost[i] = 1;

      }

      _have_lower = false;

      _stype = GEQ;

      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(),

    /// \ref supplyType().

    ///

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

    NetworkSimplex& reset() {

      // Resize vectors

      _node_num = countNodes(_graph);

      _arc_num = countArcs(_graph);

      int all_node_num = _node_num + 1;

      int max_arc_num = _arc_num + 2 * _node_num;

      _source.resize(max_arc_num);

      _target.resize(max_arc_num);

      _lower.resize(_arc_num);

      _upper.resize(_arc_num);

      _cap.resize(max_arc_num);

      _cost.resize(max_arc_num);

      _supply.resize(all_node_num);

      _flow.resize(max_arc_num);

      _pi.resize(all_node_num);

      _parent.resize(all_node_num);

      _pred.resize(all_node_num);

      _pred_dir.resize(all_node_num);

      _thread.resize(all_node_num);

      _rev_thread.resize(all_node_num);

      _succ_num.resize(all_node_num);

      _last_succ.resize(all_node_num);

      _state.resize(max_arc_num);

      // Copy the graph

      int i = 0;

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

        _node_id[n] = i;

      }

      if (_arc_mixing) {

        // Store the arcs in a mixed order

        int i = 0, j = 0;

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

          _arc_id[a] = i;

          _source[i] = _node_id[_graph.source(a)];

          _target[i] = _node_id[_graph.target(a)];

        }

      } else {

        // Store the arcs in the original order

        int i = 0;

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

          _arc_id[a] = i;

          _source[i] = _node_id[_graph.source(a)];

          _target[i] = _node_id[_graph.target(a)];

        }

      }

      // Reset parameters

      resetParams();

      return *this;

    }

    /// @}

    /// \name Query Functions

    /// The results of the algorithm can be obtained using these

    /// functions.\n

    /// The \ref run() function must be called before using them.

    /// @{

    /// \brief Return the total cost of the found flow.

    ///

    /// This function returns the total cost of the found flow.

    /// Its complexity is O(e).

    ///

    /// \note The return type of the function can be specified as a

    /// template parameter. For example,

    /// \code

    ///   ns.totalCost<double>();

    /// \endcode

    /// It is useful if the total cost cannot be stored in the \c Cost

    /// type of the algorithm, which is the default return type of the

    /// function.

    ///

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

    template <typename Number>

    Number totalCost() const {

      Number c = 0;

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

        int i = _arc_id[a];

        c += Number(_flow[i]) * Number(_cost[i]);

      }

      return c;

    }

#ifndef DOXYGEN

    Cost totalCost() const {

      return totalCost<Cost>();

    }

#endif

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

    Value flow(const Arc& a) const {

      return _flow[_arc_id[a]];

    }

    /// \brief Return the flow map (the primal solution).

    ///

    /// This function copies the flow value on each arc into the given

    /// map. The \c Value type of the algorithm must be convertible to

    /// the \c Value type of the map.

    ///

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

    template <typename FlowMap>

    void flowMap(FlowMap &map) const {

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

        map.set(a, _flow[_arc_id[a]]);

      }

    }

    /// \brief Return the potential (dual value) of the given node.

    ///

    /// This function returns the potential (dual value) of the

    /// given node.

    ///

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

    Cost potential(const Node& n) const {

      return _pi[_node_id[n]];

    }

    /// \brief Return the potential map (the dual solution).

    ///

    /// This function copies the potential (dual value) of each node

    /// into the given map.

    /// The \c Cost type of the algorithm must be convertible to the

    /// \c Value type of the map.

    ///