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_CANCEL_AND_TIGHTEN_H
20 #define LEMON_CANCEL_AND_TIGHTEN_H
22 /// \ingroup min_cost_flow
25 /// \brief Cancel and Tighten algorithm for finding a minimum cost flow.
29 #include <lemon/circulation.h>
30 #include <lemon/bellman_ford.h>
31 #include <lemon/howard.h>
32 #include <lemon/adaptors.h>
33 #include <lemon/tolerance.h>
34 #include <lemon/math.h>
36 #include <lemon/static_graph.h>
40 /// \addtogroup min_cost_flow
43 /// \brief Implementation of the Cancel and Tighten algorithm for
44 /// finding a minimum cost flow.
46 /// \ref CancelAndTighten implements the Cancel and Tighten algorithm for
47 /// finding a minimum cost flow.
49 /// \tparam Digraph The digraph type the algorithm runs on.
50 /// \tparam LowerMap The type of the lower bound map.
51 /// \tparam CapacityMap The type of the capacity (upper bound) map.
52 /// \tparam CostMap The type of the cost (length) map.
53 /// \tparam SupplyMap The type of the supply map.
56 /// - Arc capacities and costs should be \e non-negative \e integers.
57 /// - Supply values should be \e signed \e integers.
58 /// - The value types of the maps should be convertible to each other.
59 /// - \c CostMap::Value must be signed type.
61 /// \author Peter Kovacs
62 template < typename Digraph,
63 typename LowerMap = typename Digraph::template ArcMap<int>,
64 typename CapacityMap = typename Digraph::template ArcMap<int>,
65 typename CostMap = typename Digraph::template ArcMap<int>,
66 typename SupplyMap = typename Digraph::template NodeMap<int> >
67 class CancelAndTighten
69 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
71 typedef typename CapacityMap::Value Capacity;
72 typedef typename CostMap::Value Cost;
73 typedef typename SupplyMap::Value Supply;
74 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
75 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
77 typedef ResidualDigraph< const Digraph,
78 CapacityArcMap, CapacityArcMap > ResDigraph;
82 /// The type of the flow map.
83 typedef typename Digraph::template ArcMap<Capacity> FlowMap;
84 /// The type of the potential map.
85 typedef typename Digraph::template NodeMap<Cost> PotentialMap;
89 /// \brief Map adaptor class for handling residual arc costs.
91 /// Map adaptor class for handling residual arc costs.
92 class ResidualCostMap : public MapBase<typename ResDigraph::Arc, Cost>
94 typedef typename ResDigraph::Arc Arc;
98 const CostMap &_cost_map;
103 ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
106 Cost operator[](const Arc &e) const {
107 return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
110 }; //class ResidualCostMap
112 /// \brief Map adaptor class for handling reduced arc costs.
114 /// Map adaptor class for handling reduced arc costs.
115 class ReducedCostMap : public MapBase<Arc, Cost>
120 const CostMap &_cost_map;
121 const PotentialMap &_pot_map;
126 ReducedCostMap( const Digraph &gr,
127 const CostMap &cost_map,
128 const PotentialMap &pot_map ) :
129 _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
132 inline Cost operator[](const Arc &e) const {
133 return _cost_map[e] + _pot_map[_gr.source(e)]
134 - _pot_map[_gr.target(e)];
137 }; //class ReducedCostMap
139 struct BFOperationTraits {
140 static double zero() { return 0; }
142 static double infinity() {
143 return std::numeric_limits<double>::infinity();
146 static double plus(const double& left, const double& right) {
150 static bool less(const double& left, const double& right) {
151 return left + 1e-6 < right;
153 }; // class BFOperationTraits
157 // The digraph the algorithm runs on
158 const Digraph &_graph;
159 // The original lower bound map
160 const LowerMap *_lower;
161 // The modified capacity map
162 CapacityArcMap _capacity;
163 // The original cost map
164 const CostMap &_cost;
165 // The modified supply map
166 SupplyNodeMap _supply;
169 // Arc map of the current flow
172 // Node map of the current potentials
173 PotentialMap *_potential;
174 bool _local_potential;
176 // The residual digraph
177 ResDigraph *_res_graph;
178 // The residual cost map
179 ResidualCostMap _res_cost;
183 /// \brief General constructor (with lower bounds).
185 /// General constructor (with lower bounds).
187 /// \param digraph The digraph the algorithm runs on.
188 /// \param lower The lower bounds of the arcs.
189 /// \param capacity The capacities (upper bounds) of the arcs.
190 /// \param cost The cost (length) values of the arcs.
191 /// \param supply The supply values of the nodes (signed).
192 CancelAndTighten( const Digraph &digraph,
193 const LowerMap &lower,
194 const CapacityMap &capacity,
196 const SupplyMap &supply ) :
197 _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
198 _supply(digraph), _flow(NULL), _local_flow(false),
199 _potential(NULL), _local_potential(false),
200 _res_graph(NULL), _res_cost(_cost)
202 // Check the sum of supply values
204 for (NodeIt n(_graph); n != INVALID; ++n) {
205 _supply[n] = supply[n];
208 _valid_supply = sum == 0;
210 // Remove non-zero lower bounds
211 for (ArcIt e(_graph); e != INVALID; ++e) {
212 _capacity[e] = capacity[e];
214 _capacity[e] -= lower[e];
215 _supply[_graph.source(e)] -= lower[e];
216 _supply[_graph.target(e)] += lower[e];
221 /// \brief General constructor (without lower bounds).
223 /// General constructor (without lower bounds).
225 /// \param digraph The digraph the algorithm runs on.
226 /// \param capacity The capacities (upper bounds) of the arcs.
227 /// \param cost The cost (length) values of the arcs.
228 /// \param supply The supply values of the nodes (signed).
229 CancelAndTighten( const Digraph &digraph,
230 const CapacityMap &capacity,
232 const SupplyMap &supply ) :
233 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
234 _supply(supply), _flow(NULL), _local_flow(false),
235 _potential(NULL), _local_potential(false),
236 _res_graph(NULL), _res_cost(_cost)
238 // Check the sum of supply values
240 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
241 _valid_supply = sum == 0;
244 /// \brief Simple constructor (with lower bounds).
246 /// Simple constructor (with lower bounds).
248 /// \param digraph The digraph the algorithm runs on.
249 /// \param lower The lower bounds of the arcs.
250 /// \param capacity The capacities (upper bounds) of the arcs.
251 /// \param cost The cost (length) values of the arcs.
252 /// \param s The source node.
253 /// \param t The target node.
254 /// \param flow_value The required amount of flow from node \c s
255 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
256 CancelAndTighten( const Digraph &digraph,
257 const LowerMap &lower,
258 const CapacityMap &capacity,
261 Supply flow_value ) :
262 _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
263 _supply(digraph, 0), _flow(NULL), _local_flow(false),
264 _potential(NULL), _local_potential(false),
265 _res_graph(NULL), _res_cost(_cost)
267 // Remove non-zero lower bounds
268 _supply[s] = flow_value;
269 _supply[t] = -flow_value;
270 for (ArcIt e(_graph); e != INVALID; ++e) {
272 _capacity[e] -= lower[e];
273 _supply[_graph.source(e)] -= lower[e];
274 _supply[_graph.target(e)] += lower[e];
277 _valid_supply = true;
280 /// \brief Simple constructor (without lower bounds).
282 /// Simple constructor (without lower bounds).
284 /// \param digraph The digraph the algorithm runs on.
285 /// \param capacity The capacities (upper bounds) of the arcs.
286 /// \param cost The cost (length) values of the arcs.
287 /// \param s The source node.
288 /// \param t The target node.
289 /// \param flow_value The required amount of flow from node \c s
290 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
291 CancelAndTighten( const Digraph &digraph,
292 const CapacityMap &capacity,
295 Supply flow_value ) :
296 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
297 _supply(digraph, 0), _flow(NULL), _local_flow(false),
298 _potential(NULL), _local_potential(false),
299 _res_graph(NULL), _res_cost(_cost)
301 _supply[s] = flow_value;
302 _supply[t] = -flow_value;
303 _valid_supply = true;
307 ~CancelAndTighten() {
308 if (_local_flow) delete _flow;
309 if (_local_potential) delete _potential;
313 /// \brief Set the flow map.
315 /// Set the flow map.
317 /// \return \c (*this)
318 CancelAndTighten& flowMap(FlowMap &map) {
327 /// \brief Set the potential map.
329 /// Set the potential map.
331 /// \return \c (*this)
332 CancelAndTighten& potentialMap(PotentialMap &map) {
333 if (_local_potential) {
335 _local_potential = false;
341 /// \name Execution control
345 /// \brief Run the algorithm.
347 /// Run the algorithm.
349 /// \return \c true if a feasible flow can be found.
351 return init() && start();
356 /// \name Query Functions
357 /// The result of the algorithm can be obtained using these
359 /// \ref lemon::CancelAndTighten::run() "run()" must be called before
364 /// \brief Return a const reference to the arc map storing the
367 /// Return a const reference to the arc map storing the found flow.
369 /// \pre \ref run() must be called before using this function.
370 const FlowMap& flowMap() const {
374 /// \brief Return a const reference to the node map storing the
375 /// found potentials (the dual solution).
377 /// Return a const reference to the node map storing the found
378 /// potentials (the dual solution).
380 /// \pre \ref run() must be called before using this function.
381 const PotentialMap& potentialMap() const {
385 /// \brief Return the flow on the given arc.
387 /// Return the flow on the given arc.
389 /// \pre \ref run() must be called before using this function.
390 Capacity flow(const Arc& arc) const {
391 return (*_flow)[arc];
394 /// \brief Return the potential of the given node.
396 /// Return the potential of the given node.
398 /// \pre \ref run() must be called before using this function.
399 Cost potential(const Node& node) const {
400 return (*_potential)[node];
403 /// \brief Return the total cost of the found flow.
405 /// Return the total cost of the found flow. The complexity of the
406 /// function is \f$ O(e) \f$.
408 /// \pre \ref run() must be called before using this function.
409 Cost totalCost() const {
411 for (ArcIt e(_graph); e != INVALID; ++e)
412 c += (*_flow)[e] * _cost[e];
420 /// Initialize the algorithm.
422 if (!_valid_supply) return false;
424 // Initialize flow and potential maps
426 _flow = new FlowMap(_graph);
430 _potential = new PotentialMap(_graph);
431 _local_potential = true;
434 _res_graph = new ResDigraph(_graph, _capacity, *_flow);
436 // Find a feasible flow using Circulation
437 Circulation< Digraph, ConstMap<Arc, Capacity>,
438 CapacityArcMap, SupplyMap >
439 circulation( _graph, constMap<Arc>(Capacity(0)),
440 _capacity, _supply );
441 return circulation.flowMap(*_flow).run();
445 const double LIMIT_FACTOR = 0.01;
446 const int MIN_LIMIT = 3;
448 typedef typename Digraph::template NodeMap<double> FloatPotentialMap;
449 typedef typename Digraph::template NodeMap<int> LevelMap;
450 typedef typename Digraph::template NodeMap<bool> BoolNodeMap;
451 typedef typename Digraph::template NodeMap<Node> PredNodeMap;
452 typedef typename Digraph::template NodeMap<Arc> PredArcMap;
453 typedef typename ResDigraph::template ArcMap<double> ResShiftCostMap;
454 FloatPotentialMap pi(_graph);
455 LevelMap level(_graph);
456 BoolNodeMap reached(_graph);
457 BoolNodeMap processed(_graph);
458 PredNodeMap pred_node(_graph);
459 PredArcMap pred_arc(_graph);
460 int node_num = countNodes(_graph);
461 typedef std::pair<Arc, bool> pair;
462 std::vector<pair> stack(node_num);
463 std::vector<Node> proc_vector(node_num);
464 ResShiftCostMap shift_cost(*_res_graph);
466 Tolerance<double> tol;
474 // Initialize epsilon and the node potentials
476 for (ArcIt e(_graph); e != INVALID; ++e) {
477 if (_capacity[e] - (*_flow)[e] > 0 && _cost[e] < -epsilon)
479 else if ((*_flow)[e] > 0 && _cost[e] > epsilon)
482 for (NodeIt v(_graph); v != INVALID; ++v) {
487 int limit = int(LIMIT_FACTOR * node_num);
488 if (limit < MIN_LIMIT) limit = MIN_LIMIT;
490 while (epsilon * node_num >= 1) {
492 // Find and cancel cycles in the admissible digraph using DFS
493 for (NodeIt n(_graph); n != INVALID; ++n) {
495 processed[n] = false;
500 for (NodeIt start(_graph); start != INVALID; ++start) {
501 if (reached[start]) continue;
504 reached[start] = true;
505 pred_arc[start] = INVALID;
506 pred_node[start] = INVALID;
508 // Find the first admissible residual outgoing arc
509 double p = pi[start];
511 _graph.firstOut(e, start);
512 while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
513 !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
516 stack[++stack_head] = pair(e, true);
519 _graph.firstIn(e, start);
520 while ( e != INVALID && ((*_flow)[e] == 0 ||
521 !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
524 stack[++stack_head] = pair(e, false);
527 processed[start] = true;
528 proc_vector[++proc_head] = start;
532 while (stack_head >= 0) {
533 Arc se = stack[stack_head].first;
534 bool sf = stack[stack_head].second;
537 u = _graph.source(se);
538 v = _graph.target(se);
540 u = _graph.target(se);
541 v = _graph.source(se);
545 // A new node is reached
549 // Find the first admissible residual outgoing arc
552 _graph.firstOut(e, v);
553 while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
554 !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
557 stack[++stack_head] = pair(e, true);
560 _graph.firstIn(e, v);
561 while ( e != INVALID && ((*_flow)[e] == 0 ||
562 !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
564 stack[++stack_head] = pair(e, false);
570 Capacity d, delta = sf ? _capacity[se] - (*_flow)[se] :
572 for (n = u; n != v; n = pred_node[n]) {
573 d = _graph.target(pred_arc[n]) == n ?
574 _capacity[pred_arc[n]] - (*_flow)[pred_arc[n]] :
575 (*_flow)[pred_arc[n]];
583 std::cout << "CYCLE FOUND: ";
585 std::cout << _cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)];
587 std::cout << _graph.id(se) << ":" << -(_cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)]);
588 for (n = u; n != v; n = pred_node[n]) {
589 if (_graph.target(pred_arc[n]) == n)
590 std::cout << " " << _cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])];
592 std::cout << " " << -(_cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])]);
596 // Augment along the cycle
597 (*_flow)[se] = sf ? (*_flow)[se] + delta :
598 (*_flow)[se] - delta;
599 for (n = u; n != v; n = pred_node[n]) {
600 if (_graph.target(pred_arc[n]) == n)
601 (*_flow)[pred_arc[n]] += delta;
603 (*_flow)[pred_arc[n]] -= delta;
605 for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
613 // Find the next admissible residual outgoing arc
615 Arc e = stack[stack_head].first;
616 if (!stack[stack_head].second) {
621 while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
622 !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
625 stack[stack_head] = pair(e, true);
628 _graph.firstIn(e, v);
630 while ( e != INVALID && ((*_flow)[e] == 0 ||
631 !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
633 stack[stack_head] = pair(e, false);
637 while (stack_head >= 0 && stack[stack_head].first == INVALID) {
639 proc_vector[++proc_head] = v;
640 if (--stack_head >= 0) {
641 v = stack[stack_head].second ?
642 _graph.source(stack[stack_head].first) :
643 _graph.target(stack[stack_head].first);
644 // Find the next admissible residual outgoing arc
646 Arc e = stack[stack_head].first;
647 if (!stack[stack_head].second) {
652 while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
653 !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
656 stack[stack_head] = pair(e, true);
659 _graph.firstIn(e, v);
661 while ( e != INVALID && ((*_flow)[e] == 0 ||
662 !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
664 stack[stack_head] = pair(e, false);
672 // Tighten potentials and epsilon
676 for (int i = proc_head; i >= 0; --i) {
677 Node v = proc_vector[i];
680 for (InArcIt e(_graph, v); e != INVALID; ++e) {
681 Node u = _graph.source(e);
682 if ( _capacity[e] - (*_flow)[e] > 0 &&
683 tol.negative(_cost[e] + pi[u] - p) &&
684 level[u] + 1 > l ) l = level[u] + 1;
686 for (OutArcIt e(_graph, v); e != INVALID; ++e) {
687 Node u = _graph.target(e);
688 if ( (*_flow)[e] > 0 &&
689 tol.negative(-_cost[e] + pi[u] - p) &&
690 level[u] + 1 > l ) l = level[u] + 1;
697 for (ArcIt e(_graph); e != INVALID; ++e) {
698 Node u = _graph.source(e);
699 Node v = _graph.target(e);
700 if (_capacity[e] - (*_flow)[e] > 0 && level[u] - level[v] > 0) {
701 p = (_cost[e] + pi[u] - pi[v] + epsilon) /
702 (level[u] - level[v] + 1);
703 if (q < 0 || p < q) q = p;
705 else if ((*_flow)[e] > 0 && level[v] - level[u] > 0) {
706 p = (-_cost[e] - pi[u] + pi[v] + epsilon) /
707 (level[v] - level[u] + 1);
708 if (q < 0 || p < q) q = p;
711 for (NodeIt v(_graph); v != INVALID; ++v) {
712 pi[v] -= q * level[v];
717 for (ArcIt e(_graph); e != INVALID; ++e) {
718 double curr = _cost[e] + pi[_graph.source(e)]
719 - pi[_graph.target(e)];
720 if (_capacity[e] - (*_flow)[e] > 0 && curr < -epsilon)
722 else if ((*_flow)[e] > 0 && curr > epsilon)
727 // Set epsilon to the minimum cycle mean
731 StaticDigraph static_graph;
732 typename ResDigraph::template NodeMap<typename StaticDigraph::Node> node_ref(*_res_graph);
733 typename ResDigraph::template ArcMap<typename StaticDigraph::Arc> arc_ref(*_res_graph);
734 static_graph.build(*_res_graph, node_ref, arc_ref);
735 typename StaticDigraph::template NodeMap<double> static_pi(static_graph);
736 typename StaticDigraph::template ArcMap<double> static_cost(static_graph);
738 for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
739 static_cost[arc_ref[e]] = _res_cost[e];
741 Howard<StaticDigraph, typename StaticDigraph::template ArcMap<double> >
742 mmc(static_graph, static_cost);
744 epsilon = -mmc.cycleMean();
748 Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
750 epsilon = -mmc.cycleMean();
753 // Compute feasible potentials for the current epsilon
754 for (typename StaticDigraph::ArcIt e(static_graph); e != INVALID; ++e)
755 static_cost[e] += epsilon;
756 typename BellmanFord<StaticDigraph, typename StaticDigraph::template ArcMap<double> >::
757 template SetDistMap<typename StaticDigraph::template NodeMap<double> >::
758 template SetOperationTraits<BFOperationTraits>::Create
759 bf(static_graph, static_cost);
760 bf.distMap(static_pi).init(0);
762 for (NodeIt n(_graph); n != INVALID; ++n)
763 pi[n] = static_pi[node_ref[n]];
766 for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
767 shift_cost[e] = _res_cost[e] + epsilon;
768 typename BellmanFord<ResDigraph, ResShiftCostMap>::
769 template SetDistMap<FloatPotentialMap>::
770 template SetOperationTraits<BFOperationTraits>::Create
771 bf(*_res_graph, shift_cost);
772 bf.distMap(pi).init(0);
781 // std::cout << t1.realTime() << " " << t2.realTime() << " " << t3.realTime() << "\n";
783 // Handle non-zero lower bounds
785 for (ArcIt e(_graph); e != INVALID; ++e)
786 (*_flow)[e] += (*_lower)[e];
791 }; //class CancelAndTighten
797 #endif //LEMON_CANCEL_AND_TIGHTEN_H