lemon-main: comparison lemon/cost

equal deleted inserted replaced

-:e4a4a854f962
+:c18b7e8e2d31
-/* -*- C++ -*-
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
 *
-* This file is a part of LEMON, a generic C++ optimization library
+* This file is a part of LEMON, a generic C++ optimization library.
 *
-* Copyright (C) 2003-2008
+* 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
 /// finding a \ref min_cost_flow "minimum cost flow".
 ///
 /// \ref CostScaling implements a cost scaling algorithm that performs
 /// push/augment and relabel operations for finding a \ref min_cost_flow
 /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
 /// \ref goldberg97efficient, \ref bunnagel98efficient.
 /// It is a highly efficient primal-dual solution method, which
 /// can be viewed as the generalization of the \ref Preflow
 /// "preflow push-relabel" algorithm for the maximum flow problem.
 ///
 /// Most of the parameters of the problem (except for the digraph)
 /// 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
 };
 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:
 IntVector _buckets;
 IntVector _bucket_next;
 IntVector _bucket_prev;
 IntVector _rank;
 int _max_rank;
 // 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.
 // 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
 }
 _supply[_node_id[s]] =  k;
 _supply[_node_id[t]] = -k;
 return *this;
 }
 /// @}
 /// \name Execution control
 /// The algorithm can be executed using \ref run().
 for (int j = limit; j != _res_arc_num; ++j) {
 _lower[j] = 0;
 _upper[j] = INF;
 _scost[j] = 0;
 _scost[_reverse[j]] = 0;
 }
 _have_lower = false;
 return *this;
 }
 /// \brief Reset all the parameters that have been given before.
 _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);
 int fi = _arc_idf[a];
 int bi = _arc_idb[a];
 _reverse[fi] = bi;
 _reverse[bi] = fi;
 }
 // Reset parameters
 resetParams();
 return *this;
 }
 _sum_supply = 0;
 for (int i = 0; i != _root; ++i) {
 _sum_supply += _supply[i];
 }
 if (_sum_supply > 0) return INFEASIBLE;
 // Initialize vectors
 for (int i = 0; i != _res_node_num; ++i) {
 _pi[i] = 0;
 _excess[i] = _supply[i];
 }
 // Remove infinite upper bounds and check negative arcs
 const Value MAX = std::numeric_limits<Value>::max();
 int last_out;
 if (_have_lower) {
 for (int i = 0; i != _root; ++i) {
 _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;
 // 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);
 // Execute the algorithm
 switch (method) {
 case PUSH:
 startPush();
 break;
 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];
 _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];
 }
 }
 // Early termination heuristic
 bool earlyTermination() {
 const double EARLY_TERM_FACTOR = 3.0;
 // Build a static residual graph
 }
 // Global potential update heuristic
 void globalUpdate() {
 int bucket_end = _root + 1;
 // Initialize buckets
 for (int r = 0; r != _max_rank; ++r) {
 _buckets[r] = bucket_end;
 }
 Value total_excess = 0;
 for ( ; r != _max_rank; ++r) {
 while (_buckets[r] != bucket_end) {
 // Remove the first node from the current bucket
 int u = _buckets[r];
 _buckets[r] = _bucket_next[u];
 // Search the incomming arcs of u
 LargeCost pi_u = _pi[u];
 int last_out = _first_out[u+1];
 for (int a = _first_out[u]; a != last_out; ++a) {
 int ra = _reverse[a];
 // Compute the new rank of v
 LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
 int new_rank_v = old_rank_v;
 if (nrc < LargeCost(_max_rank))
 new_rank_v = r + 1 + int(nrc);
 // Change the rank of v
 if (new_rank_v < old_rank_v) {
 _rank[v] = new_rank_v;
 _next_out[v] = _first_out[v];
 // Remove v from its old bucket
 if (old_rank_v < _max_rank) {
 if (_buckets[old_rank_v] == v) {
 _buckets[old_rank_v] = _bucket_next[v];
 } else {
 _bucket_next[_bucket_prev[v]] = _bucket_next[v];
 _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
 }
 }
 // Insert v to its new bucket
 _bucket_next[v] = _buckets[new_rank_v];
 _bucket_prev[_buckets[new_rank_v]] = v;
 _buckets[new_rank_v] = v;
 }
 if (total_excess <= 0) break;
 }
 }
 if (total_excess <= 0) break;
 }
 // Relabel nodes
 for (int u = 0; u != _res_node_num; ++u) {
 int k = std::min(_rank[u], r);
 if (k > 0) {
 _pi[u] -= _epsilon * k;
 const double GLOBAL_UPDATE_FACTOR = 3.0;
 const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
 (_res_node_num + _sup_node_num * _sup_node_num));
 int next_update_limit = global_update_freq;
 int relabel_cnt = 0;
 // Perform cost scaling phases
 std::vector<int> path;
 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
 1 : _epsilon / _alpha )
 {
 // Early termination heuristic
 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
 if (earlyTermination()) break;
 }
 // Initialize current phase
 initPhase();
 // Perform partial augment and relabel operations
 while (true) {
 // Select an active node (FIFO selection)
 while (_active_nodes.size() > 0 &&
 _excess[_active_nodes.front()] <= 0) {
 const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
 (_res_node_num + _sup_node_num * _sup_node_num));
 int next_update_limit = global_update_freq;
 int relabel_cnt = 0;
 // Perform cost scaling phases
 BoolVector hyper(_res_node_num, false);
 LargeCostVector hyper_cost(_res_node_num);
 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
 1 : _epsilon / _alpha )
 {
 // Early termination heuristic
 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
 if (earlyTermination()) break;
 }
 // Initialize current phase
 initPhase();
 // Perform push and relabel operations
 while (_active_nodes.size() > 0) {
 next_node:
 // Select an active node (FIFO selection)
 n = _active_nodes.front();
 last_out = _first_out[n+1];
 pi_n = _pi[n];
 // 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 &&
 _cost[a] + pi_n - _pi[_target[a]] < 0) {
 // Push-look-ahead heuristic
 Value ahead = -_excess[t];
 int last_out_t = _first_out[t+1];
 LargeCost pi_t = _pi[t];
 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
 if (_res_cap[ta] > 0 &&
 _cost[ta] + pi_t - _pi[_target[ta]] < 0)
 ahead += _res_cap[ta];
 if (ahead >= delta) break;
 }
 if (ahead < 0) ahead = 0;
 _pi[n] -= min_red_cost + _epsilon;
 _next_out[n] = _first_out[n];
 hyper[n] = false;
 ++relabel_cnt;
 }
 // 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();
 }
 // Global update heuristic
 if (relabel_cnt >= next_update_limit) {
 globalUpdate();
 for (int u = 0; u != _res_node_num; ++u)
 hyper[u] = false;

changeset 902	a725503acfe9
parent 863	a93f1a27d831
child 919	e0cef67fe565
child 921	140c953ad5d1
child 931	f112c18bc304