[Lemon-commits] kpeter: r3511 - lemon/trunk/lemon
Lemon SVN
svn at lemon.cs.elte.hu
Wed Oct 15 14:04:12 CEST 2008
Author: kpeter
Date: Wed Oct 15 14:04:11 2008
New Revision: 3511
Modified:
lemon/trunk/lemon/cost_scaling.h
Log:
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.
Modified: lemon/trunk/lemon/cost_scaling.h
==============================================================================
--- lemon/trunk/lemon/cost_scaling.h (original)
+++ lemon/trunk/lemon/cost_scaling.h Wed Oct 15 14:04:11 2008
@@ -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<Edge, Capacity>, CapacityEdgeMap,
SupplyMap >
circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
@@ -464,35 +467,206 @@
return circulation.flowMap(*_flow).run();
}
-
- /// Execute the algorithm.
- bool start() {
+ /// 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;
+
+ // Variables
+ typename Graph::template NodeMap<Edge> pred_edge(_graph);
+ typename Graph::template NodeMap<bool> forward(_graph);
+ typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
+ typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
+ typename Graph::template NodeMap<bool> next_dir(_graph);
std::deque<Node> active_nodes;
- typename Graph::template NodeMap<bool> hyper(_graph, false);
+ std::vector<Node> 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<LargeCostMap> > ShiftCostMap;
- ShiftCostMap shift_cost(_res_cost, _epsilon);
+ ShiftCostMap shift_cost(_res_cost, 1);
BellmanFord<ResGraph, ShiftCostMap> 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);
+ 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 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<LCost>::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<CostMap> res_cost(_orig_cost);
+ BellmanFord< ResGraph, ResidualCostMap<CostMap> >
+ 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<bool> hyper(_graph, false);
+ typename Graph::template NodeMap<Edge> pred_edge(_graph);
+ typename Graph::template NodeMap<bool> forward(_graph);
+ typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
+ typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
+ typename Graph::template NodeMap<bool> next_dir(_graph);
+ std::deque<Node> 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<LargeCostMap> > ShiftCostMap;
+ ShiftCostMap shift_cost(_res_cost, 1);
+ BellmanFord<ResGraph, ShiftCostMap> 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;
}
- // 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 +682,30 @@
}
}
- // 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);
+ }
+
+ // 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;
+ }
- // Performing push and relabel operations
+ // 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<LCost>::max();
+ LCost min_red_cost = std::numeric_limits<LCost>::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<CostMap> res_cost(_orig_cost);
BellmanFord< ResGraph, ResidualCostMap<CostMap> >
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];
More information about the Lemon-commits
mailing list