RhodeCode

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

  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;

        _pi[_root] = 0;

        _excess[_root] = -_sum_supply;

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

          int ra = _reverse[a];

          _res_cap[a] = -_sum_supply + 1;

          _res_cap[ra] = 0;

          _cost[a] = 0;

          _cost[ra] = 0;

        }

      } else {

        _pi[_root] = 0;

        _excess[_root] = 0;

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

          int ra = _reverse[a];

          _res_cap[a] = 1;

          _res_cap[ra] = 0;

          _cost[a] = 0;

          _cost[ra] = 0;

        }

      }

      // Initialize delta value

      if (_factor > 1) {

        // With scaling

          Value ex = _excess[i];

          if ( ex > max_sup) max_sup =  ex;

          if (-ex > max_dem) max_dem = -ex;

        }

        max_sup = std::min(std::min(max_sup, max_dem), max_cap);

        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;

      } else {

        // Without scaling

        _delta = 1;

      }

      return OPTIMAL;

    }

    ProblemType start() {

      // Execute the algorithm

      ProblemType pt;

      if (_delta > 1)

        pt = startWithScaling();

      else

        pt = startWithoutScaling();

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

    /// By default, the so called \ref PARTIAL_AUGMENT

    /// "Partial Augment-Relabel" method is used, which proved to be

    /// the most efficient and the most robust on various test inputs.

    /// However, the other methods can be selected using the \ref run()

    /// function with the proper parameter.

    enum Method {

      /// Local push operations are used, i.e. flow is moved only on one

      /// admissible arc at once.

      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;

  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::Node, LargeCost> LargeCostNodeMap;

    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;

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

    // Data for a StaticDigraph structure

    typedef std::pair<int, int> IntPair;

    StaticDigraph _sgr;

    std::vector<IntPair> _arc_vec;

    std::vector<LargeCost> _cost_vec;

    LargeCostArcMap _cost_map;

    LargeCostNodeMap _pi_map;

  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 {

      // Initialize maps for Circulation and remove non-zero lower bounds

      ConstMap<Arc, Value> low(0);

      typedef typename Digraph::template ArcMap<Value> ValueArcMap;

      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;

      ValueArcMap cap(_graph), flow(_graph);

      ValueNodeMap sup(_graph);

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

        sup[n] = _supply[_node_id[n]];

      }

      if (_have_lower) {

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

          int j = _arc_idf[a];

          Value c = _lower[j];

          cap[a] = _upper[j] - c;

          sup[_graph.source(a)] -= c;

          sup[_graph.target(a)] += c;

        }

      } else {

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

          cap[a] = _upper[_arc_idf[a]];

        }

      }

      // Find a feasible flow using Circulation

      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>

        circ(_graph, low, cap, sup);

      if (!circ.flowMap(flow).run()) return INFEASIBLE;

      // Set residual capacities and handle GEQ supply type

      if (_sum_supply < 0) {

        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;

          sup[_graph.source(a)] -= fa;

          sup[_graph.target(a)] += fa;

        }

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

          _excess[_node_id[n]] = sup[n];

        }

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

          int u = _target[a];

          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[ra] = 0;

          _cost[a] = 0;

          _cost[ra] = 0;

        }

      }

      return OPTIMAL;

    }

    // Execute the algorithm and transform the results

    void start(Method method) {

      // Maximum path length for partial augment

      const int MAX_PATH_LENGTH = 4;

      // Execute the algorithm

      switch (method) {

        case PUSH:

          startPush();

          break;

        case AUGMENT:

          startAugment();

          break;

        case PARTIAL_AUGMENT:

          startAugment(MAX_PATH_LENGTH);

          break;

      }

      // Compute node potentials for the original costs

      _arc_vec.clear();

      _cost_vec.clear();

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

        if (_res_cap[j] > 0) {

          _arc_vec.push_back(IntPair(_source[j], _target[j]));

          _cost_vec.push_back(_scost[j]);

        }

      }

      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());

      typename BellmanFord<StaticDigraph, LargeCostArcMap>

        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);

      bf.distMap(_pi_map);

      bf.init(0);

      bf.start();

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

        }

      }

    }

    /// Execute the algorithm performing augment and relabel operations

    void startAugment(int max_length = std::numeric_limits<int>::max()) {

      // Paramters for heuristics

      // Perform cost scaling phases

      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?

                                        1 : _epsilon / _alpha )

      {

        }

        // Perform partial augment and relabel operations

        while (true) {

          // Select an active node (FIFO selection)

          while (_active_nodes.size() > 0 &&

                 _excess[_active_nodes.front()] <= 0) {

            _active_nodes.pop_front();

          }

          if (_active_nodes.size() == 0) break;

          int start = _active_nodes.front();

          // Find an augmenting path from the start node

          int tip = start;

            int u;

            for (int a = _next_out[tip]; a != last_out; ++a) {

                _next_out[tip] = a;

                tip = u;

                goto next_step;

              }

            }

            // Relabel tip node

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

              if (_res_cap[a] > 0 && rc < min_red_cost) {

                min_red_cost = rc;

              }

            }

            _pi[tip] -= min_red_cost + _epsilon;

            _next_out[tip] = _first_out[tip];

            // Step back

            if (tip != start) {

            }

          next_step: ;

          }

          // Augment along the found path (as much flow as possible)

          Value delta;

            u = v;

            delta = std::min(_res_cap[pa], _excess[u]);

            _res_cap[pa] -= delta;

            _res_cap[_reverse[pa]] += delta;

            _excess[u] -= delta;

            _excess[v] += delta;

            if (_excess[v] > 0 && _excess[v] <= delta)

              _active_nodes.push_back(v);

          }

        }

      }

    }

    /// Execute the algorithm performing push and relabel operations

    void startPush() {

      // Paramters for heuristics

      // Perform cost scaling phases

      BoolVector hyper(_res_node_num, false);

      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?

                                        1 : _epsilon / _alpha )

      {

        }

        // Perform push and relabel operations

        while (_active_nodes.size() > 0) {

          Value delta;

          int n, t, a, last_out = _res_arc_num;

          // Select an active node (FIFO selection)

          n = _active_nodes.front();

          // Perform push operations if there are admissible arcs

          if (_excess[n] > 0) {

            for (a = _next_out[n]; a != last_out; ++a) {

              if (_res_cap[a] > 0 &&

                delta = std::min(_res_cap[a], _excess[n]);

                t = _target[a];

                // Push-look-ahead heuristic

                Value ahead = -_excess[t];

                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {

                  if (_res_cap[ta] > 0 && 

                    ahead += _res_cap[ta];

                  if (ahead >= delta) break;

                }

                if (ahead < 0) ahead = 0;

                // Push flow along the arc

                  _res_cap[a] -= ahead;

                  _res_cap[_reverse[a]] += ahead;

                  _excess[n] -= ahead;

                  _excess[t] += ahead;

                  _active_nodes.push_front(t);

                  hyper[t] = true;

                  _next_out[n] = a;

                  goto next_node;

                } else {

                  _res_cap[a] -= delta;

                  _res_cap[_reverse[a]] += delta;

                  _excess[n] -= delta;

                  _excess[t] += delta;

                  if (_excess[t] > 0 && _excess[t] <= delta)

                    _active_nodes.push_back(t);

                }

                if (_excess[n] == 0) {

                  _next_out[n] = a;

                  goto remove_nodes;

                }

              }

            }

            _next_out[n] = a;

          }

          // Relabel the node if it is still active (or hyper)

          if (_excess[n] > 0 || hyper[n]) {

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

              if (_res_cap[a] > 0 && rc < min_red_cost) {

                min_red_cost = rc;

              }

            }

            _pi[n] -= min_red_cost + _epsilon;

            hyper[n] = false;

          }

          // Remove nodes that are not active nor hyper

        remove_nodes:

          while ( _active_nodes.size() > 0 &&

                  _excess[_active_nodes.front()] <= 0 &&

                  !hyper[_active_nodes.front()] ) {

            _active_nodes.pop_front();

          }

        }

      }

    }

  }; //class CostScaling

  ///@}

} //namespace lemon

#endif //LEMON_COST_SCALING_H

    enum Method {

      /// A simple cycle-canceling method, which uses the

      /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration

      /// number for detecting negative cycles in the residual network.

      SIMPLE_CYCLE_CANCELING,

      /// The "Minimum Mean Cycle-Canceling" algorithm, which is a

      /// well-known strongly polynomial method

      /// \ref goldberg89cyclecanceling. It improves along a

      /// \ref min_mean_cycle "minimum mean cycle" in each iteration.

      /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).

      MINIMUM_MEAN_CYCLE_CANCELING,

      /// The "Cancel And Tighten" algorithm, which can be viewed as an

      /// improved version of the previous method

      /// \ref goldberg89cyclecanceling.

      /// It is faster both in theory and in practice, its running time

      /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).

      CANCEL_AND_TIGHTEN

    };

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(GR);

    typedef std::vector<int> IntVector;

    typedef std::vector<double> DoubleVector;

    typedef std::vector<Value> ValueVector;

    typedef std::vector<Cost> CostVector;

  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:

    typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap;

    typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap;

  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;

    // Data structures for storing the digraph

    IntNodeMap _node_id;

    IntArcMap _arc_idf;

    IntArcMap _arc_idb;

    IntVector _first_out;

    IntVector _source;

    IntVector _target;

    IntVector _reverse;

    // Node and arc data

    ValueVector _lower;

    ValueVector _upper;

    CostVector _cost;

    ValueVector _supply;

    ValueVector _res_cap;

    CostVector _pi;

    // Data for a StaticDigraph structure

    typedef std::pair<int, int> IntPair;

    StaticDigraph _sgr;

    std::vector<IntPair> _arc_vec;

    std::vector<Cost> _cost_vec;

    IntVector _id_vec;

    CostArcMap _cost_map;

    CostNodeMap _pi_map;

  public:

          if (d < delta) delta = d;

        }

        // Augment along the cycle

        for (SPathArcIt a(cycle); a != INVALID; ++a) {

          int j = _id_vec[_sgr.id(a)];

          _res_cap[j] -= delta;

          _res_cap[_reverse[j]] += delta;

        }

        // Rebuild the residual network        

        buildResidualNetwork();

      }

    }

    // Execute the "Cancel And Tighten" method

    void startCancelAndTighten() {

      // Constants for the min mean cycle computations

      const double LIMIT_FACTOR = 1.0;

      const int MIN_LIMIT = 5;

      // Contruct auxiliary data vectors

      DoubleVector pi(_res_node_num, 0.0);

      IntVector level(_res_node_num);

      IntVector pred_node(_res_node_num);

      IntVector pred_arc(_res_node_num);

      std::vector<int> stack(_res_node_num);

      std::vector<int> proc_vector(_res_node_num);

      // Initialize epsilon

      double epsilon = 0;

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

        if (_res_cap[a] > 0 && -_cost[a] > epsilon)

          epsilon = -_cost[a];

      }

      // Start phases

      Tolerance<double> tol;

      tol.epsilon(1e-6);

      int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));

      if (limit < MIN_LIMIT) limit = MIN_LIMIT;

      int iter = limit;

      while (epsilon * _res_node_num >= 1) {

        // Find and cancel cycles in the admissible network using DFS

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

          reached[u] = false;

          processed[u] = false;

        }

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

      /// It keeps only the several best eligible arcs from the former

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

    // State constants for arcs

    enum ArcStateEnum {

      STATE_UPPER = -1,

      STATE_TREE  =  0,

      STATE_LOWER =  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;

    // Node and arc data

    ValueVector _lower;

    ValueVector _upper;

    ValueVector _cap;