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_algs
24 /// \brief Cost scaling algorithm for finding a minimum cost flow.
30 #include <lemon/core.h>
31 #include <lemon/maps.h>
32 #include <lemon/math.h>
33 #include <lemon/static_graph.h>
34 #include <lemon/circulation.h>
35 #include <lemon/bellman_ford.h>
39 /// \brief Default traits class of CostScaling algorithm.
41 /// Default traits class of CostScaling algorithm.
42 /// \tparam GR Digraph type.
43 /// \tparam V The number type used for flow amounts, capacity bounds
44 /// and supply values. By default it is \c int.
45 /// \tparam C The number type used for costs and potentials.
46 /// By default it is the same as \c V.
48 template <typename GR, typename V = int, typename C = V>
50 template < typename GR, typename V = int, typename C = V,
51 bool integer = std::numeric_limits<C>::is_integer >
53 struct CostScalingDefaultTraits
55 /// The type of the digraph
57 /// The type of the flow amounts, capacity bounds and supply values
59 /// The type of the arc costs
62 /// \brief The large cost type used for internal computations
64 /// The large cost type used for internal computations.
65 /// It is \c long \c long if the \c Cost type is integer,
66 /// otherwise it is \c double.
67 /// \c Cost must be convertible to \c LargeCost.
68 typedef double LargeCost;
71 // Default traits class for integer cost types
72 template <typename GR, typename V, typename C>
73 struct CostScalingDefaultTraits<GR, V, C, true>
78 #ifdef LEMON_HAVE_LONG_LONG
79 typedef long long LargeCost;
81 typedef long LargeCost;
86 /// \addtogroup min_cost_flow_algs
89 /// \brief Implementation of the Cost Scaling algorithm for
90 /// finding a \ref min_cost_flow "minimum cost flow".
92 /// \ref CostScaling implements a cost scaling algorithm that performs
93 /// push/augment and relabel operations for finding a minimum cost
94 /// flow. It is an efficient primal-dual solution method, which
95 /// can be viewed as the generalization of the \ref Preflow
96 /// "preflow push-relabel" algorithm for the maximum flow problem.
98 /// Most of the parameters of the problem (except for the digraph)
99 /// can be given using separate functions, and the algorithm can be
100 /// executed using the \ref run() function. If some parameters are not
101 /// specified, then default values will be used.
103 /// \tparam GR The digraph type the algorithm runs on.
104 /// \tparam V The number type used for flow amounts, capacity bounds
105 /// and supply values in the algorithm. By default it is \c int.
106 /// \tparam C The number type used for costs and potentials in the
107 /// algorithm. By default it is the same as \c V.
109 /// \warning Both number types must be signed and all input data must
111 /// \warning This algorithm does not support negative costs for such
112 /// arcs that have infinite upper bound.
114 /// \note %CostScaling provides three different internal methods,
115 /// from which the most efficient one is used by default.
116 /// For more information, see \ref Method.
118 template <typename GR, typename V, typename C, typename TR>
120 template < typename GR, typename V = int, typename C = V,
121 typename TR = CostScalingDefaultTraits<GR, V, C> >
127 /// The type of the digraph
128 typedef typename TR::Digraph Digraph;
129 /// The type of the flow amounts, capacity bounds and supply values
130 typedef typename TR::Value Value;
131 /// The type of the arc costs
132 typedef typename TR::Cost Cost;
134 /// \brief The large cost type
136 /// The large cost type used for internal computations.
137 /// Using the \ref CostScalingDefaultTraits "default traits class",
138 /// it is \c long \c long if the \c Cost type is integer,
139 /// otherwise it is \c double.
140 typedef typename TR::LargeCost LargeCost;
142 /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
147 /// \brief Problem type constants for the \c run() function.
149 /// Enum type containing the problem type constants that can be
150 /// returned by the \ref run() function of the algorithm.
152 /// The problem has no feasible solution (flow).
154 /// The problem has optimal solution (i.e. it is feasible and
155 /// bounded), and the algorithm has found optimal flow and node
156 /// potentials (primal and dual solutions).
158 /// The digraph contains an arc of negative cost and infinite
159 /// upper bound. It means that the objective function is unbounded
160 /// on that arc, however, note that it could actually be bounded
161 /// over the feasible flows, but this algroithm cannot handle
166 /// \brief Constants for selecting the internal method.
168 /// Enum type containing constants for selecting the internal method
169 /// for the \ref run() function.
171 /// \ref CostScaling provides three internal methods that differ mainly
172 /// in their base operations, which are used in conjunction with the
173 /// relabel operation.
174 /// By default, the so called \ref PARTIAL_AUGMENT
175 /// "Partial Augment-Relabel" method is used, which proved to be
176 /// the most efficient and the most robust on various test inputs.
177 /// However, the other methods can be selected using the \ref run()
178 /// function with the proper parameter.
180 /// Local push operations are used, i.e. flow is moved only on one
181 /// admissible arc at once.
183 /// Augment operations are used, i.e. flow is moved on admissible
184 /// paths from a node with excess to a node with deficit.
186 /// Partial augment operations are used, i.e. flow is moved on
187 /// admissible paths started from a node with excess, but the
188 /// lengths of these paths are limited. This method can be viewed
189 /// as a combined version of the previous two operations.
195 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
197 typedef std::vector<int> IntVector;
198 typedef std::vector<char> BoolVector;
199 typedef std::vector<Value> ValueVector;
200 typedef std::vector<Cost> CostVector;
201 typedef std::vector<LargeCost> LargeCostVector;
205 template <typename KT, typename VT>
211 VectorMap(std::vector<Value>& v) : _v(v) {}
213 const Value& operator[](const Key& key) const {
214 return _v[StaticDigraph::id(key)];
217 Value& operator[](const Key& key) {
218 return _v[StaticDigraph::id(key)];
221 void set(const Key& key, const Value& val) {
222 _v[StaticDigraph::id(key)] = val;
226 std::vector<Value>& _v;
229 typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
230 typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
234 // Data related to the underlying digraph
242 // Parameters of the problem
246 // Data structures for storing the digraph
250 IntVector _first_out;
262 ValueVector _res_cap;
263 LargeCostVector _cost;
267 std::deque<int> _active_nodes;
273 // Data for a StaticDigraph structure
274 typedef std::pair<int, int> IntPair;
276 std::vector<IntPair> _arc_vec;
277 std::vector<LargeCost> _cost_vec;
278 LargeCostArcMap _cost_map;
279 LargeCostNodeMap _pi_map;
283 /// \brief Constant for infinite upper bounds (capacities).
285 /// Constant for infinite upper bounds (capacities).
286 /// It is \c std::numeric_limits<Value>::infinity() if available,
287 /// \c std::numeric_limits<Value>::max() otherwise.
292 /// \name Named Template Parameters
295 template <typename T>
296 struct SetLargeCostTraits : public Traits {
300 /// \brief \ref named-templ-param "Named parameter" for setting
301 /// \c LargeCost type.
303 /// \ref named-templ-param "Named parameter" for setting \c LargeCost
304 /// type, which is used for internal computations in the algorithm.
305 /// \c Cost must be convertible to \c LargeCost.
306 template <typename T>
308 : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
309 typedef CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
316 /// \brief Constructor.
318 /// The constructor of the class.
320 /// \param graph The digraph the algorithm runs on.
321 CostScaling(const GR& graph) :
322 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
323 _cost_map(_cost_vec), _pi_map(_pi),
324 INF(std::numeric_limits<Value>::has_infinity ?
325 std::numeric_limits<Value>::infinity() :
326 std::numeric_limits<Value>::max())
328 // Check the number types
329 LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
330 "The flow type of CostScaling must be signed");
331 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
332 "The cost type of CostScaling must be signed");
335 _node_num = countNodes(_graph);
336 _arc_num = countArcs(_graph);
337 _res_node_num = _node_num + 1;
338 _res_arc_num = 2 * (_arc_num + _node_num);
341 _first_out.resize(_res_node_num + 1);
342 _forward.resize(_res_arc_num);
343 _source.resize(_res_arc_num);
344 _target.resize(_res_arc_num);
345 _reverse.resize(_res_arc_num);
347 _lower.resize(_res_arc_num);
348 _upper.resize(_res_arc_num);
349 _scost.resize(_res_arc_num);
350 _supply.resize(_res_node_num);
352 _res_cap.resize(_res_arc_num);
353 _cost.resize(_res_arc_num);
354 _pi.resize(_res_node_num);
355 _excess.resize(_res_node_num);
356 _next_out.resize(_res_node_num);
358 _arc_vec.reserve(_res_arc_num);
359 _cost_vec.reserve(_res_arc_num);
362 int i = 0, j = 0, k = 2 * _arc_num + _node_num;
363 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
367 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
369 for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
373 _target[j] = _node_id[_graph.runningNode(a)];
375 for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
379 _target[j] = _node_id[_graph.runningNode(a)];
392 _first_out[_res_node_num] = k;
393 for (ArcIt a(_graph); a != INVALID; ++a) {
394 int fi = _arc_idf[a];
395 int bi = _arc_idb[a];
405 /// The parameters of the algorithm can be specified using these
410 /// \brief Set the lower bounds on the arcs.
412 /// This function sets the lower bounds on the arcs.
413 /// If it is not used before calling \ref run(), the lower bounds
414 /// will be set to zero on all arcs.
416 /// \param map An arc map storing the lower 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 LowerMap>
422 CostScaling& lowerMap(const LowerMap& map) {
424 for (ArcIt a(_graph); a != INVALID; ++a) {
425 _lower[_arc_idf[a]] = map[a];
426 _lower[_arc_idb[a]] = map[a];
431 /// \brief Set the upper bounds (capacities) on the arcs.
433 /// This function sets the upper bounds (capacities) on the arcs.
434 /// If it is not used before calling \ref run(), the upper bounds
435 /// will be set to \ref INF on all arcs (i.e. the flow value will be
436 /// unbounded from above).
438 /// \param map An arc map storing the upper bounds.
439 /// Its \c Value type must be convertible to the \c Value type
440 /// of the algorithm.
442 /// \return <tt>(*this)</tt>
443 template<typename UpperMap>
444 CostScaling& upperMap(const UpperMap& map) {
445 for (ArcIt a(_graph); a != INVALID; ++a) {
446 _upper[_arc_idf[a]] = map[a];
451 /// \brief Set the costs of the arcs.
453 /// This function sets the costs of the arcs.
454 /// If it is not used before calling \ref run(), the costs
455 /// will be set to \c 1 on all arcs.
457 /// \param map An arc map storing the costs.
458 /// Its \c Value type must be convertible to the \c Cost type
459 /// of the algorithm.
461 /// \return <tt>(*this)</tt>
462 template<typename CostMap>
463 CostScaling& costMap(const CostMap& map) {
464 for (ArcIt a(_graph); a != INVALID; ++a) {
465 _scost[_arc_idf[a]] = map[a];
466 _scost[_arc_idb[a]] = -map[a];
471 /// \brief Set the supply values of the nodes.
473 /// This function sets the supply values of the nodes.
474 /// If neither this function nor \ref stSupply() is used before
475 /// calling \ref run(), the supply of each node will be set to zero.
477 /// \param map A node map storing the supply values.
478 /// Its \c Value type must be convertible to the \c Value type
479 /// of the algorithm.
481 /// \return <tt>(*this)</tt>
482 template<typename SupplyMap>
483 CostScaling& supplyMap(const SupplyMap& map) {
484 for (NodeIt n(_graph); n != INVALID; ++n) {
485 _supply[_node_id[n]] = map[n];
490 /// \brief Set single source and target nodes and a supply value.
492 /// This function sets a single source node and a single target node
493 /// and the required flow value.
494 /// If neither this function nor \ref supplyMap() is used before
495 /// calling \ref run(), the supply of each node will be set to zero.
497 /// Using this function has the same effect as using \ref supplyMap()
498 /// with such a map in which \c k is assigned to \c s, \c -k is
499 /// assigned to \c t and all other nodes have zero supply value.
501 /// \param s The source node.
502 /// \param t The target node.
503 /// \param k The required amount of flow from node \c s to node \c t
504 /// (i.e. the supply of \c s and the demand of \c t).
506 /// \return <tt>(*this)</tt>
507 CostScaling& stSupply(const Node& s, const Node& t, Value k) {
508 for (int i = 0; i != _res_node_num; ++i) {
511 _supply[_node_id[s]] = k;
512 _supply[_node_id[t]] = -k;
518 /// \name Execution control
519 /// The algorithm can be executed using \ref run().
523 /// \brief Run the algorithm.
525 /// This function runs the algorithm.
526 /// The paramters can be specified using functions \ref lowerMap(),
527 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
530 /// CostScaling<ListDigraph> cs(graph);
531 /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
532 /// .supplyMap(sup).run();
535 /// This function can be called more than once. All the parameters
536 /// that have been given are kept for the next call, unless
537 /// \ref reset() is called, thus only the modified parameters
538 /// have to be set again. See \ref reset() for examples.
539 /// However, the underlying digraph must not be modified after this
540 /// class have been constructed, since it copies and extends the graph.
542 /// \param method The internal method that will be used in the
543 /// algorithm. For more information, see \ref Method.
544 /// \param factor The cost scaling factor. It must be larger than one.
546 /// \return \c INFEASIBLE if no feasible flow exists,
547 /// \n \c OPTIMAL if the problem has optimal solution
548 /// (i.e. it is feasible and bounded), and the algorithm has found
549 /// optimal flow and node potentials (primal and dual solutions),
550 /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
551 /// and infinite upper bound. It means that the objective function
552 /// is unbounded on that arc, however, note that it could actually be
553 /// bounded over the feasible flows, but this algroithm cannot handle
556 /// \see ProblemType, Method
557 ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
559 ProblemType pt = init();
560 if (pt != OPTIMAL) return pt;
565 /// \brief Reset all the parameters that have been given before.
567 /// This function resets all the paramaters that have been given
568 /// before using functions \ref lowerMap(), \ref upperMap(),
569 /// \ref costMap(), \ref supplyMap(), \ref stSupply().
571 /// It is useful for multiple run() calls. If this function is not
572 /// used, all the parameters given before are kept for the next
574 /// However, the underlying digraph must not be modified after this
575 /// class have been constructed, since it copies and extends the graph.
579 /// CostScaling<ListDigraph> cs(graph);
582 /// cs.lowerMap(lower).upperMap(upper).costMap(cost)
583 /// .supplyMap(sup).run();
585 /// // Run again with modified cost map (reset() is not called,
586 /// // so only the cost map have to be set again)
588 /// cs.costMap(cost).run();
590 /// // Run again from scratch using reset()
591 /// // (the lower bounds will be set to zero on all arcs)
593 /// cs.upperMap(capacity).costMap(cost)
594 /// .supplyMap(sup).run();
597 /// \return <tt>(*this)</tt>
598 CostScaling& reset() {
599 for (int i = 0; i != _res_node_num; ++i) {
602 int limit = _first_out[_root];
603 for (int j = 0; j != limit; ++j) {
606 _scost[j] = _forward[j] ? 1 : -1;
608 for (int j = limit; j != _res_arc_num; ++j) {
612 _scost[_reverse[j]] = 0;
620 /// \name Query Functions
621 /// The results of the algorithm can be obtained using these
623 /// The \ref run() function must be called before using them.
627 /// \brief Return the total cost of the found flow.
629 /// This function returns the total cost of the found flow.
630 /// Its complexity is O(e).
632 /// \note The return type of the function can be specified as a
633 /// template parameter. For example,
635 /// cs.totalCost<double>();
637 /// It is useful if the total cost cannot be stored in the \c Cost
638 /// type of the algorithm, which is the default return type of the
641 /// \pre \ref run() must be called before using this function.
642 template <typename Number>
643 Number totalCost() const {
645 for (ArcIt a(_graph); a != INVALID; ++a) {
647 c += static_cast<Number>(_res_cap[i]) *
648 (-static_cast<Number>(_scost[i]));
654 Cost totalCost() const {
655 return totalCost<Cost>();
659 /// \brief Return the flow on the given arc.
661 /// This function returns the flow on the given arc.
663 /// \pre \ref run() must be called before using this function.
664 Value flow(const Arc& a) const {
665 return _res_cap[_arc_idb[a]];
668 /// \brief Return the flow map (the primal solution).
670 /// This function copies the flow value on each arc into the given
671 /// map. The \c Value type of the algorithm must be convertible to
672 /// the \c Value type of the map.
674 /// \pre \ref run() must be called before using this function.
675 template <typename FlowMap>
676 void flowMap(FlowMap &map) const {
677 for (ArcIt a(_graph); a != INVALID; ++a) {
678 map.set(a, _res_cap[_arc_idb[a]]);
682 /// \brief Return the potential (dual value) of the given node.
684 /// This function returns the potential (dual value) of the
687 /// \pre \ref run() must be called before using this function.
688 Cost potential(const Node& n) const {
689 return static_cast<Cost>(_pi[_node_id[n]]);
692 /// \brief Return the potential map (the dual solution).
694 /// This function copies the potential (dual value) of each node
695 /// into the given map.
696 /// The \c Cost type of the algorithm must be convertible to the
697 /// \c Value type of the map.
699 /// \pre \ref run() must be called before using this function.
700 template <typename PotentialMap>
701 void potentialMap(PotentialMap &map) const {
702 for (NodeIt n(_graph); n != INVALID; ++n) {
703 map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
711 // Initialize the algorithm
713 if (_res_node_num == 0) return INFEASIBLE;
715 // Check the sum of supply values
717 for (int i = 0; i != _root; ++i) {
718 _sum_supply += _supply[i];
720 if (_sum_supply > 0) return INFEASIBLE;
723 // Initialize vectors
724 for (int i = 0; i != _res_node_num; ++i) {
726 _excess[i] = _supply[i];
729 // Remove infinite upper bounds and check negative arcs
730 const Value MAX = std::numeric_limits<Value>::max();
733 for (int i = 0; i != _root; ++i) {
734 last_out = _first_out[i+1];
735 for (int j = _first_out[i]; j != last_out; ++j) {
737 Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
738 if (c >= MAX) return UNBOUNDED;
740 _excess[_target[j]] += c;
745 for (int i = 0; i != _root; ++i) {
746 last_out = _first_out[i+1];
747 for (int j = _first_out[i]; j != last_out; ++j) {
748 if (_forward[j] && _scost[j] < 0) {
750 if (c >= MAX) return UNBOUNDED;
752 _excess[_target[j]] += c;
757 Value ex, max_cap = 0;
758 for (int i = 0; i != _res_node_num; ++i) {
761 if (ex < 0) max_cap -= ex;
763 for (int j = 0; j != _res_arc_num; ++j) {
764 if (_upper[j] >= MAX) _upper[j] = max_cap;
767 // Initialize the large cost vector and the epsilon parameter
770 for (int i = 0; i != _root; ++i) {
771 last_out = _first_out[i+1];
772 for (int j = _first_out[i]; j != last_out; ++j) {
773 lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
775 if (lc > _epsilon) _epsilon = lc;
780 // Initialize maps for Circulation and remove non-zero lower bounds
781 ConstMap<Arc, Value> low(0);
782 typedef typename Digraph::template ArcMap<Value> ValueArcMap;
783 typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
784 ValueArcMap cap(_graph), flow(_graph);
785 ValueNodeMap sup(_graph);
786 for (NodeIt n(_graph); n != INVALID; ++n) {
787 sup[n] = _supply[_node_id[n]];
790 for (ArcIt a(_graph); a != INVALID; ++a) {
793 cap[a] = _upper[j] - c;
794 sup[_graph.source(a)] -= c;
795 sup[_graph.target(a)] += c;
798 for (ArcIt a(_graph); a != INVALID; ++a) {
799 cap[a] = _upper[_arc_idf[a]];
803 // Find a feasible flow using Circulation
804 Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
805 circ(_graph, low, cap, sup);
806 if (!circ.flowMap(flow).run()) return INFEASIBLE;
808 // Set residual capacities and handle GEQ supply type
809 if (_sum_supply < 0) {
810 for (ArcIt a(_graph); a != INVALID; ++a) {
812 _res_cap[_arc_idf[a]] = cap[a] - fa;
813 _res_cap[_arc_idb[a]] = fa;
814 sup[_graph.source(a)] -= fa;
815 sup[_graph.target(a)] += fa;
817 for (NodeIt n(_graph); n != INVALID; ++n) {
818 _excess[_node_id[n]] = sup[n];
820 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
822 int ra = _reverse[a];
823 _res_cap[a] = -_sum_supply + 1;
824 _res_cap[ra] = -_excess[u];
830 for (ArcIt a(_graph); a != INVALID; ++a) {
832 _res_cap[_arc_idf[a]] = cap[a] - fa;
833 _res_cap[_arc_idb[a]] = fa;
835 for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
836 int ra = _reverse[a];
847 // Execute the algorithm and transform the results
848 void start(Method method) {
849 // Maximum path length for partial augment
850 const int MAX_PATH_LENGTH = 4;
852 // Execute the algorithm
860 case PARTIAL_AUGMENT:
861 startAugment(MAX_PATH_LENGTH);
865 // Compute node potentials for the original costs
868 for (int j = 0; j != _res_arc_num; ++j) {
869 if (_res_cap[j] > 0) {
870 _arc_vec.push_back(IntPair(_source[j], _target[j]));
871 _cost_vec.push_back(_scost[j]);
874 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
876 typename BellmanFord<StaticDigraph, LargeCostArcMap>
877 ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
882 // Handle non-zero lower bounds
884 int limit = _first_out[_root];
885 for (int j = 0; j != limit; ++j) {
886 if (!_forward[j]) _res_cap[j] += _lower[j];
891 /// Execute the algorithm performing augment and relabel operations
892 void startAugment(int max_length = std::numeric_limits<int>::max()) {
893 // Paramters for heuristics
894 const int BF_HEURISTIC_EPSILON_BOUND = 1000;
895 const int BF_HEURISTIC_BOUND_FACTOR = 3;
897 // Perform cost scaling phases
898 IntVector pred_arc(_res_node_num);
899 std::vector<int> path_nodes;
900 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
901 1 : _epsilon / _alpha )
903 // "Early Termination" heuristic: use Bellman-Ford algorithm
904 // to check if the current flow is optimal
905 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
908 for (int j = 0; j != _res_arc_num; ++j) {
909 if (_res_cap[j] > 0) {
910 _arc_vec.push_back(IntPair(_source[j], _target[j]));
911 _cost_vec.push_back(_cost[j] + 1);
914 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
916 BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
919 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
920 for (int i = 0; i < K && !done; ++i)
921 done = bf.processNextWeakRound();
925 // Saturate arcs not satisfying the optimality condition
926 for (int a = 0; a != _res_arc_num; ++a) {
927 if (_res_cap[a] > 0 &&
928 _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
929 Value delta = _res_cap[a];
930 _excess[_source[a]] -= delta;
931 _excess[_target[a]] += delta;
933 _res_cap[_reverse[a]] += delta;
937 // Find active nodes (i.e. nodes with positive excess)
938 for (int u = 0; u != _res_node_num; ++u) {
939 if (_excess[u] > 0) _active_nodes.push_back(u);
942 // Initialize the next arcs
943 for (int u = 0; u != _res_node_num; ++u) {
944 _next_out[u] = _first_out[u];
947 // Perform partial augment and relabel operations
949 // Select an active node (FIFO selection)
950 while (_active_nodes.size() > 0 &&
951 _excess[_active_nodes.front()] <= 0) {
952 _active_nodes.pop_front();
954 if (_active_nodes.size() == 0) break;
955 int start = _active_nodes.front();
957 path_nodes.push_back(start);
959 // Find an augmenting path from the start node
961 while (_excess[tip] >= 0 &&
962 int(path_nodes.size()) <= max_length) {
964 LargeCost min_red_cost, rc;
965 int last_out = _sum_supply < 0 ?
966 _first_out[tip+1] : _first_out[tip+1] - 1;
967 for (int a = _next_out[tip]; a != last_out; ++a) {
968 if (_res_cap[a] > 0 &&
969 _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
974 path_nodes.push_back(tip);
980 min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
981 for (int a = _first_out[tip]; a != last_out; ++a) {
982 rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
983 if (_res_cap[a] > 0 && rc < min_red_cost) {
987 _pi[tip] -= min_red_cost + _epsilon;
989 // Reset the next arc of tip
990 _next_out[tip] = _first_out[tip];
994 path_nodes.pop_back();
995 tip = path_nodes.back();
1001 // Augment along the found path (as much flow as possible)
1003 int u, v = path_nodes.front(), pa;
1004 for (int i = 1; i < int(path_nodes.size()); ++i) {
1008 delta = std::min(_res_cap[pa], _excess[u]);
1009 _res_cap[pa] -= delta;
1010 _res_cap[_reverse[pa]] += delta;
1011 _excess[u] -= delta;
1012 _excess[v] += delta;
1013 if (_excess[v] > 0 && _excess[v] <= delta)
1014 _active_nodes.push_back(v);
1020 /// Execute the algorithm performing push and relabel operations
1022 // Paramters for heuristics
1023 const int BF_HEURISTIC_EPSILON_BOUND = 1000;
1024 const int BF_HEURISTIC_BOUND_FACTOR = 3;
1026 // Perform cost scaling phases
1027 BoolVector hyper(_res_node_num, false);
1028 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1029 1 : _epsilon / _alpha )
1031 // "Early Termination" heuristic: use Bellman-Ford algorithm
1032 // to check if the current flow is optimal
1033 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
1036 for (int j = 0; j != _res_arc_num; ++j) {
1037 if (_res_cap[j] > 0) {
1038 _arc_vec.push_back(IntPair(_source[j], _target[j]));
1039 _cost_vec.push_back(_cost[j] + 1);
1042 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
1044 BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
1047 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
1048 for (int i = 0; i < K && !done; ++i)
1049 done = bf.processNextWeakRound();
1053 // Saturate arcs not satisfying the optimality condition
1054 for (int a = 0; a != _res_arc_num; ++a) {
1055 if (_res_cap[a] > 0 &&
1056 _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1057 Value delta = _res_cap[a];
1058 _excess[_source[a]] -= delta;
1059 _excess[_target[a]] += delta;
1061 _res_cap[_reverse[a]] += delta;
1065 // Find active nodes (i.e. nodes with positive excess)
1066 for (int u = 0; u != _res_node_num; ++u) {
1067 if (_excess[u] > 0) _active_nodes.push_back(u);
1070 // Initialize the next arcs
1071 for (int u = 0; u != _res_node_num; ++u) {
1072 _next_out[u] = _first_out[u];
1075 // Perform push and relabel operations
1076 while (_active_nodes.size() > 0) {
1077 LargeCost min_red_cost, rc;
1079 int n, t, a, last_out = _res_arc_num;
1081 // Select an active node (FIFO selection)
1083 n = _active_nodes.front();
1084 last_out = _sum_supply < 0 ?
1085 _first_out[n+1] : _first_out[n+1] - 1;
1087 // Perform push operations if there are admissible arcs
1088 if (_excess[n] > 0) {
1089 for (a = _next_out[n]; a != last_out; ++a) {
1090 if (_res_cap[a] > 0 &&
1091 _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
1092 delta = std::min(_res_cap[a], _excess[n]);
1095 // Push-look-ahead heuristic
1096 Value ahead = -_excess[t];
1097 int last_out_t = _sum_supply < 0 ?
1098 _first_out[t+1] : _first_out[t+1] - 1;
1099 for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
1100 if (_res_cap[ta] > 0 &&
1101 _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
1102 ahead += _res_cap[ta];
1103 if (ahead >= delta) break;
1105 if (ahead < 0) ahead = 0;
1107 // Push flow along the arc
1108 if (ahead < delta) {
1109 _res_cap[a] -= ahead;
1110 _res_cap[_reverse[a]] += ahead;
1111 _excess[n] -= ahead;
1112 _excess[t] += ahead;
1113 _active_nodes.push_front(t);
1118 _res_cap[a] -= delta;
1119 _res_cap[_reverse[a]] += delta;
1120 _excess[n] -= delta;
1121 _excess[t] += delta;
1122 if (_excess[t] > 0 && _excess[t] <= delta)
1123 _active_nodes.push_back(t);
1126 if (_excess[n] == 0) {
1135 // Relabel the node if it is still active (or hyper)
1136 if (_excess[n] > 0 || hyper[n]) {
1137 min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
1138 for (int a = _first_out[n]; a != last_out; ++a) {
1139 rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
1140 if (_res_cap[a] > 0 && rc < min_red_cost) {
1144 _pi[n] -= min_red_cost + _epsilon;
1147 // Reset the next arc
1148 _next_out[n] = _first_out[n];
1151 // Remove nodes that are not active nor hyper
1153 while ( _active_nodes.size() > 0 &&
1154 _excess[_active_nodes.front()] <= 0 &&
1155 !hyper[_active_nodes.front()] ) {
1156 _active_nodes.pop_front();
1162 }; //class CostScaling
1168 #endif //LEMON_COST_SCALING_H