3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
19 #ifndef LEMON_COST_SCALING_H
20 #define LEMON_COST_SCALING_H
22 /// \ingroup min_cost_flow
25 /// \brief Cost scaling algorithm for finding a minimum cost flow.
28 #include <lemon/graph_adaptor.h>
29 #include <lemon/graph_utils.h>
30 #include <lemon/maps.h>
31 #include <lemon/math.h>
33 #include <lemon/circulation.h>
34 #include <lemon/bellman_ford.h>
38 /// \addtogroup min_cost_flow
41 /// \brief Implementation of the cost scaling algorithm for finding a
42 /// minimum cost flow.
44 /// \ref CostScaling implements the cost scaling algorithm performing
45 /// generalized push-relabel operations for finding a minimum cost
48 /// \tparam Graph The directed graph type the algorithm runs on.
49 /// \tparam LowerMap The type of the lower bound map.
50 /// \tparam CapacityMap The type of the capacity (upper bound) map.
51 /// \tparam CostMap The type of the cost (length) map.
52 /// \tparam SupplyMap The type of the supply map.
55 /// - Edge capacities and costs should be \e non-negative \e integers.
56 /// - Supply values should be \e signed \e integers.
57 /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
58 /// - \c CapacityMap::Value and \c SupplyMap::Value must be
59 /// convertible to each other.
60 /// - All value types must be convertible to \c CostMap::Value, which
61 /// must be signed type.
63 /// \note Edge costs are multiplied with the number of nodes during
64 /// the algorithm so overflow problems may arise more easily than with
65 /// other minimum cost flow algorithms.
66 /// If it is available, <tt>long long int</tt> type is used instead of
67 /// <tt>long int</tt> in the inside computations.
69 /// \author Peter Kovacs
71 template < typename Graph,
72 typename LowerMap = typename Graph::template EdgeMap<int>,
73 typename CapacityMap = typename Graph::template EdgeMap<int>,
74 typename CostMap = typename Graph::template EdgeMap<int>,
75 typename SupplyMap = typename Graph::template NodeMap<int> >
78 GRAPH_TYPEDEFS(typename Graph);
80 typedef typename CapacityMap::Value Capacity;
81 typedef typename CostMap::Value Cost;
82 typedef typename SupplyMap::Value Supply;
83 typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
84 typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
86 typedef ResGraphAdaptor< const Graph, Capacity,
87 CapacityEdgeMap, CapacityEdgeMap > ResGraph;
88 typedef typename ResGraph::Edge ResEdge;
90 #if defined __GNUC__ && !defined __STRICT_ANSI__
91 typedef long long int LCost;
93 typedef long int LCost;
95 typedef typename Graph::template EdgeMap<LCost> LargeCostMap;
99 /// The type of the flow map.
100 typedef CapacityEdgeMap FlowMap;
101 /// The type of the potential map.
102 typedef typename Graph::template NodeMap<LCost> PotentialMap;
106 /// \brief Map adaptor class for handling residual edge costs.
108 /// \ref ResidualCostMap is a map adaptor class for handling
109 /// residual edge costs.
110 class ResidualCostMap : public MapBase<ResEdge, LCost>
114 const LargeCostMap &_cost_map;
119 ResidualCostMap(const LargeCostMap &cost_map) :
120 _cost_map(cost_map) {}
123 LCost operator[](const ResEdge &e) const {
124 return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
127 }; //class ResidualCostMap
129 /// \brief Map adaptor class for handling reduced edge costs.
131 /// \ref ReducedCostMap is a map adaptor class for handling reduced
133 class ReducedCostMap : public MapBase<Edge, LCost>
138 const LargeCostMap &_cost_map;
139 const PotentialMap &_pot_map;
144 ReducedCostMap( const Graph &gr,
145 const LargeCostMap &cost_map,
146 const PotentialMap &pot_map ) :
147 _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
150 LCost operator[](const Edge &e) const {
151 return _cost_map[e] + _pot_map[_gr.source(e)]
152 - _pot_map[_gr.target(e)];
155 }; //class ReducedCostMap
160 static const int ALPHA = 4;
162 // Paramters for heuristics
163 static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
164 static const int BF_HEURISTIC_BOUND_FACTOR = 3;
168 // The directed graph the algorithm runs on
170 // The original lower bound map
171 const LowerMap *_lower;
172 // The modified capacity map
173 CapacityEdgeMap _capacity;
174 // The original cost map
175 const CostMap &_orig_cost;
176 // The scaled cost map
178 // The modified supply map
179 SupplyNodeMap _supply;
182 // Edge map of the current flow
184 // Node map of the current potentials
185 PotentialMap _potential;
187 // The residual graph
189 // The residual cost map
190 ResidualCostMap _res_cost;
191 // The reduced cost map
192 ReducedCostMap _red_cost;
194 SupplyNodeMap _excess;
195 // The epsilon parameter used for cost scaling
200 /// \brief General constructor of the class (with lower bounds).
202 /// General constructor of the class (with lower bounds).
204 /// \param graph The directed graph the algorithm runs on.
205 /// \param lower The lower bounds of the edges.
206 /// \param capacity The capacities (upper bounds) of the edges.
207 /// \param cost The cost (length) values of the edges.
208 /// \param supply The supply values of the nodes (signed).
209 CostScaling( const Graph &graph,
210 const LowerMap &lower,
211 const CapacityMap &capacity,
213 const SupplyMap &supply ) :
214 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
215 _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
216 _res_graph(graph, _capacity, _flow), _res_cost(_cost),
217 _red_cost(graph, _cost, _potential), _excess(graph, 0)
219 // Removing non-zero lower bounds
220 _capacity = subMap(capacity, lower);
222 for (NodeIt n(_graph); n != INVALID; ++n) {
223 Supply s = supply[n];
224 for (InEdgeIt e(_graph, n); e != INVALID; ++e)
226 for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
231 _valid_supply = sum == 0;
234 /// \brief General constructor of the class (without lower bounds).
236 /// General constructor of the class (without lower bounds).
238 /// \param graph The directed graph the algorithm runs on.
239 /// \param capacity The capacities (upper bounds) of the edges.
240 /// \param cost The cost (length) values of the edges.
241 /// \param supply The supply values of the nodes (signed).
242 CostScaling( const Graph &graph,
243 const CapacityMap &capacity,
245 const SupplyMap &supply ) :
246 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
247 _cost(graph), _supply(supply), _flow(graph, 0), _potential(graph, 0),
248 _res_graph(graph, _capacity, _flow), _res_cost(_cost),
249 _red_cost(graph, _cost, _potential), _excess(graph, 0)
251 // Checking the sum of supply values
253 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
254 _valid_supply = sum == 0;
257 /// \brief Simple constructor of the class (with lower bounds).
259 /// Simple constructor of the class (with lower bounds).
261 /// \param graph The directed graph the algorithm runs on.
262 /// \param lower The lower bounds of the edges.
263 /// \param capacity The capacities (upper bounds) of the edges.
264 /// \param cost The cost (length) values of the edges.
265 /// \param s The source node.
266 /// \param t The target node.
267 /// \param flow_value The required amount of flow from node \c s
268 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
269 CostScaling( const Graph &graph,
270 const LowerMap &lower,
271 const CapacityMap &capacity,
274 Supply flow_value ) :
275 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
276 _cost(graph), _supply(graph), _flow(graph, 0), _potential(graph, 0),
277 _res_graph(graph, _capacity, _flow), _res_cost(_cost),
278 _red_cost(graph, _cost, _potential), _excess(graph, 0)
280 // Removing nonzero lower bounds
281 _capacity = subMap(capacity, lower);
282 for (NodeIt n(_graph); n != INVALID; ++n) {
284 if (n == s) sum = flow_value;
285 if (n == t) sum = -flow_value;
286 for (InEdgeIt e(_graph, n); e != INVALID; ++e)
288 for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
292 _valid_supply = true;
295 /// \brief Simple constructor of the class (without lower bounds).
297 /// Simple constructor of the class (without lower bounds).
299 /// \param graph The directed graph the algorithm runs on.
300 /// \param capacity The capacities (upper bounds) of the edges.
301 /// \param cost The cost (length) values of the edges.
302 /// \param s The source node.
303 /// \param t The target node.
304 /// \param flow_value The required amount of flow from node \c s
305 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
306 CostScaling( const Graph &graph,
307 const CapacityMap &capacity,
310 Supply flow_value ) :
311 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
312 _cost(graph), _supply(graph, 0), _flow(graph, 0), _potential(graph, 0),
313 _res_graph(graph, _capacity, _flow), _res_cost(_cost),
314 _red_cost(graph, _cost, _potential), _excess(graph, 0)
316 _supply[s] = flow_value;
317 _supply[t] = -flow_value;
318 _valid_supply = true;
321 /// \brief Runs the algorithm.
323 /// Runs the algorithm.
325 /// \return \c true if a feasible flow can be found.
330 /// \brief Returns a const reference to the edge map storing the
333 /// Returns a const reference to the edge map storing the found flow.
335 /// \pre \ref run() must be called before using this function.
336 const FlowMap& flowMap() const {
340 /// \brief Returns a const reference to the node map storing the
341 /// found potentials (the dual solution).
343 /// Returns a const reference to the node map storing the found
344 /// potentials (the dual solution).
346 /// \pre \ref run() must be called before using this function.
347 const PotentialMap& potentialMap() const {
351 /// \brief Returns the total cost of the found flow.
353 /// Returns the total cost of the found flow. The complexity of the
354 /// function is \f$ O(e) \f$.
356 /// \pre \ref run() must be called before using this function.
357 Cost totalCost() const {
359 for (EdgeIt e(_graph); e != INVALID; ++e)
360 c += _flow[e] * _orig_cost[e];
366 /// Initializes the algorithm.
368 if (!_valid_supply) return false;
370 // Initializing the scaled cost map and the epsilon parameter
372 int node_num = countNodes(_graph);
373 for (EdgeIt e(_graph); e != INVALID; ++e) {
374 _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA;
375 if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
377 _epsilon = max_cost * node_num;
379 // Finding a feasible flow using Circulation
380 Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
382 circulation( _graph, constMap<Edge>((Capacity)0), _capacity,
384 return circulation.flowMap(_flow).run();
388 /// Executes the algorithm.
390 std::deque<Node> active_nodes;
391 typename Graph::template NodeMap<bool> hyper(_graph, false);
393 int node_num = countNodes(_graph);
394 for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ?
395 1 : _epsilon / ALPHA )
397 // Performing price refinement heuristic using Bellman-Ford
399 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
400 typedef ShiftMap<ResidualCostMap> ShiftCostMap;
401 ShiftCostMap shift_cost(_res_cost, _epsilon);
402 BellmanFord<ResGraph, ShiftCostMap> bf(_res_graph, shift_cost);
405 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
406 for (int i = 0; i < K && !done; ++i)
407 done = bf.processNextWeakRound();
409 for (NodeIt n(_graph); n != INVALID; ++n)
410 _potential[n] = bf.dist(n);
415 // Saturating edges not satisfying the optimality condition
417 for (EdgeIt e(_graph); e != INVALID; ++e) {
418 if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
419 delta = _capacity[e] - _flow[e];
420 _excess[_graph.source(e)] -= delta;
421 _excess[_graph.target(e)] += delta;
422 _flow[e] = _capacity[e];
424 if (_flow[e] > 0 && -_red_cost[e] < 0) {
425 _excess[_graph.target(e)] -= _flow[e];
426 _excess[_graph.source(e)] += _flow[e];
431 // Finding active nodes (i.e. nodes with positive excess)
432 for (NodeIt n(_graph); n != INVALID; ++n)
433 if (_excess[n] > 0) active_nodes.push_back(n);
435 // Performing push and relabel operations
436 while (active_nodes.size() > 0) {
437 Node n = active_nodes[0], t;
438 bool relabel_enabled = true;
440 // Performing push operations if there are admissible edges
441 if (_excess[n] > 0) {
442 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
443 if (_capacity[e] - _flow[e] > 0 && _red_cost[e] < 0) {
444 delta = _capacity[e] - _flow[e] <= _excess[n] ?
445 _capacity[e] - _flow[e] : _excess[n];
446 t = _graph.target(e);
448 // Push-look-ahead heuristic
449 Capacity ahead = -_excess[t];
450 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
451 if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
452 ahead += _capacity[oe] - _flow[oe];
454 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
455 if (_flow[ie] > 0 && -_red_cost[ie] < 0)
458 if (ahead < 0) ahead = 0;
460 // Pushing flow along the edge
465 active_nodes.push_front(t);
467 relabel_enabled = false;
473 if (_excess[t] > 0 && _excess[t] <= delta)
474 active_nodes.push_back(t);
477 if (_excess[n] == 0) break;
482 if (_excess[n] > 0) {
483 for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
484 if (_flow[e] > 0 && -_red_cost[e] < 0) {
485 delta = _flow[e] <= _excess[n] ? _flow[e] : _excess[n];
486 t = _graph.source(e);
488 // Push-look-ahead heuristic
489 Capacity ahead = -_excess[t];
490 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
491 if (_capacity[oe] - _flow[oe] > 0 && _red_cost[oe] < 0)
492 ahead += _capacity[oe] - _flow[oe];
494 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
495 if (_flow[ie] > 0 && -_red_cost[ie] < 0)
498 if (ahead < 0) ahead = 0;
500 // Pushing flow along the edge
505 active_nodes.push_front(t);
507 relabel_enabled = false;
513 if (_excess[t] > 0 && _excess[t] <= delta)
514 active_nodes.push_back(t);
517 if (_excess[n] == 0) break;
522 if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
523 // Performing relabel operation if the node is still active
524 LCost min_red_cost = std::numeric_limits<LCost>::max();
525 for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
526 if ( _capacity[oe] - _flow[oe] > 0 &&
527 _red_cost[oe] < min_red_cost )
528 min_red_cost = _red_cost[oe];
530 for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) {
531 if (_flow[ie] > 0 && -_red_cost[ie] < min_red_cost)
532 min_red_cost = -_red_cost[ie];
534 _potential[n] -= min_red_cost + _epsilon;
538 // Removing active nodes with non-positive excess
539 while ( active_nodes.size() > 0 &&
540 _excess[active_nodes[0]] <= 0 &&
541 !hyper[active_nodes[0]] ) {
542 active_nodes.pop_front();
547 // Handling non-zero lower bounds
549 for (EdgeIt e(_graph); e != INVALID; ++e)
550 _flow[e] += (*_lower)[e];
555 }; //class CostScaling
561 #endif //LEMON_COST_SCALING_H