1.1 --- a/lemon/cost_scaling.h Wed Mar 17 12:35:52 2010 +0100
1.2 +++ b/lemon/cost_scaling.h Sat Mar 06 14:35:12 2010 +0000
1.3 @@ -1,8 +1,8 @@
1.4 -/* -*- C++ -*-
1.5 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.6 *
1.7 - * This file is a part of LEMON, a generic C++ optimization library
1.8 + * This file is a part of LEMON, a generic C++ optimization library.
1.9 *
1.10 - * Copyright (C) 2003-2008
1.11 + * Copyright (C) 2003-2010
1.12 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.13 * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.14 *
1.15 @@ -92,7 +92,7 @@
1.16 /// \ref CostScaling implements a cost scaling algorithm that performs
1.17 /// push/augment and relabel operations for finding a \ref min_cost_flow
1.18 /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
1.19 - /// \ref goldberg97efficient, \ref bunnagel98efficient.
1.20 + /// \ref goldberg97efficient, \ref bunnagel98efficient.
1.21 /// It is a highly efficient primal-dual solution method, which
1.22 /// can be viewed as the generalization of the \ref Preflow
1.23 /// "preflow push-relabel" algorithm for the maximum flow problem.
1.24 @@ -189,7 +189,7 @@
1.25 /// Augment operations are used, i.e. flow is moved on admissible
1.26 /// paths from a node with excess to a node with deficit.
1.27 AUGMENT,
1.28 - /// Partial augment operations are used, i.e. flow is moved on
1.29 + /// Partial augment operations are used, i.e. flow is moved on
1.30 /// admissible paths started from a node with excess, but the
1.31 /// lengths of these paths are limited. This method can be viewed
1.32 /// as a combined version of the previous two operations.
1.33 @@ -208,15 +208,15 @@
1.34 // Note: vector<char> is used instead of vector<bool> for efficiency reasons
1.35
1.36 private:
1.37 -
1.38 +
1.39 template <typename KT, typename VT>
1.40 class StaticVectorMap {
1.41 public:
1.42 typedef KT Key;
1.43 typedef VT Value;
1.44 -
1.45 +
1.46 StaticVectorMap(std::vector<Value>& v) : _v(v) {}
1.47 -
1.48 +
1.49 const Value& operator[](const Key& key) const {
1.50 return _v[StaticDigraph::id(key)];
1.51 }
1.52 @@ -224,7 +224,7 @@
1.53 Value& operator[](const Key& key) {
1.54 return _v[StaticDigraph::id(key)];
1.55 }
1.56 -
1.57 +
1.58 void set(const Key& key, const Value& val) {
1.59 _v[StaticDigraph::id(key)] = val;
1.60 }
1.61 @@ -283,7 +283,7 @@
1.62 IntVector _bucket_prev;
1.63 IntVector _rank;
1.64 int _max_rank;
1.65 -
1.66 +
1.67 // Data for a StaticDigraph structure
1.68 typedef std::pair<int, int> IntPair;
1.69 StaticDigraph _sgr;
1.70 @@ -291,9 +291,9 @@
1.71 std::vector<LargeCost> _cost_vec;
1.72 LargeCostArcMap _cost_map;
1.73 LargeCostNodeMap _pi_map;
1.74 -
1.75 +
1.76 public:
1.77 -
1.78 +
1.79 /// \brief Constant for infinite upper bounds (capacities).
1.80 ///
1.81 /// Constant for infinite upper bounds (capacities).
1.82 @@ -348,7 +348,7 @@
1.83 "The flow type of CostScaling must be signed");
1.84 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.85 "The cost type of CostScaling must be signed");
1.86 -
1.87 +
1.88 // Reset data structures
1.89 reset();
1.90 }
1.91 @@ -464,7 +464,7 @@
1.92 _supply[_node_id[t]] = -k;
1.93 return *this;
1.94 }
1.95 -
1.96 +
1.97 /// @}
1.98
1.99 /// \name Execution control
1.100 @@ -566,7 +566,7 @@
1.101 _upper[j] = INF;
1.102 _scost[j] = 0;
1.103 _scost[_reverse[j]] = 0;
1.104 - }
1.105 + }
1.106 _have_lower = false;
1.107 return *this;
1.108 }
1.109 @@ -601,7 +601,7 @@
1.110 _upper.resize(_res_arc_num);
1.111 _scost.resize(_res_arc_num);
1.112 _supply.resize(_res_node_num);
1.113 -
1.114 +
1.115 _res_cap.resize(_res_arc_num);
1.116 _cost.resize(_res_arc_num);
1.117 _pi.resize(_res_node_num);
1.118 @@ -649,7 +649,7 @@
1.119 _reverse[fi] = bi;
1.120 _reverse[bi] = fi;
1.121 }
1.122 -
1.123 +
1.124 // Reset parameters
1.125 resetParams();
1.126 return *this;
1.127 @@ -758,14 +758,14 @@
1.128 _sum_supply += _supply[i];
1.129 }
1.130 if (_sum_supply > 0) return INFEASIBLE;
1.131 -
1.132 +
1.133
1.134 // Initialize vectors
1.135 for (int i = 0; i != _res_node_num; ++i) {
1.136 _pi[i] = 0;
1.137 _excess[i] = _supply[i];
1.138 }
1.139 -
1.140 +
1.141 // Remove infinite upper bounds and check negative arcs
1.142 const Value MAX = std::numeric_limits<Value>::max();
1.143 int last_out;
1.144 @@ -885,7 +885,7 @@
1.145 _cost[ra] = 0;
1.146 }
1.147 }
1.148 -
1.149 +
1.150 return OPTIMAL;
1.151 }
1.152
1.153 @@ -894,13 +894,13 @@
1.154 // Maximum path length for partial augment
1.155 const int MAX_PATH_LENGTH = 4;
1.156
1.157 - // Initialize data structures for buckets
1.158 + // Initialize data structures for buckets
1.159 _max_rank = _alpha * _res_node_num;
1.160 _buckets.resize(_max_rank);
1.161 _bucket_next.resize(_res_node_num + 1);
1.162 _bucket_prev.resize(_res_node_num + 1);
1.163 _rank.resize(_res_node_num + 1);
1.164 -
1.165 +
1.166 // Execute the algorithm
1.167 switch (method) {
1.168 case PUSH:
1.169 @@ -939,7 +939,7 @@
1.170 }
1.171 }
1.172 }
1.173 -
1.174 +
1.175 // Initialize a cost scaling phase
1.176 void initPhase() {
1.177 // Saturate arcs not satisfying the optimality condition
1.178 @@ -957,7 +957,7 @@
1.179 }
1.180 }
1.181 }
1.182 -
1.183 +
1.184 // Find active nodes (i.e. nodes with positive excess)
1.185 for (int u = 0; u != _res_node_num; ++u) {
1.186 if (_excess[u] > 0) _active_nodes.push_back(u);
1.187 @@ -968,7 +968,7 @@
1.188 _next_out[u] = _first_out[u];
1.189 }
1.190 }
1.191 -
1.192 +
1.193 // Early termination heuristic
1.194 bool earlyTermination() {
1.195 const double EARLY_TERM_FACTOR = 3.0;
1.196 @@ -998,7 +998,7 @@
1.197 // Global potential update heuristic
1.198 void globalUpdate() {
1.199 int bucket_end = _root + 1;
1.200 -
1.201 +
1.202 // Initialize buckets
1.203 for (int r = 0; r != _max_rank; ++r) {
1.204 _buckets[r] = bucket_end;
1.205 @@ -1024,7 +1024,7 @@
1.206 // Remove the first node from the current bucket
1.207 int u = _buckets[r];
1.208 _buckets[r] = _bucket_next[u];
1.209 -
1.210 +
1.211 // Search the incomming arcs of u
1.212 LargeCost pi_u = _pi[u];
1.213 int last_out = _first_out[u+1];
1.214 @@ -1039,12 +1039,12 @@
1.215 int new_rank_v = old_rank_v;
1.216 if (nrc < LargeCost(_max_rank))
1.217 new_rank_v = r + 1 + int(nrc);
1.218 -
1.219 +
1.220 // Change the rank of v
1.221 if (new_rank_v < old_rank_v) {
1.222 _rank[v] = new_rank_v;
1.223 _next_out[v] = _first_out[v];
1.224 -
1.225 +
1.226 // Remove v from its old bucket
1.227 if (old_rank_v < _max_rank) {
1.228 if (_buckets[old_rank_v] == v) {
1.229 @@ -1054,7 +1054,7 @@
1.230 _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
1.231 }
1.232 }
1.233 -
1.234 +
1.235 // Insert v to its new bucket
1.236 _bucket_next[v] = _buckets[new_rank_v];
1.237 _bucket_prev[_buckets[new_rank_v]] = v;
1.238 @@ -1072,7 +1072,7 @@
1.239 }
1.240 if (total_excess <= 0) break;
1.241 }
1.242 -
1.243 +
1.244 // Relabel nodes
1.245 for (int u = 0; u != _res_node_num; ++u) {
1.246 int k = std::min(_rank[u], r);
1.247 @@ -1092,9 +1092,9 @@
1.248 const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
1.249 (_res_node_num + _sup_node_num * _sup_node_num));
1.250 int next_update_limit = global_update_freq;
1.251 -
1.252 +
1.253 int relabel_cnt = 0;
1.254 -
1.255 +
1.256 // Perform cost scaling phases
1.257 std::vector<int> path;
1.258 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.259 @@ -1104,10 +1104,10 @@
1.260 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.261 if (earlyTermination()) break;
1.262 }
1.263 -
1.264 +
1.265 // Initialize current phase
1.266 initPhase();
1.267 -
1.268 +
1.269 // Perform partial augment and relabel operations
1.270 while (true) {
1.271 // Select an active node (FIFO selection)
1.272 @@ -1196,7 +1196,7 @@
1.273 int next_update_limit = global_update_freq;
1.274
1.275 int relabel_cnt = 0;
1.276 -
1.277 +
1.278 // Perform cost scaling phases
1.279 BoolVector hyper(_res_node_num, false);
1.280 LargeCostVector hyper_cost(_res_node_num);
1.281 @@ -1207,7 +1207,7 @@
1.282 if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
1.283 if (earlyTermination()) break;
1.284 }
1.285 -
1.286 +
1.287 // Initialize current phase
1.288 initPhase();
1.289
1.290 @@ -1222,7 +1222,7 @@
1.291 n = _active_nodes.front();
1.292 last_out = _first_out[n+1];
1.293 pi_n = _pi[n];
1.294 -
1.295 +
1.296 // Perform push operations if there are admissible arcs
1.297 if (_excess[n] > 0) {
1.298 for (a = _next_out[n]; a != last_out; ++a) {
1.299 @@ -1236,7 +1236,7 @@
1.300 int last_out_t = _first_out[t+1];
1.301 LargeCost pi_t = _pi[t];
1.302 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
1.303 - if (_res_cap[ta] > 0 &&
1.304 + if (_res_cap[ta] > 0 &&
1.305 _cost[ta] + pi_t - _pi[_target[ta]] < 0)
1.306 ahead += _res_cap[ta];
1.307 if (ahead >= delta) break;
1.308 @@ -1287,7 +1287,7 @@
1.309 hyper[n] = false;
1.310 ++relabel_cnt;
1.311 }
1.312 -
1.313 +
1.314 // Remove nodes that are not active nor hyper
1.315 remove_nodes:
1.316 while ( _active_nodes.size() > 0 &&
1.317 @@ -1295,7 +1295,7 @@
1.318 !hyper[_active_nodes.front()] ) {
1.319 _active_nodes.pop_front();
1.320 }
1.321 -
1.322 +
1.323 // Global update heuristic
1.324 if (relabel_cnt >= next_update_limit) {
1.325 globalUpdate();