1.1 --- a/lemon/cost_scaling.h Wed Oct 08 09:17:01 2008 +0000
1.2 +++ b/lemon/cost_scaling.h Wed Oct 15 12:04:11 2008 +0000
1.3 @@ -20,7 +20,6 @@
1.4 #define LEMON_COST_SCALING_H
1.5
1.6 /// \ingroup min_cost_flow
1.7 -///
1.8 /// \file
1.9 /// \brief Cost scaling algorithm for finding a minimum cost flow.
1.10
1.11 @@ -42,7 +41,7 @@
1.12 /// minimum cost flow.
1.13 ///
1.14 /// \ref CostScaling implements the cost scaling algorithm performing
1.15 - /// generalized push-relabel operations for finding a minimum cost
1.16 + /// augment/push and relabel operations for finding a minimum cost
1.17 /// flow.
1.18 ///
1.19 /// \tparam Graph The directed graph type the algorithm runs on.
1.20 @@ -116,8 +115,8 @@
1.21 _cost_map(cost_map) {}
1.22
1.23 ///\e
1.24 - typename Map::Value operator[](const ResEdge &e) const {
1.25 - return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
1.26 + inline typename Map::Value operator[](const ResEdge &e) const {
1.27 + return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
1.28 }
1.29
1.30 }; //class ResidualCostMap
1.31 @@ -142,7 +141,7 @@
1.32 _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
1.33
1.34 ///\e
1.35 - LCost operator[](const Edge &e) const {
1.36 + inline LCost operator[](const Edge &e) const {
1.37 return _cost_map[e] + _pot_map[_gr.source(e)]
1.38 - _pot_map[_gr.target(e)];
1.39 }
1.40 @@ -151,15 +150,6 @@
1.41
1.42 private:
1.43
1.44 - // Scaling factor
1.45 - static const int ALPHA = 4;
1.46 -
1.47 - // Paramters for heuristics
1.48 - static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
1.49 - static const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.50 -
1.51 - private:
1.52 -
1.53 // The directed graph the algorithm runs on
1.54 const Graph &_graph;
1.55 // The original lower bound map
1.56 @@ -191,6 +181,8 @@
1.57 SupplyNodeMap _excess;
1.58 // The epsilon parameter used for cost scaling
1.59 LCost _epsilon;
1.60 + // The scaling factor
1.61 + int _alpha;
1.62
1.63 public:
1.64
1.65 @@ -213,7 +205,7 @@
1.66 _potential(NULL), _local_potential(false), _res_cost(_cost),
1.67 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
1.68 {
1.69 - // Removing non-zero lower bounds
1.70 + // Remove non-zero lower bounds
1.71 _capacity = subMap(capacity, lower);
1.72 Supply sum = 0;
1.73 for (NodeIt n(_graph); n != INVALID; ++n) {
1.74 @@ -245,7 +237,7 @@
1.75 _potential(NULL), _local_potential(false), _res_cost(_cost),
1.76 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
1.77 {
1.78 - // Checking the sum of supply values
1.79 + // Check the sum of supply values
1.80 Supply sum = 0;
1.81 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
1.82 _valid_supply = sum == 0;
1.83 @@ -274,7 +266,7 @@
1.84 _potential(NULL), _local_potential(false), _res_cost(_cost),
1.85 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
1.86 {
1.87 - // Removing nonzero lower bounds
1.88 + // Remove nonzero lower bounds
1.89 _capacity = subMap(capacity, lower);
1.90 for (NodeIt n(_graph); n != INVALID; ++n) {
1.91 Supply sum = 0;
1.92 @@ -359,9 +351,18 @@
1.93 ///
1.94 /// Run the algorithm.
1.95 ///
1.96 + /// \param partial_augment By default the algorithm performs
1.97 + /// partial augment and relabel operations in the cost scaling
1.98 + /// phases. Set this parameter to \c false for using local push and
1.99 + /// relabel operations instead.
1.100 + ///
1.101 /// \return \c true if a feasible flow can be found.
1.102 - bool run() {
1.103 - return init() && start();
1.104 + bool run(bool partial_augment = true) {
1.105 + if (partial_augment) {
1.106 + return init() && startPartialAugment();
1.107 + } else {
1.108 + return init() && startPushRelabel();
1.109 + }
1.110 }
1.111
1.112 /// @}
1.113 @@ -433,8 +434,10 @@
1.114 /// Initialize the algorithm.
1.115 bool init() {
1.116 if (!_valid_supply) return false;
1.117 + // The scaling factor
1.118 + _alpha = 8;
1.119
1.120 - // Initializing flow and potential maps
1.121 + // Initialize flow and potential maps
1.122 if (!_flow) {
1.123 _flow = new FlowMap(_graph);
1.124 _local_flow = true;
1.125 @@ -447,16 +450,16 @@
1.126 _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
1.127 _res_graph = new ResGraph(_graph, _capacity, *_flow);
1.128
1.129 - // Initializing the scaled cost map and the epsilon parameter
1.130 + // Initialize the scaled cost map and the epsilon parameter
1.131 Cost max_cost = 0;
1.132 int node_num = countNodes(_graph);
1.133 for (EdgeIt e(_graph); e != INVALID; ++e) {
1.134 - _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA;
1.135 + _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
1.136 if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
1.137 }
1.138 _epsilon = max_cost * node_num;
1.139
1.140 - // Finding a feasible flow using Circulation
1.141 + // Find a feasible flow using Circulation
1.142 Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
1.143 SupplyMap >
1.144 circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
1.145 @@ -464,35 +467,43 @@
1.146 return circulation.flowMap(*_flow).run();
1.147 }
1.148
1.149 + /// Execute the algorithm performing partial augmentation and
1.150 + /// relabel operations.
1.151 + bool startPartialAugment() {
1.152 + // Paramters for heuristics
1.153 + const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.154 + const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.155 + // Maximum augment path length
1.156 + const int MAX_PATH_LENGTH = 4;
1.157
1.158 - /// Execute the algorithm.
1.159 - bool start() {
1.160 + // Variables
1.161 + typename Graph::template NodeMap<Edge> pred_edge(_graph);
1.162 + typename Graph::template NodeMap<bool> forward(_graph);
1.163 + typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
1.164 + typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
1.165 + typename Graph::template NodeMap<bool> next_dir(_graph);
1.166 std::deque<Node> active_nodes;
1.167 - typename Graph::template NodeMap<bool> hyper(_graph, false);
1.168 + std::vector<Node> path_nodes;
1.169
1.170 int node_num = countNodes(_graph);
1.171 - for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ?
1.172 - 1 : _epsilon / ALPHA )
1.173 + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.174 + 1 : _epsilon / _alpha )
1.175 {
1.176 - // Performing price refinement heuristic using Bellman-Ford
1.177 - // algorithm
1.178 + // "Early Termination" heuristic: use Bellman-Ford algorithm
1.179 + // to check if the current flow is optimal
1.180 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.181 typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.182 - ShiftCostMap shift_cost(_res_cost, _epsilon);
1.183 + ShiftCostMap shift_cost(_res_cost, 1);
1.184 BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.185 bf.init(0);
1.186 bool done = false;
1.187 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.188 for (int i = 0; i < K && !done; ++i)
1.189 done = bf.processNextWeakRound();
1.190 - if (done) {
1.191 - for (NodeIt n(_graph); n != INVALID; ++n)
1.192 - (*_potential)[n] = bf.dist(n);
1.193 - continue;
1.194 - }
1.195 + if (done) break;
1.196 }
1.197
1.198 - // Saturating edges not satisfying the optimality condition
1.199 + // Saturate edges not satisfying the optimality condition
1.200 Capacity delta;
1.201 for (EdgeIt e(_graph); e != INVALID; ++e) {
1.202 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.203 @@ -508,21 +519,193 @@
1.204 }
1.205 }
1.206
1.207 - // Finding active nodes (i.e. nodes with positive excess)
1.208 - for (NodeIt n(_graph); n != INVALID; ++n)
1.209 + // Find active nodes (i.e. nodes with positive excess)
1.210 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.211 if (_excess[n] > 0) active_nodes.push_back(n);
1.212 + }
1.213
1.214 - // Performing push and relabel operations
1.215 + // Initialize the next edge maps
1.216 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.217 + next_out[n] = OutEdgeIt(_graph, n);
1.218 + next_in[n] = InEdgeIt(_graph, n);
1.219 + next_dir[n] = true;
1.220 + }
1.221 +
1.222 + // Perform partial augment and relabel operations
1.223 while (active_nodes.size() > 0) {
1.224 + // Select an active node (FIFO selection)
1.225 + if (_excess[active_nodes[0]] <= 0) {
1.226 + active_nodes.pop_front();
1.227 + continue;
1.228 + }
1.229 + Node start = active_nodes[0];
1.230 + path_nodes.clear();
1.231 + path_nodes.push_back(start);
1.232 +
1.233 + // Find an augmenting path from the start node
1.234 + Node u, tip = start;
1.235 + LCost min_red_cost;
1.236 + while ( _excess[tip] >= 0 &&
1.237 + int(path_nodes.size()) <= MAX_PATH_LENGTH )
1.238 + {
1.239 + if (next_dir[tip]) {
1.240 + for (OutEdgeIt e = next_out[tip]; e != INVALID; ++e) {
1.241 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.242 + u = _graph.target(e);
1.243 + pred_edge[u] = e;
1.244 + forward[u] = true;
1.245 + next_out[tip] = e;
1.246 + tip = u;
1.247 + path_nodes.push_back(tip);
1.248 + goto next_step;
1.249 + }
1.250 + }
1.251 + next_dir[tip] = false;
1.252 + }
1.253 + for (InEdgeIt e = next_in[tip]; e != INVALID; ++e) {
1.254 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.255 + u = _graph.source(e);
1.256 + pred_edge[u] = e;
1.257 + forward[u] = false;
1.258 + next_in[tip] = e;
1.259 + tip = u;
1.260 + path_nodes.push_back(tip);
1.261 + goto next_step;
1.262 + }
1.263 + }
1.264 +
1.265 + // Relabel tip node
1.266 + min_red_cost = std::numeric_limits<LCost>::max() / 2;
1.267 + for (OutEdgeIt oe(_graph, tip); oe != INVALID; ++oe) {
1.268 + if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.269 + (*_red_cost)[oe] < min_red_cost )
1.270 + min_red_cost = (*_red_cost)[oe];
1.271 + }
1.272 + for (InEdgeIt ie(_graph, tip); ie != INVALID; ++ie) {
1.273 + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
1.274 + min_red_cost = -(*_red_cost)[ie];
1.275 + }
1.276 + (*_potential)[tip] -= min_red_cost + _epsilon;
1.277 +
1.278 + // Reset the next edge maps
1.279 + next_out[tip] = OutEdgeIt(_graph, tip);
1.280 + next_in[tip] = InEdgeIt(_graph, tip);
1.281 + next_dir[tip] = true;
1.282 +
1.283 + // Step back
1.284 + if (tip != start) {
1.285 + path_nodes.pop_back();
1.286 + tip = path_nodes[path_nodes.size()-1];
1.287 + }
1.288 +
1.289 + next_step:
1.290 + continue;
1.291 + }
1.292 +
1.293 + // Augment along the found path (as much flow as possible)
1.294 + Capacity delta;
1.295 + for (int i = 1; i < int(path_nodes.size()); ++i) {
1.296 + u = path_nodes[i];
1.297 + delta = forward[u] ?
1.298 + _capacity[pred_edge[u]] - (*_flow)[pred_edge[u]] :
1.299 + (*_flow)[pred_edge[u]];
1.300 + delta = std::min(delta, _excess[path_nodes[i-1]]);
1.301 + (*_flow)[pred_edge[u]] += forward[u] ? delta : -delta;
1.302 + _excess[path_nodes[i-1]] -= delta;
1.303 + _excess[u] += delta;
1.304 + if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
1.305 + }
1.306 + }
1.307 + }
1.308 +
1.309 + // Compute node potentials for the original costs
1.310 + ResidualCostMap<CostMap> res_cost(_orig_cost);
1.311 + BellmanFord< ResGraph, ResidualCostMap<CostMap> >
1.312 + bf(*_res_graph, res_cost);
1.313 + bf.init(0); bf.start();
1.314 + for (NodeIt n(_graph); n != INVALID; ++n)
1.315 + (*_potential)[n] = bf.dist(n);
1.316 +
1.317 + // Handle non-zero lower bounds
1.318 + if (_lower) {
1.319 + for (EdgeIt e(_graph); e != INVALID; ++e)
1.320 + (*_flow)[e] += (*_lower)[e];
1.321 + }
1.322 + return true;
1.323 + }
1.324 +
1.325 + /// Execute the algorithm performing push and relabel operations.
1.326 + bool startPushRelabel() {
1.327 + // Paramters for heuristics
1.328 + const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1.329 + const int BF_HEURISTIC_BOUND_FACTOR = 3;
1.330 +
1.331 + typename Graph::template NodeMap<bool> hyper(_graph, false);
1.332 + typename Graph::template NodeMap<Edge> pred_edge(_graph);
1.333 + typename Graph::template NodeMap<bool> forward(_graph);
1.334 + typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
1.335 + typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
1.336 + typename Graph::template NodeMap<bool> next_dir(_graph);
1.337 + std::deque<Node> active_nodes;
1.338 +
1.339 + int node_num = countNodes(_graph);
1.340 + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1.341 + 1 : _epsilon / _alpha )
1.342 + {
1.343 + // "Early Termination" heuristic: use Bellman-Ford algorithm
1.344 + // to check if the current flow is optimal
1.345 + if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1.346 + typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
1.347 + ShiftCostMap shift_cost(_res_cost, 1);
1.348 + BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
1.349 + bf.init(0);
1.350 + bool done = false;
1.351 + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
1.352 + for (int i = 0; i < K && !done; ++i)
1.353 + done = bf.processNextWeakRound();
1.354 + if (done) break;
1.355 + }
1.356 +
1.357 + // Saturate edges not satisfying the optimality condition
1.358 + Capacity delta;
1.359 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.360 + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.361 + delta = _capacity[e] - (*_flow)[e];
1.362 + _excess[_graph.source(e)] -= delta;
1.363 + _excess[_graph.target(e)] += delta;
1.364 + (*_flow)[e] = _capacity[e];
1.365 + }
1.366 + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.367 + _excess[_graph.target(e)] -= (*_flow)[e];
1.368 + _excess[_graph.source(e)] += (*_flow)[e];
1.369 + (*_flow)[e] = 0;
1.370 + }
1.371 + }
1.372 +
1.373 + // Find active nodes (i.e. nodes with positive excess)
1.374 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.375 + if (_excess[n] > 0) active_nodes.push_back(n);
1.376 + }
1.377 +
1.378 + // Initialize the next edge maps
1.379 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.380 + next_out[n] = OutEdgeIt(_graph, n);
1.381 + next_in[n] = InEdgeIt(_graph, n);
1.382 + next_dir[n] = true;
1.383 + }
1.384 +
1.385 + // Perform push and relabel operations
1.386 + while (active_nodes.size() > 0) {
1.387 + // Select an active node (FIFO selection)
1.388 Node n = active_nodes[0], t;
1.389 bool relabel_enabled = true;
1.390
1.391 - // Performing push operations if there are admissible edges
1.392 - if (_excess[n] > 0) {
1.393 - for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
1.394 + // Perform push operations if there are admissible edges
1.395 + if (_excess[n] > 0 && next_dir[n]) {
1.396 + OutEdgeIt e = next_out[n];
1.397 + for ( ; e != INVALID; ++e) {
1.398 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
1.399 - delta = _capacity[e] - (*_flow)[e] <= _excess[n] ?
1.400 - _capacity[e] - (*_flow)[e] : _excess[n];
1.401 + delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
1.402 t = _graph.target(e);
1.403
1.404 // Push-look-ahead heuristic
1.405 @@ -537,7 +720,7 @@
1.406 }
1.407 if (ahead < 0) ahead = 0;
1.408
1.409 - // Pushing flow along the edge
1.410 + // Push flow along the edge
1.411 if (ahead < delta) {
1.412 (*_flow)[e] += ahead;
1.413 _excess[n] -= ahead;
1.414 @@ -557,12 +740,18 @@
1.415 if (_excess[n] == 0) break;
1.416 }
1.417 }
1.418 + if (e != INVALID) {
1.419 + next_out[n] = e;
1.420 + } else {
1.421 + next_dir[n] = false;
1.422 + }
1.423 }
1.424
1.425 - if (_excess[n] > 0) {
1.426 - for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
1.427 + if (_excess[n] > 0 && !next_dir[n]) {
1.428 + InEdgeIt e = next_in[n];
1.429 + for ( ; e != INVALID; ++e) {
1.430 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
1.431 - delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n];
1.432 + delta = std::min((*_flow)[e], _excess[n]);
1.433 t = _graph.source(e);
1.434
1.435 // Push-look-ahead heuristic
1.436 @@ -577,7 +766,7 @@
1.437 }
1.438 if (ahead < 0) ahead = 0;
1.439
1.440 - // Pushing flow along the edge
1.441 + // Push flow along the edge
1.442 if (ahead < delta) {
1.443 (*_flow)[e] -= ahead;
1.444 _excess[n] -= ahead;
1.445 @@ -597,11 +786,12 @@
1.446 if (_excess[n] == 0) break;
1.447 }
1.448 }
1.449 + next_in[n] = e;
1.450 }
1.451
1.452 + // Relabel the node if it is still active (or hyper)
1.453 if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
1.454 - // Performing relabel operation if the node is still active
1.455 - LCost min_red_cost = std::numeric_limits<LCost>::max();
1.456 + LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
1.457 for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
1.458 if ( _capacity[oe] - (*_flow)[oe] > 0 &&
1.459 (*_red_cost)[oe] < min_red_cost )
1.460 @@ -613,9 +803,14 @@
1.461 }
1.462 (*_potential)[n] -= min_red_cost + _epsilon;
1.463 hyper[n] = false;
1.464 +
1.465 + // Reset the next edge maps
1.466 + next_out[n] = OutEdgeIt(_graph, n);
1.467 + next_in[n] = InEdgeIt(_graph, n);
1.468 + next_dir[n] = true;
1.469 }
1.470
1.471 - // Removing active nodes with non-positive excess
1.472 + // Remove nodes that are not active nor hyper
1.473 while ( active_nodes.size() > 0 &&
1.474 _excess[active_nodes[0]] <= 0 &&
1.475 !hyper[active_nodes[0]] ) {
1.476 @@ -624,7 +819,7 @@
1.477 }
1.478 }
1.479
1.480 - // Computing node potentials for the original costs
1.481 + // Compute node potentials for the original costs
1.482 ResidualCostMap<CostMap> res_cost(_orig_cost);
1.483 BellmanFord< ResGraph, ResidualCostMap<CostMap> >
1.484 bf(*_res_graph, res_cost);
1.485 @@ -632,7 +827,7 @@
1.486 for (NodeIt n(_graph); n != INVALID; ++n)
1.487 (*_potential)[n] = bf.dist(n);
1.488
1.489 - // Handling non-zero lower bounds
1.490 + // Handle non-zero lower bounds
1.491 if (_lower) {
1.492 for (EdgeIt e(_graph); e != INVALID; ++e)
1.493 (*_flow)[e] += (*_lower)[e];