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_CAPACITY_SCALING_H
20 #define LEMON_CAPACITY_SCALING_H
22 /// \ingroup min_cost_flow
25 /// \brief Capacity scaling algorithm for finding a minimum cost flow.
28 #include <lemon/bin_heap.h>
32 /// \addtogroup min_cost_flow
35 /// \brief Implementation of the capacity scaling algorithm for
36 /// finding a minimum cost flow.
38 /// \ref CapacityScaling implements the capacity scaling version
39 /// of the successive shortest path algorithm for finding a minimum
42 /// \tparam Digraph The digraph type the algorithm runs on.
43 /// \tparam LowerMap The type of the lower bound map.
44 /// \tparam CapacityMap The type of the capacity (upper bound) map.
45 /// \tparam CostMap The type of the cost (length) map.
46 /// \tparam SupplyMap The type of the supply map.
49 /// - Arc capacities and costs should be \e non-negative \e integers.
50 /// - Supply values should be \e signed \e integers.
51 /// - The value types of the maps should be convertible to each other.
52 /// - \c CostMap::Value must be signed type.
54 /// \author Peter Kovacs
55 template < typename Digraph,
56 typename LowerMap = typename Digraph::template ArcMap<int>,
57 typename CapacityMap = typename Digraph::template ArcMap<int>,
58 typename CostMap = typename Digraph::template ArcMap<int>,
59 typename SupplyMap = typename Digraph::template NodeMap<int> >
62 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
64 typedef typename CapacityMap::Value Capacity;
65 typedef typename CostMap::Value Cost;
66 typedef typename SupplyMap::Value Supply;
67 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
68 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
69 typedef typename Digraph::template NodeMap<Arc> PredMap;
73 /// The type of the flow map.
74 typedef typename Digraph::template ArcMap<Capacity> FlowMap;
75 /// The type of the potential map.
76 typedef typename Digraph::template NodeMap<Cost> PotentialMap;
80 /// \brief Special implementation of the \ref Dijkstra algorithm
81 /// for finding shortest paths in the residual network.
83 /// \ref ResidualDijkstra is a special implementation of the
84 /// \ref Dijkstra algorithm for finding shortest paths in the
85 /// residual network of the digraph with respect to the reduced arc
86 /// costs and modifying the node potentials according to the
87 /// distance of the nodes.
88 class ResidualDijkstra
90 typedef typename Digraph::template NodeMap<int> HeapCrossRef;
91 typedef BinHeap<Cost, HeapCrossRef> Heap;
95 // The digraph the algorithm runs on
96 const Digraph &_graph;
100 const CapacityArcMap &_res_cap;
101 const CostMap &_cost;
102 const SupplyNodeMap &_excess;
103 PotentialMap &_potential;
109 // The processed (i.e. permanently labeled) nodes
110 std::vector<Node> _proc_nodes;
115 ResidualDijkstra( const Digraph &digraph,
117 const CapacityArcMap &res_cap,
119 const SupplyMap &excess,
120 PotentialMap &potential,
122 _graph(digraph), _flow(flow), _res_cap(res_cap), _cost(cost),
123 _excess(excess), _potential(potential), _dist(digraph),
127 /// Run the algorithm from the given source node.
128 Node run(Node s, Capacity delta = 1) {
129 HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
130 Heap heap(heap_cross_ref);
136 while (!heap.empty() && _excess[heap.top()] > -delta) {
137 Node u = heap.top(), v;
138 Cost d = heap.prio() + _potential[u], nd;
139 _dist[u] = heap.prio();
141 _proc_nodes.push_back(u);
143 // Traversing outgoing arcs
144 for (OutArcIt e(_graph, u); e != INVALID; ++e) {
145 if (_res_cap[e] >= delta) {
146 v = _graph.target(e);
147 switch(heap.state(v)) {
149 heap.push(v, d + _cost[e] - _potential[v]);
153 nd = d + _cost[e] - _potential[v];
155 heap.decrease(v, nd);
159 case Heap::POST_HEAP:
165 // Traversing incoming arcs
166 for (InArcIt e(_graph, u); e != INVALID; ++e) {
167 if (_flow[e] >= delta) {
168 v = _graph.source(e);
169 switch(heap.state(v)) {
171 heap.push(v, d - _cost[e] - _potential[v]);
175 nd = d - _cost[e] - _potential[v];
177 heap.decrease(v, nd);
181 case Heap::POST_HEAP:
187 if (heap.empty()) return INVALID;
189 // Updating potentials of processed nodes
191 Cost t_dist = heap.prio();
192 for (int i = 0; i < int(_proc_nodes.size()); ++i)
193 _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
198 }; //class ResidualDijkstra
202 // The digraph the algorithm runs on
203 const Digraph &_graph;
204 // The original lower bound map
205 const LowerMap *_lower;
206 // The modified capacity map
207 CapacityArcMap _capacity;
208 // The original cost map
209 const CostMap &_cost;
210 // The modified supply map
211 SupplyNodeMap _supply;
214 // Arc map of the current flow
217 // Node map of the current potentials
218 PotentialMap *_potential;
219 bool _local_potential;
221 // The residual capacity map
222 CapacityArcMap _res_cap;
224 SupplyNodeMap _excess;
225 // The excess nodes (i.e. nodes with positive excess)
226 std::vector<Node> _excess_nodes;
227 // The deficit nodes (i.e. nodes with negative excess)
228 std::vector<Node> _deficit_nodes;
230 // The delta parameter used for capacity scaling
232 // The maximum number of phases
237 // Implementation of the Dijkstra algorithm for finding augmenting
238 // shortest paths in the residual network
239 ResidualDijkstra *_dijkstra;
243 /// \brief General constructor (with lower bounds).
245 /// General constructor (with lower bounds).
247 /// \param digraph The digraph the algorithm runs on.
248 /// \param lower The lower bounds of the arcs.
249 /// \param capacity The capacities (upper bounds) of the arcs.
250 /// \param cost The cost (length) values of the arcs.
251 /// \param supply The supply values of the nodes (signed).
252 CapacityScaling( const Digraph &digraph,
253 const LowerMap &lower,
254 const CapacityMap &capacity,
256 const SupplyMap &supply ) :
257 _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
258 _supply(digraph), _flow(NULL), _local_flow(false),
259 _potential(NULL), _local_potential(false),
260 _res_cap(digraph), _excess(digraph), _pred(digraph), _dijkstra(NULL)
263 for (NodeIt n(_graph); n != INVALID; ++n) {
264 _supply[n] = supply[n];
265 _excess[n] = supply[n];
268 _valid_supply = sum == 0;
269 for (ArcIt a(_graph); a != INVALID; ++a) {
270 _capacity[a] = capacity[a];
271 _res_cap[a] = capacity[a];
274 // Remove non-zero lower bounds
275 typename LowerMap::Value lcap;
276 for (ArcIt e(_graph); e != INVALID; ++e) {
277 if ((lcap = lower[e]) != 0) {
278 _capacity[e] -= lcap;
280 _supply[_graph.source(e)] -= lcap;
281 _supply[_graph.target(e)] += lcap;
282 _excess[_graph.source(e)] -= lcap;
283 _excess[_graph.target(e)] += lcap;
288 /// \brief General constructor (without lower bounds).
290 /// General constructor (without lower bounds).
292 /// \param digraph The digraph the algorithm runs on.
293 /// \param capacity The capacities (upper bounds) of the arcs.
294 /// \param cost The cost (length) values of the arcs.
295 /// \param supply The supply values of the nodes (signed).
296 CapacityScaling( const Digraph &digraph,
297 const CapacityMap &capacity,
299 const SupplyMap &supply ) :
300 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
301 _supply(supply), _flow(NULL), _local_flow(false),
302 _potential(NULL), _local_potential(false),
303 _res_cap(capacity), _excess(supply), _pred(digraph), _dijkstra(NULL)
305 // Check the sum of supply values
307 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
308 _valid_supply = sum == 0;
311 /// \brief Simple constructor (with lower bounds).
313 /// Simple constructor (with lower bounds).
315 /// \param digraph The digraph the algorithm runs on.
316 /// \param lower The lower bounds of the arcs.
317 /// \param capacity The capacities (upper bounds) of the arcs.
318 /// \param cost The cost (length) values of the arcs.
319 /// \param s The source node.
320 /// \param t The target node.
321 /// \param flow_value The required amount of flow from node \c s
322 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
323 CapacityScaling( const Digraph &digraph,
324 const LowerMap &lower,
325 const CapacityMap &capacity,
328 Supply flow_value ) :
329 _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
330 _supply(digraph, 0), _flow(NULL), _local_flow(false),
331 _potential(NULL), _local_potential(false),
332 _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
334 // Remove non-zero lower bounds
335 _supply[s] = _excess[s] = flow_value;
336 _supply[t] = _excess[t] = -flow_value;
337 typename LowerMap::Value lcap;
338 for (ArcIt e(_graph); e != INVALID; ++e) {
339 if ((lcap = lower[e]) != 0) {
340 _capacity[e] -= lcap;
342 _supply[_graph.source(e)] -= lcap;
343 _supply[_graph.target(e)] += lcap;
344 _excess[_graph.source(e)] -= lcap;
345 _excess[_graph.target(e)] += lcap;
348 _valid_supply = true;
351 /// \brief Simple constructor (without lower bounds).
353 /// Simple constructor (without lower bounds).
355 /// \param digraph The digraph the algorithm runs on.
356 /// \param capacity The capacities (upper bounds) of the arcs.
357 /// \param cost The cost (length) values of the arcs.
358 /// \param s The source node.
359 /// \param t The target node.
360 /// \param flow_value The required amount of flow from node \c s
361 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
362 CapacityScaling( const Digraph &digraph,
363 const CapacityMap &capacity,
366 Supply flow_value ) :
367 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
368 _supply(digraph, 0), _flow(NULL), _local_flow(false),
369 _potential(NULL), _local_potential(false),
370 _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL)
372 _supply[s] = _excess[s] = flow_value;
373 _supply[t] = _excess[t] = -flow_value;
374 _valid_supply = true;
379 if (_local_flow) delete _flow;
380 if (_local_potential) delete _potential;
384 /// \brief Set the flow map.
386 /// Set the flow map.
388 /// \return \c (*this)
389 CapacityScaling& flowMap(FlowMap &map) {
398 /// \brief Set the potential map.
400 /// Set the potential map.
402 /// \return \c (*this)
403 CapacityScaling& potentialMap(PotentialMap &map) {
404 if (_local_potential) {
406 _local_potential = false;
412 /// \name Execution control
416 /// \brief Run the algorithm.
418 /// This function runs the algorithm.
420 /// \param scaling Enable or disable capacity scaling.
421 /// If the maximum arc capacity and/or the amount of total supply
422 /// is rather small, the algorithm could be slightly faster without
425 /// \return \c true if a feasible flow can be found.
426 bool run(bool scaling = true) {
427 return init(scaling) && start();
432 /// \name Query Functions
433 /// The results of the algorithm can be obtained using these
435 /// \ref lemon::CapacityScaling::run() "run()" must be called before
440 /// \brief Return a const reference to the arc map storing the
443 /// Return a const reference to the arc map storing the found flow.
445 /// \pre \ref run() must be called before using this function.
446 const FlowMap& flowMap() const {
450 /// \brief Return a const reference to the node map storing the
451 /// found potentials (the dual solution).
453 /// Return a const reference to the node map storing the found
454 /// potentials (the dual solution).
456 /// \pre \ref run() must be called before using this function.
457 const PotentialMap& potentialMap() const {
461 /// \brief Return the flow on the given arc.
463 /// Return the flow on the given arc.
465 /// \pre \ref run() must be called before using this function.
466 Capacity flow(const Arc& arc) const {
467 return (*_flow)[arc];
470 /// \brief Return the potential of the given node.
472 /// Return the potential of the given node.
474 /// \pre \ref run() must be called before using this function.
475 Cost potential(const Node& node) const {
476 return (*_potential)[node];
479 /// \brief Return the total cost of the found flow.
481 /// Return the total cost of the found flow. The complexity of the
482 /// function is \f$ O(e) \f$.
484 /// \pre \ref run() must be called before using this function.
485 Cost totalCost() const {
487 for (ArcIt e(_graph); e != INVALID; ++e)
488 c += (*_flow)[e] * _cost[e];
496 /// Initialize the algorithm.
497 bool init(bool scaling) {
498 if (!_valid_supply) return false;
502 _flow = new FlowMap(_graph);
506 _potential = new PotentialMap(_graph);
507 _local_potential = true;
509 for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
510 for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
512 _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost,
513 _excess, *_potential, _pred );
515 // Initializing delta value
518 Supply max_sup = 0, max_dem = 0;
519 for (NodeIt n(_graph); n != INVALID; ++n) {
520 if ( _supply[n] > max_sup) max_sup = _supply[n];
521 if (-_supply[n] > max_dem) max_dem = -_supply[n];
523 Capacity max_cap = 0;
524 for (ArcIt e(_graph); e != INVALID; ++e) {
525 if (_capacity[e] > max_cap) max_cap = _capacity[e];
527 max_sup = std::min(std::min(max_sup, max_dem), max_cap);
529 for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
541 return startWithScaling();
543 return startWithoutScaling();
546 /// Execute the capacity scaling algorithm.
547 bool startWithScaling() {
548 // Processing capacity scaling phases
553 // Saturating all arcs not satisfying the optimality condition
554 for (ArcIt e(_graph); e != INVALID; ++e) {
555 Node u = _graph.source(e), v = _graph.target(e);
556 Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v];
557 if (c < 0 && _res_cap[e] >= _delta) {
558 _excess[u] -= _res_cap[e];
559 _excess[v] += _res_cap[e];
560 (*_flow)[e] = _capacity[e];
563 else if (c > 0 && (*_flow)[e] >= _delta) {
564 _excess[u] += (*_flow)[e];
565 _excess[v] -= (*_flow)[e];
567 _res_cap[e] = _capacity[e];
571 // Finding excess nodes and deficit nodes
572 _excess_nodes.clear();
573 _deficit_nodes.clear();
574 for (NodeIt n(_graph); n != INVALID; ++n) {
575 if (_excess[n] >= _delta) _excess_nodes.push_back(n);
576 if (_excess[n] <= -_delta) _deficit_nodes.push_back(n);
578 int next_node = 0, next_def_node = 0;
580 // Finding augmenting shortest paths
581 while (next_node < int(_excess_nodes.size())) {
582 // Checking deficit nodes
584 bool delta_deficit = false;
585 for ( ; next_def_node < int(_deficit_nodes.size());
587 if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
588 delta_deficit = true;
592 if (!delta_deficit) break;
596 s = _excess_nodes[next_node];
597 if ((t = _dijkstra->run(s, _delta)) == INVALID) {
605 // Augmenting along a shortest path from s to t.
606 Capacity d = std::min(_excess[s], -_excess[t]);
610 while ((e = _pred[u]) != INVALID) {
612 if (u == _graph.target(e)) {
614 u = _graph.source(e);
617 u = _graph.target(e);
623 while ((e = _pred[u]) != INVALID) {
624 if (u == _graph.target(e)) {
627 u = _graph.source(e);
631 u = _graph.target(e);
637 if (_excess[s] < _delta) ++next_node;
640 if (_delta == 1) break;
641 if (++phase_cnt > _phase_num / 4) factor = 2;
642 _delta = _delta <= factor ? 1 : _delta / factor;
645 // Handling non-zero lower bounds
647 for (ArcIt e(_graph); e != INVALID; ++e)
648 (*_flow)[e] += (*_lower)[e];
653 /// Execute the successive shortest path algorithm.
654 bool startWithoutScaling() {
655 // Finding excess nodes
656 for (NodeIt n(_graph); n != INVALID; ++n)
657 if (_excess[n] > 0) _excess_nodes.push_back(n);
658 if (_excess_nodes.size() == 0) return true;
661 // Finding shortest paths
663 while ( _excess[_excess_nodes[next_node]] > 0 ||
664 ++next_node < int(_excess_nodes.size()) )
667 s = _excess_nodes[next_node];
668 if ((t = _dijkstra->run(s)) == INVALID) return false;
670 // Augmenting along a shortest path from s to t
671 Capacity d = std::min(_excess[s], -_excess[t]);
675 while ((e = _pred[u]) != INVALID) {
677 if (u == _graph.target(e)) {
679 u = _graph.source(e);
682 u = _graph.target(e);
688 while ((e = _pred[u]) != INVALID) {
689 if (u == _graph.target(e)) {
692 u = _graph.source(e);
696 u = _graph.target(e);
703 // Handling non-zero lower bounds
705 for (ArcIt e(_graph); e != INVALID; ++e)
706 (*_flow)[e] += (*_lower)[e];
711 }; //class CapacityScaling
717 #endif //LEMON_CAPACITY_SCALING_H