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 /// \addtogroup min_cost_flow_algs
37 /// \brief Implementation of the Capacity Scaling algorithm for
38 /// finding a \ref min_cost_flow "minimum cost flow".
40 /// \ref CapacityScaling implements the capacity scaling version
41 /// of the successive shortest path algorithm for finding a
42 /// \ref min_cost_flow "minimum cost flow". It is an efficient dual
45 /// Most of the parameters of the problem (except for the digraph)
46 /// can be given using separate functions, and the algorithm can be
47 /// executed using the \ref run() function. If some parameters are not
48 /// specified, then default values will be used.
50 /// \tparam GR The digraph type the algorithm runs on.
51 /// \tparam V The value type used for flow amounts, capacity bounds
52 /// and supply values in the algorithm. By default it is \c int.
53 /// \tparam C The value type used for costs and potentials in the
54 /// algorithm. By default it is the same as \c V.
56 /// \warning Both value types must be signed and all input data must
58 /// \warning This algorithm does not support negative costs for such
59 /// arcs that have infinite upper bound.
60 template <typename GR, typename V = int, typename C = V>
65 /// The type of the flow amounts, capacity bounds and supply values
67 /// The type of the arc costs
72 /// \brief Problem type constants for the \c run() function.
74 /// Enum type containing the problem type constants that can be
75 /// returned by the \ref run() function of the algorithm.
77 /// The problem has no feasible solution (flow).
79 /// The problem has optimal solution (i.e. it is feasible and
80 /// bounded), and the algorithm has found optimal flow and node
81 /// potentials (primal and dual solutions).
83 /// The digraph contains an arc of negative cost and infinite
84 /// upper bound. It means that the objective function is unbounded
85 /// on that arc, however note that it could actually be bounded
86 /// over the feasible flows, but this algroithm cannot handle
93 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
95 typedef std::vector<Arc> ArcVector;
96 typedef std::vector<Node> NodeVector;
97 typedef std::vector<int> IntVector;
98 typedef std::vector<bool> BoolVector;
99 typedef std::vector<Value> ValueVector;
100 typedef std::vector<Cost> CostVector;
104 // Data related to the underlying digraph
111 // Parameters of the problem
115 // Data structures for storing the digraph
119 IntVector _first_out;
131 ValueVector _res_cap;
134 IntVector _excess_nodes;
135 IntVector _deficit_nodes;
143 /// \brief Constant for infinite upper bounds (capacities).
145 /// Constant for infinite upper bounds (capacities).
146 /// It is \c std::numeric_limits<Value>::infinity() if available,
147 /// \c std::numeric_limits<Value>::max() otherwise.
152 // Special implementation of the Dijkstra algorithm for finding
153 // shortest paths in the residual network of the digraph with
154 // respect to the reduced arc costs and modifying the node
155 // potentials according to the found distance labels.
156 class ResidualDijkstra
158 typedef RangeMap<int> HeapCrossRef;
159 typedef BinHeap<Cost, HeapCrossRef> Heap;
164 const IntVector &_first_out;
165 const IntVector &_target;
166 const CostVector &_cost;
167 const ValueVector &_res_cap;
168 const ValueVector &_excess;
172 IntVector _proc_nodes;
177 ResidualDijkstra(CapacityScaling& cs) :
178 _node_num(cs._node_num), _first_out(cs._first_out),
179 _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap),
180 _excess(cs._excess), _pi(cs._pi), _pred(cs._pred),
184 int run(int s, Value delta = 1) {
185 HeapCrossRef heap_cross_ref(_node_num, Heap::PRE_HEAP);
186 Heap heap(heap_cross_ref);
192 while (!heap.empty() && _excess[heap.top()] > -delta) {
193 int u = heap.top(), v;
194 Cost d = heap.prio() + _pi[u], dn;
195 _dist[u] = heap.prio();
196 _proc_nodes.push_back(u);
199 // Traverse outgoing residual arcs
200 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
201 if (_res_cap[a] < delta) continue;
203 switch (heap.state(v)) {
205 heap.push(v, d + _cost[a] - _pi[v]);
209 dn = d + _cost[a] - _pi[v];
211 heap.decrease(v, dn);
215 case Heap::POST_HEAP:
220 if (heap.empty()) return -1;
222 // Update potentials of processed nodes
224 Cost dt = heap.prio();
225 for (int i = 0; i < int(_proc_nodes.size()); ++i) {
226 _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
232 }; //class ResidualDijkstra
236 /// \brief Constructor.
238 /// The constructor of the class.
240 /// \param graph The digraph the algorithm runs on.
241 CapacityScaling(const GR& graph) :
242 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
243 INF(std::numeric_limits<Value>::has_infinity ?
244 std::numeric_limits<Value>::infinity() :
245 std::numeric_limits<Value>::max())
247 // Check the value types
248 LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
249 "The flow type of CapacityScaling must be signed");
250 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
251 "The cost type of CapacityScaling must be signed");
254 _node_num = countNodes(_graph);
255 _arc_num = countArcs(_graph);
256 _res_arc_num = 2 * (_arc_num + _node_num);
260 _first_out.resize(_node_num + 1);
261 _forward.resize(_res_arc_num);
262 _source.resize(_res_arc_num);
263 _target.resize(_res_arc_num);
264 _reverse.resize(_res_arc_num);
266 _lower.resize(_res_arc_num);
267 _upper.resize(_res_arc_num);
268 _cost.resize(_res_arc_num);
269 _supply.resize(_node_num);
271 _res_cap.resize(_res_arc_num);
272 _pi.resize(_node_num);
273 _excess.resize(_node_num);
274 _pred.resize(_node_num);
277 int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
278 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
282 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
284 for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
288 _target[j] = _node_id[_graph.runningNode(a)];
290 for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
294 _target[j] = _node_id[_graph.runningNode(a)];
307 _first_out[_node_num] = k;
308 for (ArcIt a(_graph); a != INVALID; ++a) {
309 int fi = _arc_idf[a];
310 int bi = _arc_idb[a];
320 /// The parameters of the algorithm can be specified using these
325 /// \brief Set the lower bounds on the arcs.
327 /// This function sets the lower bounds on the arcs.
328 /// If it is not used before calling \ref run(), the lower bounds
329 /// will be set to zero on all arcs.
331 /// \param map An arc map storing the lower bounds.
332 /// Its \c Value type must be convertible to the \c Value type
333 /// of the algorithm.
335 /// \return <tt>(*this)</tt>
336 template <typename LowerMap>
337 CapacityScaling& lowerMap(const LowerMap& map) {
339 for (ArcIt a(_graph); a != INVALID; ++a) {
340 _lower[_arc_idf[a]] = map[a];
341 _lower[_arc_idb[a]] = map[a];
346 /// \brief Set the upper bounds (capacities) on the arcs.
348 /// This function sets the upper bounds (capacities) on the arcs.
349 /// If it is not used before calling \ref run(), the upper bounds
350 /// will be set to \ref INF on all arcs (i.e. the flow value will be
351 /// unbounded from above on each arc).
353 /// \param map An arc map storing the upper bounds.
354 /// Its \c Value type must be convertible to the \c Value type
355 /// of the algorithm.
357 /// \return <tt>(*this)</tt>
358 template<typename UpperMap>
359 CapacityScaling& upperMap(const UpperMap& map) {
360 for (ArcIt a(_graph); a != INVALID; ++a) {
361 _upper[_arc_idf[a]] = map[a];
366 /// \brief Set the costs of the arcs.
368 /// This function sets the costs of the arcs.
369 /// If it is not used before calling \ref run(), the costs
370 /// will be set to \c 1 on all arcs.
372 /// \param map An arc map storing the costs.
373 /// Its \c Value type must be convertible to the \c Cost type
374 /// of the algorithm.
376 /// \return <tt>(*this)</tt>
377 template<typename CostMap>
378 CapacityScaling& costMap(const CostMap& map) {
379 for (ArcIt a(_graph); a != INVALID; ++a) {
380 _cost[_arc_idf[a]] = map[a];
381 _cost[_arc_idb[a]] = -map[a];
386 /// \brief Set the supply values of the nodes.
388 /// This function sets the supply values of the nodes.
389 /// If neither this function nor \ref stSupply() is used before
390 /// calling \ref run(), the supply of each node will be set to zero.
392 /// \param map A node map storing the supply values.
393 /// Its \c Value type must be convertible to the \c Value type
394 /// of the algorithm.
396 /// \return <tt>(*this)</tt>
397 template<typename SupplyMap>
398 CapacityScaling& supplyMap(const SupplyMap& map) {
399 for (NodeIt n(_graph); n != INVALID; ++n) {
400 _supply[_node_id[n]] = map[n];
405 /// \brief Set single source and target nodes and a supply value.
407 /// This function sets a single source node and a single target node
408 /// and the required flow value.
409 /// If neither this function nor \ref supplyMap() is used before
410 /// calling \ref run(), the supply of each node will be set to zero.
412 /// Using this function has the same effect as using \ref supplyMap()
413 /// with such a map in which \c k is assigned to \c s, \c -k is
414 /// assigned to \c t and all other nodes have zero supply value.
416 /// \param s The source node.
417 /// \param t The target node.
418 /// \param k The required amount of flow from node \c s to node \c t
419 /// (i.e. the supply of \c s and the demand of \c t).
421 /// \return <tt>(*this)</tt>
422 CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
423 for (int i = 0; i != _node_num; ++i) {
426 _supply[_node_id[s]] = k;
427 _supply[_node_id[t]] = -k;
433 /// \name Execution control
437 /// \brief Run the algorithm.
439 /// This function runs the algorithm.
440 /// The paramters can be specified using functions \ref lowerMap(),
441 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
444 /// CapacityScaling<ListDigraph> cs(graph);
445 /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
446 /// .supplyMap(sup).run();
449 /// This function can be called more than once. All the parameters
450 /// that have been given are kept for the next call, unless
451 /// \ref reset() is called, thus only the modified parameters
452 /// have to be set again. See \ref reset() for examples.
453 /// However the underlying digraph must not be modified after this
454 /// class have been constructed, since it copies the digraph.
456 /// \param scaling Enable or disable capacity scaling.
457 /// If the maximum upper bound and/or the amount of total supply
458 /// is rather small, the algorithm could be slightly faster without
461 /// \return \c INFEASIBLE if no feasible flow exists,
462 /// \n \c OPTIMAL if the problem has optimal solution
463 /// (i.e. it is feasible and bounded), and the algorithm has found
464 /// optimal flow and node potentials (primal and dual solutions),
465 /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
466 /// and infinite upper bound. It means that the objective function
467 /// is unbounded on that arc, however note that it could actually be
468 /// bounded over the feasible flows, but this algroithm cannot handle
472 ProblemType run(bool scaling = true) {
473 ProblemType pt = init(scaling);
474 if (pt != OPTIMAL) return pt;
478 /// \brief Reset all the parameters that have been given before.
480 /// This function resets all the paramaters that have been given
481 /// before using functions \ref lowerMap(), \ref upperMap(),
482 /// \ref costMap(), \ref supplyMap(), \ref stSupply().
484 /// It is useful for multiple run() calls. If this function is not
485 /// used, all the parameters given before are kept for the next
487 /// However the underlying digraph must not be modified after this
488 /// class have been constructed, since it copies and extends the graph.
492 /// CapacityScaling<ListDigraph> cs(graph);
495 /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
496 /// .supplyMap(sup).run();
498 /// // Run again with modified cost map (reset() is not called,
499 /// // so only the cost map have to be set again)
501 /// cs.costMap(cost).run();
503 /// // Run again from scratch using reset()
504 /// // (the lower bounds will be set to zero on all arcs)
506 /// cs.upperMap(capacity).costMap(cost)
507 /// .supplyMap(sup).run();
510 /// \return <tt>(*this)</tt>
511 CapacityScaling& reset() {
512 for (int i = 0; i != _node_num; ++i) {
515 for (int j = 0; j != _res_arc_num; ++j) {
518 _cost[j] = _forward[j] ? 1 : -1;
526 /// \name Query Functions
527 /// The results of the algorithm can be obtained using these
529 /// The \ref run() function must be called before using them.
533 /// \brief Return the total cost of the found flow.
535 /// This function returns the total cost of the found flow.
536 /// Its complexity is O(e).
538 /// \note The return type of the function can be specified as a
539 /// template parameter. For example,
541 /// cs.totalCost<double>();
543 /// It is useful if the total cost cannot be stored in the \c Cost
544 /// type of the algorithm, which is the default return type of the
547 /// \pre \ref run() must be called before using this function.
548 template <typename Number>
549 Number totalCost() const {
551 for (ArcIt a(_graph); a != INVALID; ++a) {
553 c += static_cast<Number>(_res_cap[i]) *
554 (-static_cast<Number>(_cost[i]));
560 Cost totalCost() const {
561 return totalCost<Cost>();
565 /// \brief Return the flow on the given arc.
567 /// This function returns the flow on the given arc.
569 /// \pre \ref run() must be called before using this function.
570 Value flow(const Arc& a) const {
571 return _res_cap[_arc_idb[a]];
574 /// \brief Return the flow map (the primal solution).
576 /// This function copies the flow value on each arc into the given
577 /// map. The \c Value type of the algorithm must be convertible to
578 /// the \c Value type of the map.
580 /// \pre \ref run() must be called before using this function.
581 template <typename FlowMap>
582 void flowMap(FlowMap &map) const {
583 for (ArcIt a(_graph); a != INVALID; ++a) {
584 map.set(a, _res_cap[_arc_idb[a]]);
588 /// \brief Return the potential (dual value) of the given node.
590 /// This function returns the potential (dual value) of the
593 /// \pre \ref run() must be called before using this function.
594 Cost potential(const Node& n) const {
595 return _pi[_node_id[n]];
598 /// \brief Return the potential map (the dual solution).
600 /// This function copies the potential (dual value) of each node
601 /// into the given map.
602 /// The \c Cost type of the algorithm must be convertible to the
603 /// \c Value type of the map.
605 /// \pre \ref run() must be called before using this function.
606 template <typename PotentialMap>
607 void potentialMap(PotentialMap &map) const {
608 for (NodeIt n(_graph); n != INVALID; ++n) {
609 map.set(n, _pi[_node_id[n]]);
617 // Initialize the algorithm
618 ProblemType init(bool scaling) {
619 if (_node_num == 0) return INFEASIBLE;
621 // Check the sum of supply values
623 for (int i = 0; i != _root; ++i) {
624 _sum_supply += _supply[i];
626 if (_sum_supply > 0) return INFEASIBLE;
629 for (int i = 0; i != _root; ++i) {
631 _excess[i] = _supply[i];
634 // Remove non-zero lower bounds
636 for (int i = 0; i != _root; ++i) {
637 for (int j = _first_out[i]; j != _first_out[i+1]; ++j) {
641 _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF;
643 _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF;
646 _excess[_target[j]] += c;
653 for (int j = 0; j != _res_arc_num; ++j) {
654 _res_cap[j] = _forward[j] ? _upper[j] : 0;
658 // Handle negative costs
659 for (int u = 0; u != _root; ++u) {
660 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
661 Value rc = _res_cap[a];
662 if (_cost[a] < 0 && rc > 0) {
663 if (rc == INF) return UNBOUNDED;
665 _excess[_target[a]] += rc;
667 _res_cap[_reverse[a]] += rc;
672 // Handle GEQ supply type
673 if (_sum_supply < 0) {
675 _excess[_root] = -_sum_supply;
676 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
678 if (_excess[u] < 0) {
679 _res_cap[a] = -_excess[u] + 1;
683 _res_cap[_reverse[a]] = 0;
685 _cost[_reverse[a]] = 0;
690 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
692 _res_cap[_reverse[a]] = 0;
694 _cost[_reverse[a]] = 0;
698 // Initialize delta value
701 Value max_sup = 0, max_dem = 0;
702 for (int i = 0; i != _node_num; ++i) {
703 if ( _excess[i] > max_sup) max_sup = _excess[i];
704 if (-_excess[i] > max_dem) max_dem = -_excess[i];
707 for (int j = 0; j != _res_arc_num; ++j) {
708 if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
710 max_sup = std::min(std::min(max_sup, max_dem), max_cap);
712 for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2)
722 ProblemType start() {
723 // Execute the algorithm
726 pt = startWithScaling();
728 pt = startWithoutScaling();
730 // Handle non-zero lower bounds
732 for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) {
733 if (!_forward[j]) _res_cap[j] += _lower[j];
737 // Shift potentials if necessary
738 Cost pr = _pi[_root];
739 if (_sum_supply < 0 || pr > 0) {
740 for (int i = 0; i != _node_num; ++i) {
748 // Execute the capacity scaling algorithm
749 ProblemType startWithScaling() {
750 // Process capacity scaling phases
754 ResidualDijkstra _dijkstra(*this);
756 // Saturate all arcs not satisfying the optimality condition
757 for (int u = 0; u != _node_num; ++u) {
758 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) {
760 Cost c = _cost[a] + _pi[u] - _pi[v];
761 Value rc = _res_cap[a];
762 if (c < 0 && rc >= _delta) {
766 _res_cap[_reverse[a]] += rc;
771 // Find excess nodes and deficit nodes
772 _excess_nodes.clear();
773 _deficit_nodes.clear();
774 for (int u = 0; u != _node_num; ++u) {
775 if (_excess[u] >= _delta) _excess_nodes.push_back(u);
776 if (_excess[u] <= -_delta) _deficit_nodes.push_back(u);
778 int next_node = 0, next_def_node = 0;
780 // Find augmenting shortest paths
781 while (next_node < int(_excess_nodes.size())) {
782 // Check deficit nodes
784 bool delta_deficit = false;
785 for ( ; next_def_node < int(_deficit_nodes.size());
787 if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
788 delta_deficit = true;
792 if (!delta_deficit) break;
795 // Run Dijkstra in the residual network
796 s = _excess_nodes[next_node];
797 if ((t = _dijkstra.run(s, _delta)) == -1) {
805 // Augment along a shortest path from s to t
806 Value d = std::min(_excess[s], -_excess[t]);
810 while ((a = _pred[u]) != -1) {
811 if (_res_cap[a] < d) d = _res_cap[a];
816 while ((a = _pred[u]) != -1) {
818 _res_cap[_reverse[a]] += d;
824 if (_excess[s] < _delta) ++next_node;
827 if (_delta == 1) break;
828 if (++phase_cnt == _phase_num / 4) factor = 2;
829 _delta = _delta <= factor ? 1 : _delta / factor;
835 // Execute the successive shortest path algorithm
836 ProblemType startWithoutScaling() {
838 _excess_nodes.clear();
839 for (int i = 0; i != _node_num; ++i) {
840 if (_excess[i] > 0) _excess_nodes.push_back(i);
842 if (_excess_nodes.size() == 0) return OPTIMAL;
845 // Find shortest paths
847 ResidualDijkstra _dijkstra(*this);
848 while ( _excess[_excess_nodes[next_node]] > 0 ||
849 ++next_node < int(_excess_nodes.size()) )
851 // Run Dijkstra in the residual network
852 s = _excess_nodes[next_node];
853 if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
855 // Augment along a shortest path from s to t
856 Value d = std::min(_excess[s], -_excess[t]);
860 while ((a = _pred[u]) != -1) {
861 if (_res_cap[a] < d) d = _res_cap[a];
866 while ((a = _pred[u]) != -1) {
868 _res_cap[_reverse[a]] += d;
878 }; //class CapacityScaling
884 #endif //LEMON_CAPACITY_SCALING_H