# HG changeset patch # User kpeter # Date 1224072251 0 # Node ID c51b320bc51c9a5f1add3c6137d67f6a70cf3be4 # Parent dc4dd5fc0e2526a23ab23cade6d8ccbd131459c2 Major improvement in the cost scaling algorithm - Add a new variant that use the partial augment-relabel method. - Use this method instead of push-relabel by default. - Use the "Early Termination" heuristic instead of "Price Refinement". Using the new method and heuristic the algorithm proved to be 2-2.5 times faster on all input files. diff -r dc4dd5fc0e25 -r c51b320bc51c lemon/cost_scaling.h --- a/lemon/cost_scaling.h Wed Oct 08 09:17:01 2008 +0000 +++ b/lemon/cost_scaling.h Wed Oct 15 12:04:11 2008 +0000 @@ -20,7 +20,6 @@ #define LEMON_COST_SCALING_H /// \ingroup min_cost_flow -/// /// \file /// \brief Cost scaling algorithm for finding a minimum cost flow. @@ -42,7 +41,7 @@ /// minimum cost flow. /// /// \ref CostScaling implements the cost scaling algorithm performing - /// generalized push-relabel operations for finding a minimum cost + /// augment/push and relabel operations for finding a minimum cost /// flow. /// /// \tparam Graph The directed graph type the algorithm runs on. @@ -116,8 +115,8 @@ _cost_map(cost_map) {} ///\e - typename Map::Value operator[](const ResEdge &e) const { - return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e]; + inline typename Map::Value operator[](const ResEdge &e) const { + return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e]; } }; //class ResidualCostMap @@ -142,7 +141,7 @@ _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {} ///\e - LCost operator[](const Edge &e) const { + inline LCost operator[](const Edge &e) const { return _cost_map[e] + _pot_map[_gr.source(e)] - _pot_map[_gr.target(e)]; } @@ -151,15 +150,6 @@ private: - // Scaling factor - static const int ALPHA = 4; - - // Paramters for heuristics - static const int BF_HEURISTIC_EPSILON_BOUND = 5000; - static const int BF_HEURISTIC_BOUND_FACTOR = 3; - - private: - // The directed graph the algorithm runs on const Graph &_graph; // The original lower bound map @@ -191,6 +181,8 @@ SupplyNodeMap _excess; // The epsilon parameter used for cost scaling LCost _epsilon; + // The scaling factor + int _alpha; public: @@ -213,7 +205,7 @@ _potential(NULL), _local_potential(false), _res_cost(_cost), _res_graph(NULL), _red_cost(NULL), _excess(graph, 0) { - // Removing non-zero lower bounds + // Remove non-zero lower bounds _capacity = subMap(capacity, lower); Supply sum = 0; for (NodeIt n(_graph); n != INVALID; ++n) { @@ -245,7 +237,7 @@ _potential(NULL), _local_potential(false), _res_cost(_cost), _res_graph(NULL), _red_cost(NULL), _excess(graph, 0) { - // Checking the sum of supply values + // Check the sum of supply values Supply sum = 0; for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; _valid_supply = sum == 0; @@ -274,7 +266,7 @@ _potential(NULL), _local_potential(false), _res_cost(_cost), _res_graph(NULL), _red_cost(NULL), _excess(graph, 0) { - // Removing nonzero lower bounds + // Remove nonzero lower bounds _capacity = subMap(capacity, lower); for (NodeIt n(_graph); n != INVALID; ++n) { Supply sum = 0; @@ -359,9 +351,18 @@ /// /// Run the algorithm. /// + /// \param partial_augment By default the algorithm performs + /// partial augment and relabel operations in the cost scaling + /// phases. Set this parameter to \c false for using local push and + /// relabel operations instead. + /// /// \return \c true if a feasible flow can be found. - bool run() { - return init() && start(); + bool run(bool partial_augment = true) { + if (partial_augment) { + return init() && startPartialAugment(); + } else { + return init() && startPushRelabel(); + } } /// @} @@ -433,8 +434,10 @@ /// Initialize the algorithm. bool init() { if (!_valid_supply) return false; + // The scaling factor + _alpha = 8; - // Initializing flow and potential maps + // Initialize flow and potential maps if (!_flow) { _flow = new FlowMap(_graph); _local_flow = true; @@ -447,16 +450,16 @@ _red_cost = new ReducedCostMap(_graph, _cost, *_potential); _res_graph = new ResGraph(_graph, _capacity, *_flow); - // Initializing the scaled cost map and the epsilon parameter + // Initialize the scaled cost map and the epsilon parameter Cost max_cost = 0; int node_num = countNodes(_graph); for (EdgeIt e(_graph); e != INVALID; ++e) { - _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA; + _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha; if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e]; } _epsilon = max_cost * node_num; - // Finding a feasible flow using Circulation + // Find a feasible flow using Circulation Circulation< Graph, ConstMap, CapacityEdgeMap, SupplyMap > circulation( _graph, constMap(Capacity(0)), _capacity, @@ -464,35 +467,43 @@ return circulation.flowMap(*_flow).run(); } + /// Execute the algorithm performing partial augmentation and + /// relabel operations. + bool startPartialAugment() { + // Paramters for heuristics + const int BF_HEURISTIC_EPSILON_BOUND = 1000; + const int BF_HEURISTIC_BOUND_FACTOR = 3; + // Maximum augment path length + const int MAX_PATH_LENGTH = 4; - /// Execute the algorithm. - bool start() { + // Variables + typename Graph::template NodeMap pred_edge(_graph); + typename Graph::template NodeMap forward(_graph); + typename Graph::template NodeMap next_out(_graph); + typename Graph::template NodeMap next_in(_graph); + typename Graph::template NodeMap next_dir(_graph); std::deque active_nodes; - typename Graph::template NodeMap hyper(_graph, false); + std::vector path_nodes; int node_num = countNodes(_graph); - for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ? - 1 : _epsilon / ALPHA ) + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? + 1 : _epsilon / _alpha ) { - // Performing price refinement heuristic using Bellman-Ford - // algorithm + // "Early Termination" heuristic: use Bellman-Ford algorithm + // to check if the current flow is optimal if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { typedef ShiftMap< ResidualCostMap > ShiftCostMap; - ShiftCostMap shift_cost(_res_cost, _epsilon); + ShiftCostMap shift_cost(_res_cost, 1); BellmanFord bf(*_res_graph, shift_cost); bf.init(0); bool done = false; int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num)); for (int i = 0; i < K && !done; ++i) done = bf.processNextWeakRound(); - if (done) { - for (NodeIt n(_graph); n != INVALID; ++n) - (*_potential)[n] = bf.dist(n); - continue; - } + if (done) break; } - // Saturating edges not satisfying the optimality condition + // Saturate edges not satisfying the optimality condition Capacity delta; for (EdgeIt e(_graph); e != INVALID; ++e) { if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { @@ -508,21 +519,193 @@ } } - // Finding active nodes (i.e. nodes with positive excess) - for (NodeIt n(_graph); n != INVALID; ++n) + // Find active nodes (i.e. nodes with positive excess) + for (NodeIt n(_graph); n != INVALID; ++n) { if (_excess[n] > 0) active_nodes.push_back(n); + } - // Performing push and relabel operations + // Initialize the next edge maps + for (NodeIt n(_graph); n != INVALID; ++n) { + next_out[n] = OutEdgeIt(_graph, n); + next_in[n] = InEdgeIt(_graph, n); + next_dir[n] = true; + } + + // Perform partial augment and relabel operations while (active_nodes.size() > 0) { + // Select an active node (FIFO selection) + if (_excess[active_nodes[0]] <= 0) { + active_nodes.pop_front(); + continue; + } + Node start = active_nodes[0]; + path_nodes.clear(); + path_nodes.push_back(start); + + // Find an augmenting path from the start node + Node u, tip = start; + LCost min_red_cost; + while ( _excess[tip] >= 0 && + int(path_nodes.size()) <= MAX_PATH_LENGTH ) + { + if (next_dir[tip]) { + for (OutEdgeIt e = next_out[tip]; e != INVALID; ++e) { + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { + u = _graph.target(e); + pred_edge[u] = e; + forward[u] = true; + next_out[tip] = e; + tip = u; + path_nodes.push_back(tip); + goto next_step; + } + } + next_dir[tip] = false; + } + for (InEdgeIt e = next_in[tip]; e != INVALID; ++e) { + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { + u = _graph.source(e); + pred_edge[u] = e; + forward[u] = false; + next_in[tip] = e; + tip = u; + path_nodes.push_back(tip); + goto next_step; + } + } + + // Relabel tip node + min_red_cost = std::numeric_limits::max() / 2; + for (OutEdgeIt oe(_graph, tip); oe != INVALID; ++oe) { + if ( _capacity[oe] - (*_flow)[oe] > 0 && + (*_red_cost)[oe] < min_red_cost ) + min_red_cost = (*_red_cost)[oe]; + } + for (InEdgeIt ie(_graph, tip); ie != INVALID; ++ie) { + if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost) + min_red_cost = -(*_red_cost)[ie]; + } + (*_potential)[tip] -= min_red_cost + _epsilon; + + // Reset the next edge maps + next_out[tip] = OutEdgeIt(_graph, tip); + next_in[tip] = InEdgeIt(_graph, tip); + next_dir[tip] = true; + + // Step back + if (tip != start) { + path_nodes.pop_back(); + tip = path_nodes[path_nodes.size()-1]; + } + + next_step: + continue; + } + + // Augment along the found path (as much flow as possible) + Capacity delta; + for (int i = 1; i < int(path_nodes.size()); ++i) { + u = path_nodes[i]; + delta = forward[u] ? + _capacity[pred_edge[u]] - (*_flow)[pred_edge[u]] : + (*_flow)[pred_edge[u]]; + delta = std::min(delta, _excess[path_nodes[i-1]]); + (*_flow)[pred_edge[u]] += forward[u] ? delta : -delta; + _excess[path_nodes[i-1]] -= delta; + _excess[u] += delta; + if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u); + } + } + } + + // Compute node potentials for the original costs + ResidualCostMap res_cost(_orig_cost); + BellmanFord< ResGraph, ResidualCostMap > + bf(*_res_graph, res_cost); + bf.init(0); bf.start(); + for (NodeIt n(_graph); n != INVALID; ++n) + (*_potential)[n] = bf.dist(n); + + // Handle non-zero lower bounds + if (_lower) { + for (EdgeIt e(_graph); e != INVALID; ++e) + (*_flow)[e] += (*_lower)[e]; + } + return true; + } + + /// Execute the algorithm performing push and relabel operations. + bool startPushRelabel() { + // Paramters for heuristics + const int BF_HEURISTIC_EPSILON_BOUND = 1000; + const int BF_HEURISTIC_BOUND_FACTOR = 3; + + typename Graph::template NodeMap hyper(_graph, false); + typename Graph::template NodeMap pred_edge(_graph); + typename Graph::template NodeMap forward(_graph); + typename Graph::template NodeMap next_out(_graph); + typename Graph::template NodeMap next_in(_graph); + typename Graph::template NodeMap next_dir(_graph); + std::deque active_nodes; + + int node_num = countNodes(_graph); + for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? + 1 : _epsilon / _alpha ) + { + // "Early Termination" heuristic: use Bellman-Ford algorithm + // to check if the current flow is optimal + if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { + typedef ShiftMap< ResidualCostMap > ShiftCostMap; + ShiftCostMap shift_cost(_res_cost, 1); + BellmanFord bf(*_res_graph, shift_cost); + bf.init(0); + bool done = false; + int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num)); + for (int i = 0; i < K && !done; ++i) + done = bf.processNextWeakRound(); + if (done) break; + } + + // Saturate edges not satisfying the optimality condition + Capacity delta; + for (EdgeIt e(_graph); e != INVALID; ++e) { + if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { + delta = _capacity[e] - (*_flow)[e]; + _excess[_graph.source(e)] -= delta; + _excess[_graph.target(e)] += delta; + (*_flow)[e] = _capacity[e]; + } + if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { + _excess[_graph.target(e)] -= (*_flow)[e]; + _excess[_graph.source(e)] += (*_flow)[e]; + (*_flow)[e] = 0; + } + } + + // Find active nodes (i.e. nodes with positive excess) + for (NodeIt n(_graph); n != INVALID; ++n) { + if (_excess[n] > 0) active_nodes.push_back(n); + } + + // Initialize the next edge maps + for (NodeIt n(_graph); n != INVALID; ++n) { + next_out[n] = OutEdgeIt(_graph, n); + next_in[n] = InEdgeIt(_graph, n); + next_dir[n] = true; + } + + // Perform push and relabel operations + while (active_nodes.size() > 0) { + // Select an active node (FIFO selection) Node n = active_nodes[0], t; bool relabel_enabled = true; - // Performing push operations if there are admissible edges - if (_excess[n] > 0) { - for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { + // Perform push operations if there are admissible edges + if (_excess[n] > 0 && next_dir[n]) { + OutEdgeIt e = next_out[n]; + for ( ; e != INVALID; ++e) { if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) { - delta = _capacity[e] - (*_flow)[e] <= _excess[n] ? - _capacity[e] - (*_flow)[e] : _excess[n]; + delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]); t = _graph.target(e); // Push-look-ahead heuristic @@ -537,7 +720,7 @@ } if (ahead < 0) ahead = 0; - // Pushing flow along the edge + // Push flow along the edge if (ahead < delta) { (*_flow)[e] += ahead; _excess[n] -= ahead; @@ -557,12 +740,18 @@ if (_excess[n] == 0) break; } } + if (e != INVALID) { + next_out[n] = e; + } else { + next_dir[n] = false; + } } - if (_excess[n] > 0) { - for (InEdgeIt e(_graph, n); e != INVALID; ++e) { + if (_excess[n] > 0 && !next_dir[n]) { + InEdgeIt e = next_in[n]; + for ( ; e != INVALID; ++e) { if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) { - delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n]; + delta = std::min((*_flow)[e], _excess[n]); t = _graph.source(e); // Push-look-ahead heuristic @@ -577,7 +766,7 @@ } if (ahead < 0) ahead = 0; - // Pushing flow along the edge + // Push flow along the edge if (ahead < delta) { (*_flow)[e] -= ahead; _excess[n] -= ahead; @@ -597,11 +786,12 @@ if (_excess[n] == 0) break; } } + next_in[n] = e; } + // Relabel the node if it is still active (or hyper) if (relabel_enabled && (_excess[n] > 0 || hyper[n])) { - // Performing relabel operation if the node is still active - LCost min_red_cost = std::numeric_limits::max(); + LCost min_red_cost = std::numeric_limits::max() / 2; for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) { if ( _capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < min_red_cost ) @@ -613,9 +803,14 @@ } (*_potential)[n] -= min_red_cost + _epsilon; hyper[n] = false; + + // Reset the next edge maps + next_out[n] = OutEdgeIt(_graph, n); + next_in[n] = InEdgeIt(_graph, n); + next_dir[n] = true; } - // Removing active nodes with non-positive excess + // Remove nodes that are not active nor hyper while ( active_nodes.size() > 0 && _excess[active_nodes[0]] <= 0 && !hyper[active_nodes[0]] ) { @@ -624,7 +819,7 @@ } } - // Computing node potentials for the original costs + // Compute node potentials for the original costs ResidualCostMap res_cost(_orig_cost); BellmanFord< ResGraph, ResidualCostMap > bf(*_res_graph, res_cost); @@ -632,7 +827,7 @@ for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = bf.dist(n); - // Handling non-zero lower bounds + // Handle non-zero lower bounds if (_lower) { for (EdgeIt e(_graph); e != INVALID; ++e) (*_flow)[e] += (*_lower)[e];