|
1 /* -*- C++ -*- |
|
2 * |
|
3 * This file is a part of LEMON, a generic C++ optimization library |
|
4 * |
|
5 * Copyright (C) 2003-2008 |
|
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
|
7 * (Egervary Research Group on Combinatorial Optimization, EGRES). |
|
8 * |
|
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. |
|
12 * |
|
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 |
|
15 * purpose. |
|
16 * |
|
17 */ |
|
18 |
|
19 #ifndef LEMON_CYCLE_CANCELING_H |
|
20 #define LEMON_CYCLE_CANCELING_H |
|
21 |
|
22 /// \ingroup min_cost_flow |
|
23 /// |
|
24 /// \file |
|
25 /// \brief Cycle-canceling algorithm for finding a minimum cost flow. |
|
26 |
|
27 #include <vector> |
|
28 #include <lemon/adaptors.h> |
|
29 #include <lemon/path.h> |
|
30 |
|
31 #include <lemon/circulation.h> |
|
32 #include <lemon/bellman_ford.h> |
|
33 #include <lemon/howard.h> |
|
34 |
|
35 namespace lemon { |
|
36 |
|
37 /// \addtogroup min_cost_flow |
|
38 /// @{ |
|
39 |
|
40 /// \brief Implementation of a cycle-canceling algorithm for |
|
41 /// finding a minimum cost flow. |
|
42 /// |
|
43 /// \ref CycleCanceling implements a cycle-canceling algorithm for |
|
44 /// finding a minimum cost flow. |
|
45 /// |
|
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. |
|
51 /// |
|
52 /// \warning |
|
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. |
|
57 /// |
|
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 |
|
64 /// parameter. |
|
65 /// |
|
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> > |
|
72 class CycleCanceling |
|
73 { |
|
74 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
|
75 |
|
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; |
|
81 |
|
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; |
|
88 |
|
89 public: |
|
90 |
|
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; |
|
95 |
|
96 private: |
|
97 |
|
98 /// \brief Map adaptor class for handling residual arc costs. |
|
99 /// |
|
100 /// Map adaptor class for handling residual arc costs. |
|
101 class ResidualCostMap : public MapBase<ResArc, Cost> |
|
102 { |
|
103 private: |
|
104 |
|
105 const CostMap &_cost_map; |
|
106 |
|
107 public: |
|
108 |
|
109 ///\e |
|
110 ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {} |
|
111 |
|
112 ///\e |
|
113 Cost operator[](const ResArc &e) const { |
|
114 return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e]; |
|
115 } |
|
116 |
|
117 }; //class ResidualCostMap |
|
118 |
|
119 private: |
|
120 |
|
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; |
|
127 |
|
128 private: |
|
129 |
|
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; |
|
140 bool _valid_supply; |
|
141 |
|
142 // Arc map of the current flow |
|
143 FlowMap *_flow; |
|
144 bool _local_flow; |
|
145 // Node map of the current potentials |
|
146 PotentialMap *_potential; |
|
147 bool _local_potential; |
|
148 |
|
149 // The residual digraph |
|
150 ResDigraph *_res_graph; |
|
151 // The residual cost map |
|
152 ResidualCostMap _res_cost; |
|
153 |
|
154 public: |
|
155 |
|
156 /// \brief General constructor (with lower bounds). |
|
157 /// |
|
158 /// General constructor (with lower bounds). |
|
159 /// |
|
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, |
|
168 const CostMap &cost, |
|
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) |
|
174 { |
|
175 // Check the sum of supply values |
|
176 Supply sum = 0; |
|
177 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
178 _supply[n] = supply[n]; |
|
179 sum += _supply[n]; |
|
180 } |
|
181 _valid_supply = sum == 0; |
|
182 |
|
183 // Remove non-zero lower bounds |
|
184 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
185 _capacity[e] = capacity[e]; |
|
186 if (lower[e] != 0) { |
|
187 _capacity[e] -= lower[e]; |
|
188 _supply[_graph.source(e)] -= lower[e]; |
|
189 _supply[_graph.target(e)] += lower[e]; |
|
190 } |
|
191 } |
|
192 } |
|
193 /* |
|
194 /// \brief General constructor (without lower bounds). |
|
195 /// |
|
196 /// General constructor (without lower bounds). |
|
197 /// |
|
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, |
|
204 const CostMap &cost, |
|
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), |
|
209 _res_cost(_cost) |
|
210 { |
|
211 // Check the sum of supply values |
|
212 Supply sum = 0; |
|
213 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; |
|
214 _valid_supply = sum == 0; |
|
215 } |
|
216 |
|
217 /// \brief Simple constructor (with lower bounds). |
|
218 /// |
|
219 /// Simple constructor (with lower bounds). |
|
220 /// |
|
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, |
|
232 const CostMap &cost, |
|
233 Node s, Node t, |
|
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), |
|
238 _res_cost(_cost) |
|
239 { |
|
240 // Remove non-zero lower bounds |
|
241 _supply[s] = flow_value; |
|
242 _supply[t] = -flow_value; |
|
243 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
244 if (lower[e] != 0) { |
|
245 _capacity[e] -= lower[e]; |
|
246 _supply[_graph.source(e)] -= lower[e]; |
|
247 _supply[_graph.target(e)] += lower[e]; |
|
248 } |
|
249 } |
|
250 _valid_supply = true; |
|
251 } |
|
252 |
|
253 /// \brief Simple constructor (without lower bounds). |
|
254 /// |
|
255 /// Simple constructor (without lower bounds). |
|
256 /// |
|
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, |
|
266 const CostMap &cost, |
|
267 Node s, Node t, |
|
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), |
|
272 _res_cost(_cost) |
|
273 { |
|
274 _supply[s] = flow_value; |
|
275 _supply[t] = -flow_value; |
|
276 _valid_supply = true; |
|
277 } |
|
278 */ |
|
279 /// Destructor. |
|
280 ~CycleCanceling() { |
|
281 if (_local_flow) delete _flow; |
|
282 if (_local_potential) delete _potential; |
|
283 delete _res_graph; |
|
284 } |
|
285 |
|
286 /// \brief Set the flow map. |
|
287 /// |
|
288 /// Set the flow map. |
|
289 /// |
|
290 /// \return \c (*this) |
|
291 CycleCanceling& flowMap(FlowMap &map) { |
|
292 if (_local_flow) { |
|
293 delete _flow; |
|
294 _local_flow = false; |
|
295 } |
|
296 _flow = ↦ |
|
297 return *this; |
|
298 } |
|
299 |
|
300 /// \brief Set the potential map. |
|
301 /// |
|
302 /// Set the potential map. |
|
303 /// |
|
304 /// \return \c (*this) |
|
305 CycleCanceling& potentialMap(PotentialMap &map) { |
|
306 if (_local_potential) { |
|
307 delete _potential; |
|
308 _local_potential = false; |
|
309 } |
|
310 _potential = ↦ |
|
311 return *this; |
|
312 } |
|
313 |
|
314 /// \name Execution control |
|
315 |
|
316 /// @{ |
|
317 |
|
318 /// \brief Run the algorithm. |
|
319 /// |
|
320 /// Run the algorithm. |
|
321 /// |
|
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. |
|
325 /// |
|
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); |
|
329 } |
|
330 |
|
331 /// @} |
|
332 |
|
333 /// \name Query Functions |
|
334 /// The result of the algorithm can be obtained using these |
|
335 /// functions.\n |
|
336 /// \ref lemon::CycleCanceling::run() "run()" must be called before |
|
337 /// using them. |
|
338 |
|
339 /// @{ |
|
340 |
|
341 /// \brief Return a const reference to the arc map storing the |
|
342 /// found flow. |
|
343 /// |
|
344 /// Return a const reference to the arc map storing the found flow. |
|
345 /// |
|
346 /// \pre \ref run() must be called before using this function. |
|
347 const FlowMap& flowMap() const { |
|
348 return *_flow; |
|
349 } |
|
350 |
|
351 /// \brief Return a const reference to the node map storing the |
|
352 /// found potentials (the dual solution). |
|
353 /// |
|
354 /// Return a const reference to the node map storing the found |
|
355 /// potentials (the dual solution). |
|
356 /// |
|
357 /// \pre \ref run() must be called before using this function. |
|
358 const PotentialMap& potentialMap() const { |
|
359 return *_potential; |
|
360 } |
|
361 |
|
362 /// \brief Return the flow on the given arc. |
|
363 /// |
|
364 /// Return the flow on the given arc. |
|
365 /// |
|
366 /// \pre \ref run() must be called before using this function. |
|
367 Capacity flow(const Arc& arc) const { |
|
368 return (*_flow)[arc]; |
|
369 } |
|
370 |
|
371 /// \brief Return the potential of the given node. |
|
372 /// |
|
373 /// Return the potential of the given node. |
|
374 /// |
|
375 /// \pre \ref run() must be called before using this function. |
|
376 Cost potential(const Node& node) const { |
|
377 return (*_potential)[node]; |
|
378 } |
|
379 |
|
380 /// \brief Return the total cost of the found flow. |
|
381 /// |
|
382 /// Return the total cost of the found flow. The complexity of the |
|
383 /// function is \f$ O(e) \f$. |
|
384 /// |
|
385 /// \pre \ref run() must be called before using this function. |
|
386 Cost totalCost() const { |
|
387 Cost c = 0; |
|
388 for (ArcIt e(_graph); e != INVALID; ++e) |
|
389 c += (*_flow)[e] * _cost[e]; |
|
390 return c; |
|
391 } |
|
392 |
|
393 /// @} |
|
394 |
|
395 private: |
|
396 |
|
397 /// Initialize the algorithm. |
|
398 bool init() { |
|
399 if (!_valid_supply) return false; |
|
400 |
|
401 // Initializing flow and potential maps |
|
402 if (!_flow) { |
|
403 _flow = new FlowMap(_graph); |
|
404 _local_flow = true; |
|
405 } |
|
406 if (!_potential) { |
|
407 _potential = new PotentialMap(_graph); |
|
408 _local_potential = true; |
|
409 } |
|
410 |
|
411 _res_graph = new ResDigraph(_graph, _capacity, *_flow); |
|
412 |
|
413 // Finding a feasible flow using Circulation |
|
414 Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap, |
|
415 SupplyMap > |
|
416 circulation( _graph, constMap<Arc>(Capacity(0)), _capacity, |
|
417 _supply ); |
|
418 return circulation.flowMap(*_flow).run(); |
|
419 } |
|
420 |
|
421 bool start(bool min_mean_cc) { |
|
422 if (min_mean_cc) |
|
423 startMinMean(); |
|
424 else |
|
425 start(); |
|
426 |
|
427 // Handling non-zero lower bounds |
|
428 if (_lower) { |
|
429 for (ArcIt e(_graph); e != INVALID; ++e) |
|
430 (*_flow)[e] += (*_lower)[e]; |
|
431 } |
|
432 return true; |
|
433 } |
|
434 |
|
435 /// \brief Execute the algorithm using \ref BellmanFord. |
|
436 /// |
|
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. |
|
440 void start() { |
|
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); |
|
445 |
|
446 int length_bound = BF_FIRST_LIMIT; |
|
447 bool optimal = false; |
|
448 while (!optimal) { |
|
449 BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost); |
|
450 bf.predMap(pred); |
|
451 bf.init(0); |
|
452 int iter_num = 0; |
|
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()) { |
|
461 real_iter_num = i; |
|
462 break; |
|
463 } |
|
464 } |
|
465 if (real_iter_num < curr_iter_num) { |
|
466 // Optimal flow is found |
|
467 optimal = true; |
|
468 // Setting node potentials |
|
469 for (NodeIt n(_graph); n != INVALID; ++n) |
|
470 (*_potential)[n] = bf.dist(n); |
|
471 break; |
|
472 } else { |
|
473 // Searching for node disjoint negative cycles |
|
474 for (ResNodeIt n(*_res_graph); n != INVALID; ++n) |
|
475 visited[n] = 0; |
|
476 int id = 0; |
|
477 for (ResNodeIt n(*_res_graph); n != INVALID; ++n) { |
|
478 if (visited[n] > 0) continue; |
|
479 visited[n] = ++id; |
|
480 ResNode u = pred[n] == INVALID ? |
|
481 INVALID : _res_graph->source(pred[n]); |
|
482 while (u != INVALID && visited[u] == 0) { |
|
483 visited[u] = id; |
|
484 u = pred[u] == INVALID ? |
|
485 INVALID : _res_graph->source(pred[u]); |
|
486 } |
|
487 if (u != INVALID && visited[u] == id) { |
|
488 // Finding the negative cycle |
|
489 cycle_found = true; |
|
490 cycle.clear(); |
|
491 ResArc e = pred[u]; |
|
492 cycle.push_back(e); |
|
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); |
|
498 } |
|
499 |
|
500 // Augmenting along the cycle |
|
501 for (int i = 0; i < int(cycle.size()); ++i) |
|
502 _res_graph->augment(cycle[i], d); |
|
503 } |
|
504 } |
|
505 } |
|
506 |
|
507 if (!cycle_found) |
|
508 length_bound = length_bound * BF_LIMIT_FACTOR / 100; |
|
509 } |
|
510 } |
|
511 } |
|
512 |
|
513 /// \brief Execute the algorithm using \ref Howard. |
|
514 /// |
|
515 /// Execute the algorithm using \ref Howard for negative |
|
516 /// cycle detection. |
|
517 void startMinMean() { |
|
518 typedef Path<ResDigraph> ResPath; |
|
519 Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost); |
|
520 ResPath cycle; |
|
521 |
|
522 mmc.cycle(cycle); |
|
523 if (mmc.findMinMean()) { |
|
524 while (mmc.cycleLength() < 0) { |
|
525 // Finding the cycle |
|
526 mmc.findCycle(); |
|
527 |
|
528 // Finding the largest flow amount that can be augmented |
|
529 // along the cycle |
|
530 Capacity delta = 0; |
|
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); |
|
534 } |
|
535 |
|
536 // Augmenting along the cycle |
|
537 for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) |
|
538 _res_graph->augment(e, delta); |
|
539 |
|
540 // Finding the minimum cycle mean for the modified residual |
|
541 // digraph |
|
542 if (!mmc.findMinMean()) break; |
|
543 } |
|
544 } |
|
545 |
|
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); |
|
551 } |
|
552 |
|
553 }; //class CycleCanceling |
|
554 |
|
555 ///@} |
|
556 |
|
557 } //namespace lemon |
|
558 |
|
559 #endif //LEMON_CYCLE_CANCELING_H |