RhodeCode

      PUSH,

      /// Augment operations are used, i.e. flow is moved on admissible

      /// paths from a node with excess to a node with deficit.

      AUGMENT,

      /// Partial augment operations are used, i.e. flow is moved on

      /// admissible paths started from a node with excess, but the

      /// lengths of these paths are limited. This method can be viewed

      /// as a combined version of the previous two operations.

      PARTIAL_AUGMENT

    };

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(GR);

    typedef std::vector<int> IntVector;

    typedef std::vector<Value> ValueVector;

    typedef std::vector<Cost> CostVector;

    typedef std::vector<LargeCost> LargeCostVector;

    typedef std::vector<char> BoolVector;

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

  private:

    template <typename KT, typename VT>

    class StaticVectorMap {

    public:

      typedef KT Key;

      typedef VT Value;

      StaticVectorMap(std::vector<Value>& v) : _v(v) {}

      const Value& operator[](const Key& key) const {

        return _v[StaticDigraph::id(key)];

      }

      Value& operator[](const Key& key) {

        return _v[StaticDigraph::id(key)];

      }

      void set(const Key& key, const Value& val) {

        _v[StaticDigraph::id(key)] = val;

      }

    private:

      std::vector<Value>& _v;

    };

    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;

  private:

    // Data related to the underlying digraph

    const GR &_graph;

    int _node_num;

    int _arc_num;

    int _res_node_num;

    int _res_arc_num;

    int _root;

    // Parameters of the problem

    bool _have_lower;

    Value _sum_supply;

    int _sup_node_num;

    // 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 _scost;

    ValueVector _supply;

    ValueVector _res_cap;

    LargeCostVector _cost;

    LargeCostVector _pi;

    ValueVector _excess;

    IntVector _next_out;

    std::deque<int> _active_nodes;

    // Data for scaling

    LargeCost _epsilon;

    int _alpha;

    IntVector _buckets;

    IntVector _bucket_next;

    IntVector _bucket_prev;

    IntVector _rank;

    int _max_rank;

  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;

  public:

    /// \name Named Template Parameters

    /// @{

    template <typename T>

    struct SetLargeCostTraits : public Traits {

      typedef T LargeCost;

    };

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

    /// \c LargeCost type.

    ///

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

    /// type, which is used for internal computations in the algorithm.

    /// \c Cost must be convertible to \c LargeCost.

    template <typename T>

    struct SetLargeCost

      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {

      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;

    };

    /// @}

  protected:

    CostScaling() {}

  public:

    /// \brief Constructor.

    ///

    /// The constructor of the class.

    ///

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

    CostScaling(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 CostScaling must be signed");

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

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

    CostScaling& 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).

    ///

      _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. By default, 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()

    CostScaling& reset() {

      // Resize vectors

      _node_num = countNodes(_graph);

      _arc_num = countArcs(_graph);

      _res_node_num = _node_num + 1;

      _res_arc_num = 2 * (_arc_num + _node_num);

      _root = _node_num;

      _first_out.resize(_res_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);

      _scost.resize(_res_arc_num);

      _supply.resize(_res_node_num);

      _res_cap.resize(_res_arc_num);

      _cost.resize(_res_arc_num);

      _pi.resize(_res_node_num);

      _excess.resize(_res_node_num);

      _next_out.resize(_res_node_num);

      // Copy the graph

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

      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)];

        }

        _forward[j] = false;

        _source[j] = i;

        _target[j] = _root;

        _reverse[j] = k;

        _forward[k] = true;

        _source[k] = _root;

        _target[k] = i;

        _reverse[k] = j;

        ++j; ++k;

      }

      _first_out[i] = j;

      _first_out[_res_node_num] = k;

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

        int fi = _arc_idf[a];

        int bi = _arc_idb[a];

        _reverse[fi] = bi;

        _reverse[bi] = fi;

      }

      // Reset parameters

      resetParams();

      return *this;

    }

    /// @}

    /// \name Query Functions

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

          int ra = _reverse[a];

          _res_cap[a] = -_sum_supply + 1;

          _res_cap[ra] = -_excess[u];

          _cost[a] = 0;

          _cost[ra] = 0;

          _excess[u] = 0;

        }

      } else {

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

          Value fa = flow[a];

          _res_cap[_arc_idf[a]] = cap[a] - fa;

          _res_cap[_arc_idb[a]] = fa;

        }

        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {

          int ra = _reverse[a];

          _res_cap[a] = 0;

          _res_cap[ra] = 0;

          _cost[a] = 0;

          _cost[ra] = 0;

        }

      }

      // Initialize data structures for buckets

      _max_rank = _alpha * _res_node_num;

      _buckets.resize(_max_rank);

      _bucket_next.resize(_res_node_num + 1);

      _bucket_prev.resize(_res_node_num + 1);

      _rank.resize(_res_node_num + 1);

      return OPTIMAL;

    }

    // Execute the algorithm and transform the results

    void start(Method method) {

      const int MAX_PARTIAL_PATH_LENGTH = 4;

      switch (method) {

        case PUSH:

          startPush();

          break;

        case AUGMENT:

          startAugment(_res_node_num - 1);

          break;

        case PARTIAL_AUGMENT:

          startAugment(MAX_PARTIAL_PATH_LENGTH);

          break;

      }

        }

      }

      // Handle non-zero lower bounds

      if (_have_lower) {

        int limit = _first_out[_root];

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

          if (!_forward[j]) _res_cap[j] += _lower[j];

        }

      }

    }

    // Initialize a cost scaling phase

    void initPhase() {

      // Saturate arcs not satisfying the optimality condition

      for (int u = 0; u != _res_node_num; ++u) {

        int last_out = _first_out[u+1];

        LargeCost pi_u = _pi[u];

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

          Value delta = _res_cap[a];

          if (delta > 0) {

            int v = _target[a];

            if (_cost[a] + pi_u - _pi[v] < 0) {

              _excess[u] -= delta;

              _excess[v] += delta;

              _res_cap[a] = 0;

              _res_cap[_reverse[a]] += delta;

            }

          }

        }

      }

      // Find active nodes (i.e. nodes with positive excess)

      for (int u = 0; u != _res_node_num; ++u) {

        if (_excess[u] > 0) _active_nodes.push_back(u);

      }

      // Initialize the next arcs

      for (int u = 0; u != _res_node_num; ++u) {

        _next_out[u] = _first_out[u];

      }

    }

    // Price (potential) refinement heuristic

    bool priceRefinement() {

      // Stack for stroing the topological order

      IntVector stack(_res_node_num);

      int stack_top;