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_algs
25 /// \brief Capacity Scaling algorithm for finding a minimum cost flow.
29 #include <lemon/core.h>
30 #include <lemon/bin_heap.h>
34 /// \brief Default traits class of CapacityScaling algorithm.
36 /// Default traits class of CapacityScaling algorithm.
37 /// \tparam GR Digraph type.
38 /// \tparam V The value type used for flow amounts, capacity bounds
39 /// and supply values. By default it is \c int.
40 /// \tparam C The value type used for costs and potentials.
41 /// By default it is the same as \c V.
42 template <typename GR, typename V = int, typename C = V>
43 struct CapacityScalingDefaultTraits
45 /// The type of the digraph
47 /// The type of the flow amounts, capacity bounds and supply values
49 /// The type of the arc costs
52 /// \brief The type of the heap used for internal Dijkstra computations.
54 /// The type of the heap used for internal Dijkstra computations.
55 /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
56 /// its priority type must be \c Cost and its cross reference type
57 /// must be \ref RangeMap "RangeMap<int>".
58 typedef BinHeap<Cost, RangeMap<int> > Heap;
61 /// \addtogroup min_cost_flow_algs
64 /// \brief Implementation of the Capacity Scaling algorithm for
65 /// finding a \ref min_cost_flow "minimum cost flow".
67 /// \ref CapacityScaling implements the capacity scaling version
68 /// of the successive shortest path algorithm for finding a
69 /// \ref min_cost_flow "minimum cost flow". It is an efficient dual
72 /// Most of the parameters of the problem (except for the digraph)
73 /// can be given using separate functions, and the algorithm can be
74 /// executed using the \ref run() function. If some parameters are not
75 /// specified, then default values will be used.
77 /// \tparam GR The digraph type the algorithm runs on.
78 /// \tparam V The value type used for flow amounts, capacity bounds
79 /// and supply values in the algorithm. By default it is \c int.
80 /// \tparam C The value type used for costs and potentials in the
81 /// algorithm. By default it is the same as \c V.
83 /// \warning Both value types must be signed and all input data must
85 /// \warning This algorithm does not support negative costs for such
86 /// arcs that have infinite upper bound.
88 template <typename GR, typename V, typename C, typename TR>
90 template < typename GR, typename V = int, typename C = V,
91 typename TR = CapacityScalingDefaultTraits<GR, V, C> >
97 /// The type of the digraph
98 typedef typename TR::Digraph Digraph;
99 /// The type of the flow amounts, capacity bounds and supply values
100 typedef typename TR::Value Value;
101 /// The type of the arc costs
102 typedef typename TR::Cost Cost;
104 /// The type of the heap used for internal Dijkstra computations
105 typedef typename TR::Heap Heap;
107 /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
112 /// \brief Problem type constants for the \c run() function.
114 /// Enum type containing the problem type constants that can be
115 /// returned by the \ref run() function of the algorithm.
117 /// The problem has no feasible solution (flow).
119 /// The problem has optimal solution (i.e. it is feasible and
120 /// bounded), and the algorithm has found optimal flow and node
121 /// potentials (primal and dual solutions).
123 /// The digraph contains an arc of negative cost and infinite
124 /// upper bound. It means that the objective function is unbounded
125 /// on that arc, however note that it could actually be bounded
126 /// over the feasible flows, but this algroithm cannot handle
133 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
135 typedef std::vector<int> IntVector;
136 typedef std::vector<char> BoolVector;
137 typedef std::vector<Value> ValueVector;
138 typedef std::vector<Cost> CostVector;
142 // Data related to the underlying digraph
149 // Parameters of the problem
153 // Data structures for storing the digraph
157 IntVector _first_out;
169 ValueVector _res_cap;
172 IntVector _excess_nodes;
173 IntVector _deficit_nodes;
181 /// \brief Constant for infinite upper bounds (capacities).
183 /// Constant for infinite upper bounds (capacities).
184 /// It is \c std::numeric_limits<Value>::infinity() if available,
185 /// \c std::numeric_limits<Value>::max() otherwise.
190 // Special implementation of the Dijkstra algorithm for finding
191 // shortest paths in the residual network of the digraph with
192 // respect to the reduced arc costs and modifying the node
193 // potentials according to the found distance labels.
194 class ResidualDijkstra
200 const IntVector &_first_out;
201 const IntVector &_target;
202 const CostVector &_cost;
203 const ValueVector &_res_cap;
204 const ValueVector &_excess;
208 IntVector _proc_nodes;
213 ResidualDijkstra(CapacityScaling& cs) :
214 _node_num(cs._node_num), _geq(cs._sum_supply < 0),
215 _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
216 _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
217 _pred(cs._pred), _dist(cs._node_num)
220 int run(int s, Value delta = 1) {
221 RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
222 Heap heap(heap_cross_ref);
228 while (!heap.empty() && _excess[heap.top()] > -delta) {
229 int u = heap.top(), v;
230 Cost d = heap.prio() + _pi[u], dn;
231 _dist[u] = heap.prio();
232 _proc_nodes.push_back(u);
235 // Traverse outgoing residual arcs
236 int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
237 for (int a = _first_out[u]; a != last_out; ++a) {
238 if (_res_cap[a] < delta) continue;
240 switch (heap.state(v)) {
242 heap.push(v, d + _cost[a] - _pi[v]);
246 dn = d + _cost[a] - _pi[v];
248 heap.decrease(v, dn);
252 case Heap::POST_HEAP:
257 if (heap.empty()) return -1;
259 // Update potentials of processed nodes
261 Cost dt = heap.prio();
262 for (int i = 0; i < int(_proc_nodes.size()); ++i) {
263 _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
269 }; //class ResidualDijkstra
273 /// \name Named Template Parameters
276 template <typename T>
277 struct SetHeapTraits : public Traits {
281 /// \brief \ref named-templ-param "Named parameter" for setting
284 /// \ref named-templ-param "Named parameter" for setting \c Heap
285 /// type, which is used for internal Dijkstra computations.
286 /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
287 /// its priority type must be \c Cost and its cross reference type
288 /// must be \ref RangeMap "RangeMap<int>".
289 template <typename T>
291 : public CapacityScaling<GR, V, C, SetHeapTraits<T> > {
292 typedef CapacityScaling<GR, V, C, SetHeapTraits<T> > Create;
299 /// \brief Constructor.
301 /// The constructor of the class.
303 /// \param graph The digraph the algorithm runs on.
304 CapacityScaling(const GR& graph) :
305 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
306 INF(std::numeric_limits<Value>::has_infinity ?
307 std::numeric_limits<Value>::infinity() :
308 std::numeric_limits<Value>::max())
310 // Check the value types
311 LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
312 "The flow type of CapacityScaling must be signed");
313 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
314 "The cost type of CapacityScaling must be signed");
317 _node_num = countNodes(_graph);
318 _arc_num = countArcs(_graph);
319 _res_arc_num = 2 * (_arc_num + _node_num);
323 _first_out.resize(_node_num + 1);
324 _forward.resize(_res_arc_num);
325 _source.resize(_res_arc_num);
326 _target.resize(_res_arc_num);
327 _reverse.resize(_res_arc_num);
329 _lower.resize(_res_arc_num);
330 _upper.resize(_res_arc_num);
331 _cost.resize(_res_arc_num);
332 _supply.resize(_node_num);
334 _res_cap.resize(_res_arc_num);
335 _pi.resize(_node_num);
336 _excess.resize(_node_num);
337 _pred.resize(_node_num);
340 int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
341 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
345 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
347 for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
351 _target[j] = _node_id[_graph.runningNode(a)];
353 for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
357 _target[j] = _node_id[_graph.runningNode(a)];
370 _first_out[_node_num] = k;
371 for (ArcIt a(_graph); a != INVALID; ++a) {
372 int fi = _arc_idf[a];
373 int bi = _arc_idb[a];
383 /// The parameters of the algorithm can be specified using these
388 /// \brief Set the lower bounds on the arcs.
390 /// This function sets the lower bounds on the arcs.
391 /// If it is not used before calling \ref run(), the lower bounds
392 /// will be set to zero on all arcs.
394 /// \param map An arc map storing the lower bounds.
395 /// Its \c Value type must be convertible to the \c Value type
396 /// of the algorithm.
398 /// \return <tt>(*this)</tt>
399 template <typename LowerMap>
400 CapacityScaling& lowerMap(const LowerMap& map) {
402 for (ArcIt a(_graph); a != INVALID; ++a) {
403 _lower[_arc_idf[a]] = map[a];
404 _lower[_arc_idb[a]] = map[a];
409 /// \brief Set the upper bounds (capacities) on the arcs.
411 /// This function sets the upper bounds (capacities) on the arcs.
412 /// If it is not used before calling \ref run(), the upper bounds
413 /// will be set to \ref INF on all arcs (i.e. the flow value will be
414 /// unbounded from above on each arc).
416 /// \param map An arc map storing the upper bounds.
417 /// Its \c Value type must be convertible to the \c Value type
418 /// of the algorithm.
420 /// \return <tt>(*this)</tt>
421 template<typename UpperMap>
422 CapacityScaling& upperMap(const UpperMap& map) {
423 for (ArcIt a(_graph); a != INVALID; ++a) {
424 _upper[_arc_idf[a]] = map[a];
429 /// \brief Set the costs of the arcs.
431 /// This function sets the costs of the arcs.
432 /// If it is not used before calling \ref run(), the costs
433 /// will be set to \c 1 on all arcs.
435 /// \param map An arc map storing the costs.
436 /// Its \c Value type must be convertible to the \c Cost type
437 /// of the algorithm.
439 /// \return <tt>(*this)</tt>
440 template<typename CostMap>
441 CapacityScaling& costMap(const CostMap& map) {
442 for (ArcIt a(_graph); a != INVALID; ++a) {
443 _cost[_arc_idf[a]] = map[a];
444 _cost[_arc_idb[a]] = -map[a];
449 /// \brief Set the supply values of the nodes.
451 /// This function sets the supply values of the nodes.
452 /// If neither this function nor \ref stSupply() is used before
453 /// calling \ref run(), the supply of each node will be set to zero.
455 /// \param map A node map storing the supply values.
456 /// Its \c Value type must be convertible to the \c Value type
457 /// of the algorithm.
459 /// \return <tt>(*this)</tt>
460 template<typename SupplyMap>
461 CapacityScaling& supplyMap(const SupplyMap& map) {
462 for (NodeIt n(_graph); n != INVALID; ++n) {
463 _supply[_node_id[n]] = map[n];
468 /// \brief Set single source and target nodes and a supply value.
470 /// This function sets a single source node and a single target node
471 /// and the required flow value.
472 /// If neither this function nor \ref supplyMap() is used before
473 /// calling \ref run(), the supply of each node will be set to zero.
475 /// Using this function has the same effect as using \ref supplyMap()
476 /// with such a map in which \c k is assigned to \c s, \c -k is
477 /// assigned to \c t and all other nodes have zero supply value.
479 /// \param s The source node.
480 /// \param t The target node.
481 /// \param k The required amount of flow from node \c s to node \c t
482 /// (i.e. the supply of \c s and the demand of \c t).
484 /// \return <tt>(*this)</tt>
485 CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
486 for (int i = 0; i != _node_num; ++i) {
489 _supply[_node_id[s]] = k;
490 _supply[_node_id[t]] = -k;
496 /// \name Execution control
497 /// The algorithm can be executed using \ref run().
501 /// \brief Run the algorithm.
503 /// This function runs the algorithm.
504 /// The paramters can be specified using functions \ref lowerMap(),
505 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
508 /// CapacityScaling<ListDigraph> cs(graph);
509 /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
510 /// .supplyMap(sup).run();
513 /// This function can be called more than once. All the parameters
514 /// that have been given are kept for the next call, unless
515 /// \ref reset() is called, thus only the modified parameters
516 /// have to be set again. See \ref reset() for examples.
517 /// However the underlying digraph must not be modified after this
518 /// class have been constructed, since it copies and extends the graph.
520 /// \param factor The capacity scaling factor. It must be larger than
521 /// one to use scaling. If it is less or equal to one, then scaling
522 /// will be disabled.
524 /// \return \c INFEASIBLE if no feasible flow exists,
525 /// \n \c OPTIMAL if the problem has optimal solution
526 /// (i.e. it is feasible and bounded), and the algorithm has found
527 /// optimal flow and node potentials (primal and dual solutions),
528 /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
529 /// and infinite upper bound. It means that the objective function
530 /// is unbounded on that arc, however note that it could actually be
531 /// bounded over the feasible flows, but this algroithm cannot handle
535 ProblemType run(int factor = 4) {
537 ProblemType pt = init();
538 if (pt != OPTIMAL) return pt;
542 /// \brief Reset all the parameters that have been given before.
544 /// This function resets all the paramaters that have been given
545 /// before using functions \ref lowerMap(), \ref upperMap(),
546 /// \ref costMap(), \ref supplyMap(), \ref stSupply().
548 /// It is useful for multiple run() calls. If this function is not
549 /// used, all the parameters given before are kept for the next
551 /// However, the underlying digraph must not be modified after this
552 /// class have been constructed, since it copies and extends the graph.
556 /// CapacityScaling<ListDigraph> cs(graph);
559 /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
560 /// .supplyMap(sup).run();
562 /// // Run again with modified cost map (reset() is not called,
563 /// // so only the cost map have to be set again)
565 /// cs.costMap(cost).run();
567 /// // Run again from scratch using reset()
568 /// // (the lower bounds will be set to zero on all arcs)
570 /// cs.upperMap(capacity).costMap(cost)
571 /// .supplyMap(sup).run();
574 /// \return <tt>(*this)</tt>
575 CapacityScaling& reset() {
576 for (int i = 0; i != _node_num; ++i) {
579 for (int j = 0; j != _res_arc_num; ++j) {
582 _cost[j] = _forward[j] ? 1 : -1;
590 /// \name Query Functions
591 /// The results of the algorithm can be obtained using these
593 /// The \ref run() function must be called before using them.
597 /// \brief Return the total cost of the found flow.
599 /// This function returns the total cost of the found flow.
600 /// Its complexity is O(e).
602 /// \note The return type of the function can be specified as a
603 /// template parameter. For example,
605 /// cs.totalCost<double>();
607 /// It is useful if the total cost cannot be stored in the \c Cost
608 /// type of the algorithm, which is the default return type of the
611 /// \pre \ref run() must be called before using this function.
612 template <typename Number>
613 Number totalCost() const {
615 for (ArcIt a(_graph); a != INVALID; ++a) {
617 c += static_cast<Number>(_res_cap[i]) *
618 (-static_cast<Number>(_cost[i]));
624 Cost totalCost() const {
625 return totalCost<Cost>();
629 /// \brief Return the flow on the given arc.
631 /// This function returns the flow on the given arc.
633 /// \pre \ref run() must be called before using this function.
634 Value flow(const Arc& a) const {
635 return _res_cap[_arc_idb[a]];
638 /// \brief Return the flow map (the primal solution).
640 /// This function copies the flow value on each arc into the given
641 /// map. The \c Value type of the algorithm must be convertible to
642 /// the \c Value type of the map.
644 /// \pre \ref run() must be called before using this function.
645 template <typename FlowMap>
646 void flowMap(FlowMap &map) const {
647 for (ArcIt a(_graph); a != INVALID; ++a) {
648 map.set(a, _res_cap[_arc_idb[a]]);
652 /// \brief Return the potential (dual value) of the given node.
654 /// This function returns the potential (dual value) of the
657 /// \pre \ref run() must be called before using this function.
658 Cost potential(const Node& n) const {
659 return _pi[_node_id[n]];
662 /// \brief Return the potential map (the dual solution).
664 /// This function copies the potential (dual value) of each node
665 /// into the given map.
666 /// The \c Cost type of the algorithm must be convertible to the
667 /// \c Value type of the map.
669 /// \pre \ref run() must be called before using this function.
670 template <typename PotentialMap>
671 void potentialMap(PotentialMap &map) const {
672 for (NodeIt n(_graph); n != INVALID; ++n) {
673 map.set(n, _pi[_node_id[n]]);
681 // Initialize the algorithm
683 if (_node_num == 0) return INFEASIBLE;
685 // Check the sum of supply values
687 for (int i = 0; i != _root; ++i) {
688 _sum_supply += _supply[i];
690 if (_sum_supply > 0) return INFEASIBLE;
692 // Initialize vectors
693 for (int i = 0; i != _root; ++i) {
695 _excess[i] = _supply[i];
698 // Remove non-zero lower bounds
699 const Value MAX = std::numeric_limits<Value>::max();
702 for (int i = 0; i != _root; ++i) {
703 last_out = _first_out[i+1];
704 for (int j = _first_out[i]; j != last_out; ++j) {
708 _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
710 _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
713 _excess[_target[j]] += c;
720 for (int j = 0; j != _res_arc_num; ++j) {
721 _res_cap[j] = _forward[j] ? _upper[j] : 0;
725 // Handle negative costs
726 for (int i = 0; i != _root; ++i) {
727 last_out = _first_out[i+1] - 1;
728 for (int j = _first_out[i]; j != last_out; ++j) {
729 Value rc = _res_cap[j];
730 if (_cost[j] < 0 && rc > 0) {
731 if (rc >= MAX) return UNBOUNDED;
733 _excess[_target[j]] += rc;
735 _res_cap[_reverse[j]] += rc;
740 // Handle GEQ supply type
741 if (_sum_supply < 0) {
743 _excess[_root] = -_sum_supply;
744 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
745 int ra = _reverse[a];
746 _res_cap[a] = -_sum_supply + 1;
754 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
755 int ra = _reverse[a];
763 // Initialize delta value
766 Value max_sup = 0, max_dem = 0;
767 for (int i = 0; i != _node_num; ++i) {
768 Value ex = _excess[i];
769 if ( ex > max_sup) max_sup = ex;
770 if (-ex > max_dem) max_dem = -ex;
773 for (int j = 0; j != _res_arc_num; ++j) {
774 if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
776 max_sup = std::min(std::min(max_sup, max_dem), max_cap);
777 for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
786 ProblemType start() {
787 // Execute the algorithm
790 pt = startWithScaling();
792 pt = startWithoutScaling();
794 // Handle non-zero lower bounds
796 int limit = _first_out[_root];
797 for (int j = 0; j != limit; ++j) {
798 if (!_forward[j]) _res_cap[j] += _lower[j];
802 // Shift potentials if necessary
803 Cost pr = _pi[_root];
804 if (_sum_supply < 0 || pr > 0) {
805 for (int i = 0; i != _node_num; ++i) {
813 // Execute the capacity scaling algorithm
814 ProblemType startWithScaling() {
815 // Perform capacity scaling phases
817 ResidualDijkstra _dijkstra(*this);
819 // Saturate all arcs not satisfying the optimality condition
821 for (int u = 0; u != _node_num; ++u) {
822 last_out = _sum_supply < 0 ?
823 _first_out[u+1] : _first_out[u+1] - 1;
824 for (int a = _first_out[u]; a != last_out; ++a) {
826 Cost c = _cost[a] + _pi[u] - _pi[v];
827 Value rc = _res_cap[a];
828 if (c < 0 && rc >= _delta) {
832 _res_cap[_reverse[a]] += rc;
837 // Find excess nodes and deficit nodes
838 _excess_nodes.clear();
839 _deficit_nodes.clear();
840 for (int u = 0; u != _node_num; ++u) {
841 Value ex = _excess[u];
842 if (ex >= _delta) _excess_nodes.push_back(u);
843 if (ex <= -_delta) _deficit_nodes.push_back(u);
845 int next_node = 0, next_def_node = 0;
847 // Find augmenting shortest paths
848 while (next_node < int(_excess_nodes.size())) {
849 // Check deficit nodes
851 bool delta_deficit = false;
852 for ( ; next_def_node < int(_deficit_nodes.size());
854 if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
855 delta_deficit = true;
859 if (!delta_deficit) break;
862 // Run Dijkstra in the residual network
863 s = _excess_nodes[next_node];
864 if ((t = _dijkstra.run(s, _delta)) == -1) {
872 // Augment along a shortest path from s to t
873 Value d = std::min(_excess[s], -_excess[t]);
877 while ((a = _pred[u]) != -1) {
878 if (_res_cap[a] < d) d = _res_cap[a];
883 while ((a = _pred[u]) != -1) {
885 _res_cap[_reverse[a]] += d;
891 if (_excess[s] < _delta) ++next_node;
894 if (_delta == 1) break;
895 _delta = _delta <= _factor ? 1 : _delta / _factor;
901 // Execute the successive shortest path algorithm
902 ProblemType startWithoutScaling() {
904 _excess_nodes.clear();
905 for (int i = 0; i != _node_num; ++i) {
906 if (_excess[i] > 0) _excess_nodes.push_back(i);
908 if (_excess_nodes.size() == 0) return OPTIMAL;
911 // Find shortest paths
913 ResidualDijkstra _dijkstra(*this);
914 while ( _excess[_excess_nodes[next_node]] > 0 ||
915 ++next_node < int(_excess_nodes.size()) )
917 // Run Dijkstra in the residual network
918 s = _excess_nodes[next_node];
919 if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
921 // Augment along a shortest path from s to t
922 Value d = std::min(_excess[s], -_excess[t]);
926 while ((a = _pred[u]) != -1) {
927 if (_res_cap[a] < d) d = _res_cap[a];
932 while ((a = _pred[u]) != -1) {
934 _res_cap[_reverse[a]] += d;
944 }; //class CapacityScaling
950 #endif //LEMON_CAPACITY_SCALING_H