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
24 /// \brief Cost scaling algorithm for finding a minimum cost flow.
27 #include <lemon/graph_adaptor.h>
28 #include <lemon/graph_utils.h>
29 #include <lemon/maps.h>
30 #include <lemon/math.h>
32 #include <lemon/circulation.h>
33 #include <lemon/bellman_ford.h>
37 /// \addtogroup min_cost_flow
40 /// \brief Implementation of the cost scaling algorithm for finding a
41 /// minimum cost flow.
43 /// \ref CostScaling implements the cost scaling algorithm performing
44 /// augment/push and relabel operations for finding a minimum cost
47 /// \tparam Graph The directed graph type the algorithm runs on.
48 /// \tparam LowerMap The type of the lower bound map.
49 /// \tparam CapacityMap The type of the capacity (upper bound) map.
50 /// \tparam CostMap The type of the cost (length) map.
51 /// \tparam SupplyMap The type of the supply map.
54 /// - Edge capacities and costs should be \e non-negative \e integers.
55 /// - Supply values should be \e signed \e integers.
56 /// - The value types of the maps should be convertible to each other.
57 /// - \c CostMap::Value must be signed type.
59 /// \note Edge costs are multiplied with the number of nodes during
60 /// the algorithm so overflow problems may arise more easily than with
61 /// other minimum cost flow algorithms.
62 /// If it is available, <tt>long long int</tt> type is used instead of
63 /// <tt>long int</tt> in the inside computations.
65 /// \author Peter Kovacs
66 template < typename Graph,
67 typename LowerMap = typename Graph::template EdgeMap<int>,
68 typename CapacityMap = typename Graph::template EdgeMap<int>,
69 typename CostMap = typename Graph::template EdgeMap<int>,
70 typename SupplyMap = typename Graph::template NodeMap<int> >
73 GRAPH_TYPEDEFS(typename Graph);
75 typedef typename CapacityMap::Value Capacity;
76 typedef typename CostMap::Value Cost;
77 typedef typename SupplyMap::Value Supply;
78 typedef typename Graph::template EdgeMap<Capacity> CapacityEdgeMap;
79 typedef typename Graph::template NodeMap<Supply> SupplyNodeMap;
81 typedef ResGraphAdaptor< const Graph, Capacity,
82 CapacityEdgeMap, CapacityEdgeMap > ResGraph;
83 typedef typename ResGraph::Edge ResEdge;
85 #if defined __GNUC__ && !defined __STRICT_ANSI__
86 typedef long long int LCost;
88 typedef long int LCost;
90 typedef typename Graph::template EdgeMap<LCost> LargeCostMap;
94 /// The type of the flow map.
95 typedef typename Graph::template EdgeMap<Capacity> FlowMap;
96 /// The type of the potential map.
97 typedef typename Graph::template NodeMap<LCost> PotentialMap;
101 /// \brief Map adaptor class for handling residual edge costs.
103 /// Map adaptor class for handling residual edge costs.
104 template <typename Map>
105 class ResidualCostMap : public MapBase<ResEdge, typename Map::Value>
109 const Map &_cost_map;
114 ResidualCostMap(const Map &cost_map) :
115 _cost_map(cost_map) {}
118 inline typename Map::Value operator[](const ResEdge &e) const {
119 return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
122 }; //class ResidualCostMap
124 /// \brief Map adaptor class for handling reduced edge costs.
126 /// Map adaptor class for handling reduced edge costs.
127 class ReducedCostMap : public MapBase<Edge, LCost>
132 const LargeCostMap &_cost_map;
133 const PotentialMap &_pot_map;
138 ReducedCostMap( const Graph &gr,
139 const LargeCostMap &cost_map,
140 const PotentialMap &pot_map ) :
141 _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
144 inline LCost operator[](const Edge &e) const {
145 return _cost_map[e] + _pot_map[_gr.source(e)]
146 - _pot_map[_gr.target(e)];
149 }; //class ReducedCostMap
153 // The directed graph the algorithm runs on
155 // The original lower bound map
156 const LowerMap *_lower;
157 // The modified capacity map
158 CapacityEdgeMap _capacity;
159 // The original cost map
160 const CostMap &_orig_cost;
161 // The scaled cost map
163 // The modified supply map
164 SupplyNodeMap _supply;
167 // Edge map of the current flow
170 // Node map of the current potentials
171 PotentialMap *_potential;
172 bool _local_potential;
174 // The residual cost map
175 ResidualCostMap<LargeCostMap> _res_cost;
176 // The residual graph
177 ResGraph *_res_graph;
178 // The reduced cost map
179 ReducedCostMap *_red_cost;
181 SupplyNodeMap _excess;
182 // The epsilon parameter used for cost scaling
184 // The scaling factor
189 /// \brief General constructor (with lower bounds).
191 /// General constructor (with lower bounds).
193 /// \param graph The directed graph the algorithm runs on.
194 /// \param lower The lower bounds of the edges.
195 /// \param capacity The capacities (upper bounds) of the edges.
196 /// \param cost The cost (length) values of the edges.
197 /// \param supply The supply values of the nodes (signed).
198 CostScaling( const Graph &graph,
199 const LowerMap &lower,
200 const CapacityMap &capacity,
202 const SupplyMap &supply ) :
203 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
204 _cost(graph), _supply(graph), _flow(NULL), _local_flow(false),
205 _potential(NULL), _local_potential(false), _res_cost(_cost),
206 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
208 // Remove non-zero lower bounds
209 _capacity = subMap(capacity, lower);
211 for (NodeIt n(_graph); n != INVALID; ++n) {
212 Supply s = supply[n];
213 for (InEdgeIt e(_graph, n); e != INVALID; ++e)
215 for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
220 _valid_supply = sum == 0;
223 /// \brief General constructor (without lower bounds).
225 /// General constructor (without lower bounds).
227 /// \param graph The directed graph the algorithm runs on.
228 /// \param capacity The capacities (upper bounds) of the edges.
229 /// \param cost The cost (length) values of the edges.
230 /// \param supply The supply values of the nodes (signed).
231 CostScaling( const Graph &graph,
232 const CapacityMap &capacity,
234 const SupplyMap &supply ) :
235 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
236 _cost(graph), _supply(supply), _flow(NULL), _local_flow(false),
237 _potential(NULL), _local_potential(false), _res_cost(_cost),
238 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
240 // Check the sum of supply values
242 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
243 _valid_supply = sum == 0;
246 /// \brief Simple constructor (with lower bounds).
248 /// Simple constructor (with lower bounds).
250 /// \param graph The directed graph the algorithm runs on.
251 /// \param lower The lower bounds of the edges.
252 /// \param capacity The capacities (upper bounds) of the edges.
253 /// \param cost The cost (length) values of the edges.
254 /// \param s The source node.
255 /// \param t The target node.
256 /// \param flow_value The required amount of flow from node \c s
257 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
258 CostScaling( const Graph &graph,
259 const LowerMap &lower,
260 const CapacityMap &capacity,
263 Supply flow_value ) :
264 _graph(graph), _lower(&lower), _capacity(graph), _orig_cost(cost),
265 _cost(graph), _supply(graph), _flow(NULL), _local_flow(false),
266 _potential(NULL), _local_potential(false), _res_cost(_cost),
267 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
269 // Remove nonzero lower bounds
270 _capacity = subMap(capacity, lower);
271 for (NodeIt n(_graph); n != INVALID; ++n) {
273 if (n == s) sum = flow_value;
274 if (n == t) sum = -flow_value;
275 for (InEdgeIt e(_graph, n); e != INVALID; ++e)
277 for (OutEdgeIt e(_graph, n); e != INVALID; ++e)
281 _valid_supply = true;
284 /// \brief Simple constructor (without lower bounds).
286 /// Simple constructor (without lower bounds).
288 /// \param graph The directed graph the algorithm runs on.
289 /// \param capacity The capacities (upper bounds) of the edges.
290 /// \param cost The cost (length) values of the edges.
291 /// \param s The source node.
292 /// \param t The target node.
293 /// \param flow_value The required amount of flow from node \c s
294 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
295 CostScaling( const Graph &graph,
296 const CapacityMap &capacity,
299 Supply flow_value ) :
300 _graph(graph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
301 _cost(graph), _supply(graph, 0), _flow(NULL), _local_flow(false),
302 _potential(NULL), _local_potential(false), _res_cost(_cost),
303 _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
305 _supply[s] = flow_value;
306 _supply[t] = -flow_value;
307 _valid_supply = true;
312 if (_local_flow) delete _flow;
313 if (_local_potential) delete _potential;
318 /// \brief Set the flow map.
320 /// Set the flow map.
322 /// \return \c (*this)
323 CostScaling& flowMap(FlowMap &map) {
332 /// \brief Set the potential map.
334 /// Set the potential map.
336 /// \return \c (*this)
337 CostScaling& potentialMap(PotentialMap &map) {
338 if (_local_potential) {
340 _local_potential = false;
346 /// \name Execution control
350 /// \brief Run the algorithm.
352 /// Run the algorithm.
354 /// \param partial_augment By default the algorithm performs
355 /// partial augment and relabel operations in the cost scaling
356 /// phases. Set this parameter to \c false for using local push and
357 /// relabel operations instead.
359 /// \return \c true if a feasible flow can be found.
360 bool run(bool partial_augment = true) {
361 if (partial_augment) {
362 return init() && startPartialAugment();
364 return init() && startPushRelabel();
370 /// \name Query Functions
371 /// The result of the algorithm can be obtained using these
373 /// \ref lemon::CostScaling::run() "run()" must be called before
378 /// \brief Return a const reference to the edge map storing the
381 /// Return a const reference to the edge map storing the found flow.
383 /// \pre \ref run() must be called before using this function.
384 const FlowMap& flowMap() const {
388 /// \brief Return a const reference to the node map storing the
389 /// found potentials (the dual solution).
391 /// Return a const reference to the node map storing the found
392 /// potentials (the dual solution).
394 /// \pre \ref run() must be called before using this function.
395 const PotentialMap& potentialMap() const {
399 /// \brief Return the flow on the given edge.
401 /// Return the flow on the given edge.
403 /// \pre \ref run() must be called before using this function.
404 Capacity flow(const Edge& edge) const {
405 return (*_flow)[edge];
408 /// \brief Return the potential of the given node.
410 /// Return the potential of the given node.
412 /// \pre \ref run() must be called before using this function.
413 Cost potential(const Node& node) const {
414 return (*_potential)[node];
417 /// \brief Return the total cost of the found flow.
419 /// Return the total cost of the found flow. The complexity of the
420 /// function is \f$ O(e) \f$.
422 /// \pre \ref run() must be called before using this function.
423 Cost totalCost() const {
425 for (EdgeIt e(_graph); e != INVALID; ++e)
426 c += (*_flow)[e] * _orig_cost[e];
434 /// Initialize the algorithm.
436 if (!_valid_supply) return false;
437 // The scaling factor
440 // Initialize flow and potential maps
442 _flow = new FlowMap(_graph);
446 _potential = new PotentialMap(_graph);
447 _local_potential = true;
450 _red_cost = new ReducedCostMap(_graph, _cost, *_potential);
451 _res_graph = new ResGraph(_graph, _capacity, *_flow);
453 // Initialize the scaled cost map and the epsilon parameter
455 int node_num = countNodes(_graph);
456 for (EdgeIt e(_graph); e != INVALID; ++e) {
457 _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
458 if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
460 _epsilon = max_cost * node_num;
462 // Find a feasible flow using Circulation
463 Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
465 circulation( _graph, constMap<Edge>(Capacity(0)), _capacity,
467 return circulation.flowMap(*_flow).run();
470 /// Execute the algorithm performing partial augmentation and
471 /// relabel operations.
472 bool startPartialAugment() {
473 // Paramters for heuristics
474 const int BF_HEURISTIC_EPSILON_BOUND = 1000;
475 const int BF_HEURISTIC_BOUND_FACTOR = 3;
476 // Maximum augment path length
477 const int MAX_PATH_LENGTH = 4;
480 typename Graph::template NodeMap<Edge> pred_edge(_graph);
481 typename Graph::template NodeMap<bool> forward(_graph);
482 typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
483 typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
484 typename Graph::template NodeMap<bool> next_dir(_graph);
485 std::deque<Node> active_nodes;
486 std::vector<Node> path_nodes;
488 int node_num = countNodes(_graph);
489 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
490 1 : _epsilon / _alpha )
492 // "Early Termination" heuristic: use Bellman-Ford algorithm
493 // to check if the current flow is optimal
494 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
495 typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
496 ShiftCostMap shift_cost(_res_cost, 1);
497 BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
500 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
501 for (int i = 0; i < K && !done; ++i)
502 done = bf.processNextWeakRound();
506 // Saturate edges not satisfying the optimality condition
508 for (EdgeIt e(_graph); e != INVALID; ++e) {
509 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
510 delta = _capacity[e] - (*_flow)[e];
511 _excess[_graph.source(e)] -= delta;
512 _excess[_graph.target(e)] += delta;
513 (*_flow)[e] = _capacity[e];
515 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
516 _excess[_graph.target(e)] -= (*_flow)[e];
517 _excess[_graph.source(e)] += (*_flow)[e];
522 // Find active nodes (i.e. nodes with positive excess)
523 for (NodeIt n(_graph); n != INVALID; ++n) {
524 if (_excess[n] > 0) active_nodes.push_back(n);
527 // Initialize the next edge maps
528 for (NodeIt n(_graph); n != INVALID; ++n) {
529 next_out[n] = OutEdgeIt(_graph, n);
530 next_in[n] = InEdgeIt(_graph, n);
534 // Perform partial augment and relabel operations
535 while (active_nodes.size() > 0) {
536 // Select an active node (FIFO selection)
537 if (_excess[active_nodes[0]] <= 0) {
538 active_nodes.pop_front();
541 Node start = active_nodes[0];
543 path_nodes.push_back(start);
545 // Find an augmenting path from the start node
548 while ( _excess[tip] >= 0 &&
549 int(path_nodes.size()) <= MAX_PATH_LENGTH )
552 for (OutEdgeIt e = next_out[tip]; e != INVALID; ++e) {
553 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
554 u = _graph.target(e);
559 path_nodes.push_back(tip);
563 next_dir[tip] = false;
565 for (InEdgeIt e = next_in[tip]; e != INVALID; ++e) {
566 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
567 u = _graph.source(e);
572 path_nodes.push_back(tip);
578 min_red_cost = std::numeric_limits<LCost>::max() / 2;
579 for (OutEdgeIt oe(_graph, tip); oe != INVALID; ++oe) {
580 if ( _capacity[oe] - (*_flow)[oe] > 0 &&
581 (*_red_cost)[oe] < min_red_cost )
582 min_red_cost = (*_red_cost)[oe];
584 for (InEdgeIt ie(_graph, tip); ie != INVALID; ++ie) {
585 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
586 min_red_cost = -(*_red_cost)[ie];
588 (*_potential)[tip] -= min_red_cost + _epsilon;
590 // Reset the next edge maps
591 next_out[tip] = OutEdgeIt(_graph, tip);
592 next_in[tip] = InEdgeIt(_graph, tip);
593 next_dir[tip] = true;
597 path_nodes.pop_back();
598 tip = path_nodes[path_nodes.size()-1];
605 // Augment along the found path (as much flow as possible)
607 for (int i = 1; i < int(path_nodes.size()); ++i) {
610 _capacity[pred_edge[u]] - (*_flow)[pred_edge[u]] :
611 (*_flow)[pred_edge[u]];
612 delta = std::min(delta, _excess[path_nodes[i-1]]);
613 (*_flow)[pred_edge[u]] += forward[u] ? delta : -delta;
614 _excess[path_nodes[i-1]] -= delta;
616 if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
621 // Compute node potentials for the original costs
622 ResidualCostMap<CostMap> res_cost(_orig_cost);
623 BellmanFord< ResGraph, ResidualCostMap<CostMap> >
624 bf(*_res_graph, res_cost);
625 bf.init(0); bf.start();
626 for (NodeIt n(_graph); n != INVALID; ++n)
627 (*_potential)[n] = bf.dist(n);
629 // Handle non-zero lower bounds
631 for (EdgeIt e(_graph); e != INVALID; ++e)
632 (*_flow)[e] += (*_lower)[e];
637 /// Execute the algorithm performing push and relabel operations.
638 bool startPushRelabel() {
639 // Paramters for heuristics
640 const int BF_HEURISTIC_EPSILON_BOUND = 1000;
641 const int BF_HEURISTIC_BOUND_FACTOR = 3;
643 typename Graph::template NodeMap<bool> hyper(_graph, false);
644 typename Graph::template NodeMap<Edge> pred_edge(_graph);
645 typename Graph::template NodeMap<bool> forward(_graph);
646 typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
647 typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
648 typename Graph::template NodeMap<bool> next_dir(_graph);
649 std::deque<Node> active_nodes;
651 int node_num = countNodes(_graph);
652 for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
653 1 : _epsilon / _alpha )
655 // "Early Termination" heuristic: use Bellman-Ford algorithm
656 // to check if the current flow is optimal
657 if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
658 typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
659 ShiftCostMap shift_cost(_res_cost, 1);
660 BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
663 int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
664 for (int i = 0; i < K && !done; ++i)
665 done = bf.processNextWeakRound();
669 // Saturate edges not satisfying the optimality condition
671 for (EdgeIt e(_graph); e != INVALID; ++e) {
672 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
673 delta = _capacity[e] - (*_flow)[e];
674 _excess[_graph.source(e)] -= delta;
675 _excess[_graph.target(e)] += delta;
676 (*_flow)[e] = _capacity[e];
678 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
679 _excess[_graph.target(e)] -= (*_flow)[e];
680 _excess[_graph.source(e)] += (*_flow)[e];
685 // Find active nodes (i.e. nodes with positive excess)
686 for (NodeIt n(_graph); n != INVALID; ++n) {
687 if (_excess[n] > 0) active_nodes.push_back(n);
690 // Initialize the next edge maps
691 for (NodeIt n(_graph); n != INVALID; ++n) {
692 next_out[n] = OutEdgeIt(_graph, n);
693 next_in[n] = InEdgeIt(_graph, n);
697 // Perform push and relabel operations
698 while (active_nodes.size() > 0) {
699 // Select an active node (FIFO selection)
700 Node n = active_nodes[0], t;
701 bool relabel_enabled = true;
703 // Perform push operations if there are admissible edges
704 if (_excess[n] > 0 && next_dir[n]) {
705 OutEdgeIt e = next_out[n];
706 for ( ; e != INVALID; ++e) {
707 if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
708 delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
709 t = _graph.target(e);
711 // Push-look-ahead heuristic
712 Capacity ahead = -_excess[t];
713 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
714 if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
715 ahead += _capacity[oe] - (*_flow)[oe];
717 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
718 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
719 ahead += (*_flow)[ie];
721 if (ahead < 0) ahead = 0;
723 // Push flow along the edge
725 (*_flow)[e] += ahead;
728 active_nodes.push_front(t);
730 relabel_enabled = false;
733 (*_flow)[e] += delta;
736 if (_excess[t] > 0 && _excess[t] <= delta)
737 active_nodes.push_back(t);
740 if (_excess[n] == 0) break;
750 if (_excess[n] > 0 && !next_dir[n]) {
751 InEdgeIt e = next_in[n];
752 for ( ; e != INVALID; ++e) {
753 if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
754 delta = std::min((*_flow)[e], _excess[n]);
755 t = _graph.source(e);
757 // Push-look-ahead heuristic
758 Capacity ahead = -_excess[t];
759 for (OutEdgeIt oe(_graph, t); oe != INVALID; ++oe) {
760 if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
761 ahead += _capacity[oe] - (*_flow)[oe];
763 for (InEdgeIt ie(_graph, t); ie != INVALID; ++ie) {
764 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
765 ahead += (*_flow)[ie];
767 if (ahead < 0) ahead = 0;
769 // Push flow along the edge
771 (*_flow)[e] -= ahead;
774 active_nodes.push_front(t);
776 relabel_enabled = false;
779 (*_flow)[e] -= delta;
782 if (_excess[t] > 0 && _excess[t] <= delta)
783 active_nodes.push_back(t);
786 if (_excess[n] == 0) break;
792 // Relabel the node if it is still active (or hyper)
793 if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
794 LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
795 for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
796 if ( _capacity[oe] - (*_flow)[oe] > 0 &&
797 (*_red_cost)[oe] < min_red_cost )
798 min_red_cost = (*_red_cost)[oe];
800 for (InEdgeIt ie(_graph, n); ie != INVALID; ++ie) {
801 if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
802 min_red_cost = -(*_red_cost)[ie];
804 (*_potential)[n] -= min_red_cost + _epsilon;
807 // Reset the next edge maps
808 next_out[n] = OutEdgeIt(_graph, n);
809 next_in[n] = InEdgeIt(_graph, n);
813 // Remove nodes that are not active nor hyper
814 while ( active_nodes.size() > 0 &&
815 _excess[active_nodes[0]] <= 0 &&
816 !hyper[active_nodes[0]] ) {
817 active_nodes.pop_front();
822 // Compute node potentials for the original costs
823 ResidualCostMap<CostMap> res_cost(_orig_cost);
824 BellmanFord< ResGraph, ResidualCostMap<CostMap> >
825 bf(*_res_graph, res_cost);
826 bf.init(0); bf.start();
827 for (NodeIt n(_graph); n != INVALID; ++n)
828 (*_potential)[n] = bf.dist(n);
830 // Handle non-zero lower bounds
832 for (EdgeIt e(_graph); e != INVALID; ++e)
833 (*_flow)[e] += (*_lower)[e];
838 }; //class CostScaling
844 #endif //LEMON_COST_SCALING_H