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_CYCLE_CANCELING_H
20 #define LEMON_CYCLE_CANCELING_H
22 /// \ingroup min_cost_flow
25 /// \brief Cycle-canceling algorithm for finding a minimum cost flow.
28 #include <lemon/adaptors.h>
29 #include <lemon/path.h>
31 #include <lemon/circulation.h>
32 #include <lemon/bellman_ford.h>
33 #include <lemon/howard.h>
37 /// \addtogroup min_cost_flow
40 /// \brief Implementation of a cycle-canceling algorithm for
41 /// finding a minimum cost flow.
43 /// \ref CycleCanceling implements a cycle-canceling algorithm for
44 /// finding a minimum cost flow.
46 /// \tparam Digraph The digraph type the algorithm runs on.
47 /// \tparam LowerMap The type of the lower bound map.
48 /// \tparam CapacityMap The type of the capacity (upper bound) map.
49 /// \tparam CostMap The type of the cost (length) map.
50 /// \tparam SupplyMap The type of the supply map.
53 /// - Arc capacities and costs should be \e non-negative \e integers.
54 /// - Supply values should be \e signed \e integers.
55 /// - The value types of the maps should be convertible to each other.
56 /// - \c CostMap::Value must be signed type.
58 /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
59 /// used for negative cycle detection with limited iteration number.
60 /// However \ref CycleCanceling also provides the "Minimum Mean
61 /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial,
62 /// but rather slower in practice.
63 /// To use this version of the algorithm, call \ref run() with \c true
66 /// \author Peter Kovacs
67 template < typename Digraph,
68 typename LowerMap = typename Digraph::template ArcMap<int>,
69 typename CapacityMap = typename Digraph::template ArcMap<int>,
70 typename CostMap = typename Digraph::template ArcMap<int>,
71 typename SupplyMap = typename Digraph::template NodeMap<int> >
74 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
76 typedef typename CapacityMap::Value Capacity;
77 typedef typename CostMap::Value Cost;
78 typedef typename SupplyMap::Value Supply;
79 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
80 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
82 typedef ResidualDigraph< const Digraph,
83 CapacityArcMap, CapacityArcMap > ResDigraph;
84 typedef typename ResDigraph::Node ResNode;
85 typedef typename ResDigraph::NodeIt ResNodeIt;
86 typedef typename ResDigraph::Arc ResArc;
87 typedef typename ResDigraph::ArcIt ResArcIt;
91 /// The type of the flow map.
92 typedef typename Digraph::template ArcMap<Capacity> FlowMap;
93 /// The type of the potential map.
94 typedef typename Digraph::template NodeMap<Cost> PotentialMap;
98 /// \brief Map adaptor class for handling residual arc costs.
100 /// Map adaptor class for handling residual arc costs.
101 class ResidualCostMap : public MapBase<ResArc, Cost>
105 const CostMap &_cost_map;
110 ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
113 Cost operator[](const ResArc &e) const {
114 return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
117 }; //class ResidualCostMap
121 // The maximum number of iterations for the first execution of the
122 // Bellman-Ford algorithm. It should be at least 2.
123 static const int BF_FIRST_LIMIT = 2;
124 // The iteration limit for the Bellman-Ford algorithm is multiplied
125 // by BF_LIMIT_FACTOR/100 in every round.
126 static const int BF_LIMIT_FACTOR = 150;
130 // The digraph the algorithm runs on
131 const Digraph &_graph;
132 // The original lower bound map
133 const LowerMap *_lower;
134 // The modified capacity map
135 CapacityArcMap _capacity;
136 // The original cost map
137 const CostMap &_cost;
138 // The modified supply map
139 SupplyNodeMap _supply;
142 // Arc map of the current flow
145 // Node map of the current potentials
146 PotentialMap *_potential;
147 bool _local_potential;
149 // The residual digraph
150 ResDigraph *_res_graph;
151 // The residual cost map
152 ResidualCostMap _res_cost;
156 /// \brief General constructor (with lower bounds).
158 /// General constructor (with lower bounds).
160 /// \param digraph The digraph the algorithm runs on.
161 /// \param lower The lower bounds of the arcs.
162 /// \param capacity The capacities (upper bounds) of the arcs.
163 /// \param cost The cost (length) values of the arcs.
164 /// \param supply The supply values of the nodes (signed).
165 CycleCanceling( const Digraph &digraph,
166 const LowerMap &lower,
167 const CapacityMap &capacity,
169 const SupplyMap &supply ) :
170 _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
171 _supply(digraph), _flow(NULL), _local_flow(false),
172 _potential(NULL), _local_potential(false),
173 _res_graph(NULL), _res_cost(_cost)
175 // Check the sum of supply values
177 for (NodeIt n(_graph); n != INVALID; ++n) {
178 _supply[n] = supply[n];
181 _valid_supply = sum == 0;
183 // Remove non-zero lower bounds
184 for (ArcIt e(_graph); e != INVALID; ++e) {
185 _capacity[e] = capacity[e];
187 _capacity[e] -= lower[e];
188 _supply[_graph.source(e)] -= lower[e];
189 _supply[_graph.target(e)] += lower[e];
194 /// \brief General constructor (without lower bounds).
196 /// General constructor (without lower bounds).
198 /// \param digraph The digraph the algorithm runs on.
199 /// \param capacity The capacities (upper bounds) of the arcs.
200 /// \param cost The cost (length) values of the arcs.
201 /// \param supply The supply values of the nodes (signed).
202 CycleCanceling( const Digraph &digraph,
203 const CapacityMap &capacity,
205 const SupplyMap &supply ) :
206 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
207 _supply(supply), _flow(NULL), _local_flow(false),
208 _potential(NULL), _local_potential(false), _res_graph(NULL),
211 // Check the sum of supply values
213 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
214 _valid_supply = sum == 0;
217 /// \brief Simple constructor (with lower bounds).
219 /// Simple constructor (with lower bounds).
221 /// \param digraph The digraph the algorithm runs on.
222 /// \param lower The lower bounds of the arcs.
223 /// \param capacity The capacities (upper bounds) of the arcs.
224 /// \param cost The cost (length) values of the arcs.
225 /// \param s The source node.
226 /// \param t The target node.
227 /// \param flow_value The required amount of flow from node \c s
228 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
229 CycleCanceling( const Digraph &digraph,
230 const LowerMap &lower,
231 const CapacityMap &capacity,
234 Supply flow_value ) :
235 _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
236 _supply(digraph, 0), _flow(NULL), _local_flow(false),
237 _potential(NULL), _local_potential(false), _res_graph(NULL),
240 // Remove non-zero lower bounds
241 _supply[s] = flow_value;
242 _supply[t] = -flow_value;
243 for (ArcIt e(_graph); e != INVALID; ++e) {
245 _capacity[e] -= lower[e];
246 _supply[_graph.source(e)] -= lower[e];
247 _supply[_graph.target(e)] += lower[e];
250 _valid_supply = true;
253 /// \brief Simple constructor (without lower bounds).
255 /// Simple constructor (without lower bounds).
257 /// \param digraph The digraph the algorithm runs on.
258 /// \param capacity The capacities (upper bounds) of the arcs.
259 /// \param cost The cost (length) values of the arcs.
260 /// \param s The source node.
261 /// \param t The target node.
262 /// \param flow_value The required amount of flow from node \c s
263 /// to node \c t (i.e. the supply of \c s and the demand of \c t).
264 CycleCanceling( const Digraph &digraph,
265 const CapacityMap &capacity,
268 Supply flow_value ) :
269 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
270 _supply(digraph, 0), _flow(NULL), _local_flow(false),
271 _potential(NULL), _local_potential(false), _res_graph(NULL),
274 _supply[s] = flow_value;
275 _supply[t] = -flow_value;
276 _valid_supply = true;
281 if (_local_flow) delete _flow;
282 if (_local_potential) delete _potential;
286 /// \brief Set the flow map.
288 /// Set the flow map.
290 /// \return \c (*this)
291 CycleCanceling& flowMap(FlowMap &map) {
300 /// \brief Set the potential map.
302 /// Set the potential map.
304 /// \return \c (*this)
305 CycleCanceling& potentialMap(PotentialMap &map) {
306 if (_local_potential) {
308 _local_potential = false;
314 /// \name Execution control
318 /// \brief Run the algorithm.
320 /// Run the algorithm.
322 /// \param min_mean_cc Set this parameter to \c true to run the
323 /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly
324 /// polynomial, but rather slower in practice.
326 /// \return \c true if a feasible flow can be found.
327 bool run(bool min_mean_cc = false) {
328 return init() && start(min_mean_cc);
333 /// \name Query Functions
334 /// The result of the algorithm can be obtained using these
336 /// \ref lemon::CycleCanceling::run() "run()" must be called before
341 /// \brief Return a const reference to the arc map storing the
344 /// Return a const reference to the arc map storing the found flow.
346 /// \pre \ref run() must be called before using this function.
347 const FlowMap& flowMap() const {
351 /// \brief Return a const reference to the node map storing the
352 /// found potentials (the dual solution).
354 /// Return a const reference to the node map storing the found
355 /// potentials (the dual solution).
357 /// \pre \ref run() must be called before using this function.
358 const PotentialMap& potentialMap() const {
362 /// \brief Return the flow on the given arc.
364 /// Return the flow on the given arc.
366 /// \pre \ref run() must be called before using this function.
367 Capacity flow(const Arc& arc) const {
368 return (*_flow)[arc];
371 /// \brief Return the potential of the given node.
373 /// Return the potential of the given node.
375 /// \pre \ref run() must be called before using this function.
376 Cost potential(const Node& node) const {
377 return (*_potential)[node];
380 /// \brief Return the total cost of the found flow.
382 /// Return the total cost of the found flow. The complexity of the
383 /// function is \f$ O(e) \f$.
385 /// \pre \ref run() must be called before using this function.
386 Cost totalCost() const {
388 for (ArcIt e(_graph); e != INVALID; ++e)
389 c += (*_flow)[e] * _cost[e];
397 /// Initialize the algorithm.
399 if (!_valid_supply) return false;
401 // Initializing flow and potential maps
403 _flow = new FlowMap(_graph);
407 _potential = new PotentialMap(_graph);
408 _local_potential = true;
411 _res_graph = new ResDigraph(_graph, _capacity, *_flow);
413 // Finding a feasible flow using Circulation
414 Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
416 circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
418 return circulation.flowMap(*_flow).run();
421 bool start(bool min_mean_cc) {
427 // Handling non-zero lower bounds
429 for (ArcIt e(_graph); e != INVALID; ++e)
430 (*_flow)[e] += (*_lower)[e];
435 /// \brief Execute the algorithm using \ref BellmanFord.
437 /// Execute the algorithm using the \ref BellmanFord
438 /// "Bellman-Ford" algorithm for negative cycle detection with
439 /// successively larger limit for the number of iterations.
441 typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph);
442 typename ResDigraph::template NodeMap<int> visited(*_res_graph);
443 std::vector<ResArc> cycle;
444 int node_num = countNodes(_graph);
446 int length_bound = BF_FIRST_LIMIT;
447 bool optimal = false;
449 BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
453 bool cycle_found = false;
454 while (!cycle_found) {
455 int curr_iter_num = iter_num + length_bound <= node_num ?
456 length_bound : node_num - iter_num;
457 iter_num += curr_iter_num;
458 int real_iter_num = curr_iter_num;
459 for (int i = 0; i < curr_iter_num; ++i) {
460 if (bf.processNextWeakRound()) {
465 if (real_iter_num < curr_iter_num) {
466 // Optimal flow is found
468 // Setting node potentials
469 for (NodeIt n(_graph); n != INVALID; ++n)
470 (*_potential)[n] = bf.dist(n);
473 // Searching for node disjoint negative cycles
474 for (ResNodeIt n(*_res_graph); n != INVALID; ++n)
477 for (ResNodeIt n(*_res_graph); n != INVALID; ++n) {
478 if (visited[n] > 0) continue;
480 ResNode u = pred[n] == INVALID ?
481 INVALID : _res_graph->source(pred[n]);
482 while (u != INVALID && visited[u] == 0) {
484 u = pred[u] == INVALID ?
485 INVALID : _res_graph->source(pred[u]);
487 if (u != INVALID && visited[u] == id) {
488 // Finding the negative cycle
493 Capacity d = _res_graph->residualCapacity(e);
494 while (_res_graph->source(e) != u) {
495 cycle.push_back(e = pred[_res_graph->source(e)]);
496 if (_res_graph->residualCapacity(e) < d)
497 d = _res_graph->residualCapacity(e);
500 // Augmenting along the cycle
501 for (int i = 0; i < int(cycle.size()); ++i)
502 _res_graph->augment(cycle[i], d);
508 length_bound = length_bound * BF_LIMIT_FACTOR / 100;
513 /// \brief Execute the algorithm using \ref Howard.
515 /// Execute the algorithm using \ref Howard for negative
517 void startMinMean() {
518 typedef Path<ResDigraph> ResPath;
519 Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
523 if (mmc.findMinMean()) {
524 while (mmc.cycleLength() < 0) {
528 // Finding the largest flow amount that can be augmented
531 for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) {
532 if (delta == 0 || _res_graph->residualCapacity(e) < delta)
533 delta = _res_graph->residualCapacity(e);
536 // Augmenting along the cycle
537 for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e)
538 _res_graph->augment(e, delta);
540 // Finding the minimum cycle mean for the modified residual
542 if (!mmc.findMinMean()) break;
546 // Computing node potentials
547 BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
548 bf.init(0); bf.start();
549 for (NodeIt n(_graph); n != INVALID; ++n)
550 (*_potential)[n] = bf.dist(n);
553 }; //class CycleCanceling
559 #endif //LEMON_CYCLE_CANCELING_H