lemon-0.x: comparison lemon/cost

equal deleted inserted replaced

-:42fe7654f163
+:64c49446f513
 #ifndef LEMON_COST_SCALING_H
 #define LEMON_COST_SCALING_H
 /// \ingroup min_cost_flow
-///
 /// \file
 /// \brief Cost scaling algorithm for finding a minimum cost flow.
 #include <deque>
 #include <lemon/graph_adaptor.h>
 /// \brief Implementation of the cost scaling algorithm for finding a
 /// 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.
 /// \tparam LowerMap The type of the lower bound map.
 /// \tparam CapacityMap The type of the capacity (upper bound) map.
 ///\e
 ResidualCostMap(const Map &cost_map) :
 _cost_map(cost_map) {}
 ///\e
-typename Map::Value operator[](const ResEdge &e) const {
+inline typename Map::Value operator[](const ResEdge &e) const {
-return ResGraph::forward(e) ?  _cost_map[e] : -_cost_map[e];
+return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
 }
 }; //class ResidualCostMap
 /// \brief Map adaptor class for handling reduced edge costs.
 const LargeCostMap &cost_map,
 const PotentialMap &pot_map ) :
 _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)];
 }
 }; //class ReducedCostMap
-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;
 ReducedCostMap *_red_cost;
 // The excess map
 SupplyNodeMap _excess;
 // The epsilon parameter used for cost scaling
 LCost _epsilon;
+// The scaling factor
+int _alpha;
 public:
 /// \brief General constructor (with lower bounds).
 ///
 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
 _cost(graph), _supply(graph), _flow(NULL), _local_flow(false),
 _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) {
 Supply s = supply[n];
 for (InEdgeIt e(_graph, n); e != INVALID; ++e)
 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
 _cost(graph), _supply(supply), _flow(NULL), _local_flow(false),
 _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;
 }
 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
 _cost(graph), _supply(graph), _flow(NULL), _local_flow(false),
 _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;
 if (n == s) sum =  flow_value;
 if (n == t) sum = -flow_value;
 /// \brief Run the algorithm.
 ///
 /// 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() {
+bool run(bool partial_augment = true) {
-return init() && start();
+if (partial_augment) {
+return init() && startPartialAugment();
+} else {
+return init() && startPushRelabel();
+}
 }
 /// @}
 /// \name Query Functions
 private:
 /// Initialize the algorithm.
 bool init() {
 if (!_valid_supply) return false;
+// The scaling factor
-// Initializing flow and potential maps
+_alpha = 8;
+// Initialize flow and potential maps
 if (!_flow) {
 _flow = new FlowMap(_graph);
 _local_flow = true;
 }
 if (!_potential) {
 }
 _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,
 _supply );
 return circulation.flowMap(*_flow).run();
 }
+/// Execute the algorithm performing partial augmentation and
-/// Execute the algorithm.
+/// relabel operations.
-bool start() {
+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 ?
+for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
-1 : _epsilon / ALPHA )
+1 : _epsilon / _alpha )
 {
-// Performing price refinement heuristic using Bellman-Ford
+// "Early Termination" heuristic: use Bellman-Ford algorithm
-// 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) {
+if (done) break;
-for (NodeIt n(_graph); n != INVALID; ++n)
+}
-(*_potential)[n] = bf.dist(n);
-continue;
+// Saturate edges not satisfying the optimality condition
-}
-}
-// Saturating 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.source(e)] += (*_flow)[e];
 (*_flow)[e] = 0;
 }
 }
-// Finding active nodes (i.e. nodes with positive excess)
+// Find active nodes (i.e. nodes with positive excess)
-for (NodeIt n(_graph); n != INVALID; ++n)
+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<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;
+}
+// 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
+// Perform push operations if there are admissible edges
-if (_excess[n] > 0) {
+if (_excess[n] > 0 && next_dir[n]) {
-for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
+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] ?
+delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
-_capacity[e] - (*_flow)[e] : _excess[n];
 t = _graph.target(e);
 // Push-look-ahead heuristic
 Capacity ahead = -_excess[t];
 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
 ahead += (*_flow)[ie];
 }
 if (ahead < 0) ahead = 0;
-// Pushing flow along the edge
+// Push flow along the edge
 if (ahead < delta) {
 (*_flow)[e] += ahead;
 _excess[n] -= ahead;
 _excess[t] += ahead;
 active_nodes.push_front(t);
 }
 if (_excess[n] == 0) break;
 }
 }
-}
+if (e != INVALID) {
+next_out[n] = e;
-if (_excess[n] > 0) {
+} else {
-for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
+next_dir[n] = false;
+}
+}
+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
 Capacity ahead = -_excess[t];
 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
 ahead += (*_flow)[ie];
 }
 if (ahead < 0) ahead = 0;
-// Pushing flow along the edge
+// Push flow along the edge
 if (ahead < delta) {
 (*_flow)[e] -= ahead;
 _excess[n] -= ahead;
 _excess[t] += ahead;
 active_nodes.push_front(t);
 }
 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() / 2;
-LCost min_red_cost = std::numeric_limits<LCost>::max();
 for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
 if ( _capacity[oe] - (*_flow)[oe] > 0 &&
 (*_red_cost)[oe] < min_red_cost )
 min_red_cost = (*_red_cost)[oe];
 }
 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
 min_red_cost = -(*_red_cost)[ie];
 }
 (*_potential)[n] -= min_red_cost + _epsilon;
 hyper[n] = false;
-}
+// Reset the next edge maps
-// Removing active nodes with non-positive excess
+next_out[n] = OutEdgeIt(_graph, n);
+next_in[n] = InEdgeIt(_graph, n);
+next_dir[n] = true;
+}
+// Remove nodes that are not active nor hyper
 while ( active_nodes.size() > 0 &&
 _excess[active_nodes[0]] <= 0 &&
 !hyper[active_nodes[0]] ) {
 active_nodes.pop_front();
 }
 }
 }
-// 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);
 bf.init(0); bf.start();
 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];
 }
 return true;

changeset 2626	324cfbf66a12
parent 2623	90defb96ee61
child 2629	84354c78b068