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_CANCEL_AND_TIGHTEN_H |
|
20 #define LEMON_CANCEL_AND_TIGHTEN_H |
|
21 |
|
22 /// \ingroup min_cost_flow |
|
23 /// |
|
24 /// \file |
|
25 /// \brief Cancel and Tighten algorithm for finding a minimum cost flow. |
|
26 |
|
27 #include <vector> |
|
28 |
|
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> |
|
35 |
|
36 #include <lemon/static_graph.h> |
|
37 |
|
38 namespace lemon { |
|
39 |
|
40 /// \addtogroup min_cost_flow |
|
41 /// @{ |
|
42 |
|
43 /// \brief Implementation of the Cancel and Tighten algorithm for |
|
44 /// finding a minimum cost flow. |
|
45 /// |
|
46 /// \ref CancelAndTighten implements the Cancel and Tighten algorithm for |
|
47 /// finding a minimum cost flow. |
|
48 /// |
|
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. |
|
54 /// |
|
55 /// \warning |
|
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. |
|
60 /// |
|
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 |
|
68 { |
|
69 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
|
70 |
|
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; |
|
76 |
|
77 typedef ResidualDigraph< const Digraph, |
|
78 CapacityArcMap, CapacityArcMap > ResDigraph; |
|
79 |
|
80 public: |
|
81 |
|
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; |
|
86 |
|
87 private: |
|
88 |
|
89 /// \brief Map adaptor class for handling residual arc costs. |
|
90 /// |
|
91 /// Map adaptor class for handling residual arc costs. |
|
92 class ResidualCostMap : public MapBase<typename ResDigraph::Arc, Cost> |
|
93 { |
|
94 typedef typename ResDigraph::Arc Arc; |
|
95 |
|
96 private: |
|
97 |
|
98 const CostMap &_cost_map; |
|
99 |
|
100 public: |
|
101 |
|
102 ///\e |
|
103 ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {} |
|
104 |
|
105 ///\e |
|
106 Cost operator[](const Arc &e) const { |
|
107 return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e]; |
|
108 } |
|
109 |
|
110 }; //class ResidualCostMap |
|
111 |
|
112 /// \brief Map adaptor class for handling reduced arc costs. |
|
113 /// |
|
114 /// Map adaptor class for handling reduced arc costs. |
|
115 class ReducedCostMap : public MapBase<Arc, Cost> |
|
116 { |
|
117 private: |
|
118 |
|
119 const Digraph &_gr; |
|
120 const CostMap &_cost_map; |
|
121 const PotentialMap &_pot_map; |
|
122 |
|
123 public: |
|
124 |
|
125 ///\e |
|
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) {} |
|
130 |
|
131 ///\e |
|
132 inline Cost operator[](const Arc &e) const { |
|
133 return _cost_map[e] + _pot_map[_gr.source(e)] |
|
134 - _pot_map[_gr.target(e)]; |
|
135 } |
|
136 |
|
137 }; //class ReducedCostMap |
|
138 |
|
139 struct BFOperationTraits { |
|
140 static double zero() { return 0; } |
|
141 |
|
142 static double infinity() { |
|
143 return std::numeric_limits<double>::infinity(); |
|
144 } |
|
145 |
|
146 static double plus(const double& left, const double& right) { |
|
147 return left + right; |
|
148 } |
|
149 |
|
150 static bool less(const double& left, const double& right) { |
|
151 return left + 1e-6 < right; |
|
152 } |
|
153 }; // class BFOperationTraits |
|
154 |
|
155 private: |
|
156 |
|
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; |
|
167 bool _valid_supply; |
|
168 |
|
169 // Arc map of the current flow |
|
170 FlowMap *_flow; |
|
171 bool _local_flow; |
|
172 // Node map of the current potentials |
|
173 PotentialMap *_potential; |
|
174 bool _local_potential; |
|
175 |
|
176 // The residual digraph |
|
177 ResDigraph *_res_graph; |
|
178 // The residual cost map |
|
179 ResidualCostMap _res_cost; |
|
180 |
|
181 public: |
|
182 |
|
183 /// \brief General constructor (with lower bounds). |
|
184 /// |
|
185 /// General constructor (with lower bounds). |
|
186 /// |
|
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, |
|
195 const CostMap &cost, |
|
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) |
|
201 { |
|
202 // Check the sum of supply values |
|
203 Supply sum = 0; |
|
204 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
205 _supply[n] = supply[n]; |
|
206 sum += _supply[n]; |
|
207 } |
|
208 _valid_supply = sum == 0; |
|
209 |
|
210 // Remove non-zero lower bounds |
|
211 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
212 _capacity[e] = capacity[e]; |
|
213 if (lower[e] != 0) { |
|
214 _capacity[e] -= lower[e]; |
|
215 _supply[_graph.source(e)] -= lower[e]; |
|
216 _supply[_graph.target(e)] += lower[e]; |
|
217 } |
|
218 } |
|
219 } |
|
220 /* |
|
221 /// \brief General constructor (without lower bounds). |
|
222 /// |
|
223 /// General constructor (without lower bounds). |
|
224 /// |
|
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, |
|
231 const CostMap &cost, |
|
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) |
|
237 { |
|
238 // Check the sum of supply values |
|
239 Supply sum = 0; |
|
240 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; |
|
241 _valid_supply = sum == 0; |
|
242 } |
|
243 |
|
244 /// \brief Simple constructor (with lower bounds). |
|
245 /// |
|
246 /// Simple constructor (with lower bounds). |
|
247 /// |
|
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, |
|
259 const CostMap &cost, |
|
260 Node s, Node t, |
|
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) |
|
266 { |
|
267 // Remove non-zero lower bounds |
|
268 _supply[s] = flow_value; |
|
269 _supply[t] = -flow_value; |
|
270 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
271 if (lower[e] != 0) { |
|
272 _capacity[e] -= lower[e]; |
|
273 _supply[_graph.source(e)] -= lower[e]; |
|
274 _supply[_graph.target(e)] += lower[e]; |
|
275 } |
|
276 } |
|
277 _valid_supply = true; |
|
278 } |
|
279 |
|
280 /// \brief Simple constructor (without lower bounds). |
|
281 /// |
|
282 /// Simple constructor (without lower bounds). |
|
283 /// |
|
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, |
|
293 const CostMap &cost, |
|
294 Node s, Node t, |
|
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) |
|
300 { |
|
301 _supply[s] = flow_value; |
|
302 _supply[t] = -flow_value; |
|
303 _valid_supply = true; |
|
304 } |
|
305 */ |
|
306 /// Destructor. |
|
307 ~CancelAndTighten() { |
|
308 if (_local_flow) delete _flow; |
|
309 if (_local_potential) delete _potential; |
|
310 delete _res_graph; |
|
311 } |
|
312 |
|
313 /// \brief Set the flow map. |
|
314 /// |
|
315 /// Set the flow map. |
|
316 /// |
|
317 /// \return \c (*this) |
|
318 CancelAndTighten& flowMap(FlowMap &map) { |
|
319 if (_local_flow) { |
|
320 delete _flow; |
|
321 _local_flow = false; |
|
322 } |
|
323 _flow = ↦ |
|
324 return *this; |
|
325 } |
|
326 |
|
327 /// \brief Set the potential map. |
|
328 /// |
|
329 /// Set the potential map. |
|
330 /// |
|
331 /// \return \c (*this) |
|
332 CancelAndTighten& potentialMap(PotentialMap &map) { |
|
333 if (_local_potential) { |
|
334 delete _potential; |
|
335 _local_potential = false; |
|
336 } |
|
337 _potential = ↦ |
|
338 return *this; |
|
339 } |
|
340 |
|
341 /// \name Execution control |
|
342 |
|
343 /// @{ |
|
344 |
|
345 /// \brief Run the algorithm. |
|
346 /// |
|
347 /// Run the algorithm. |
|
348 /// |
|
349 /// \return \c true if a feasible flow can be found. |
|
350 bool run() { |
|
351 return init() && start(); |
|
352 } |
|
353 |
|
354 /// @} |
|
355 |
|
356 /// \name Query Functions |
|
357 /// The result of the algorithm can be obtained using these |
|
358 /// functions.\n |
|
359 /// \ref lemon::CancelAndTighten::run() "run()" must be called before |
|
360 /// using them. |
|
361 |
|
362 /// @{ |
|
363 |
|
364 /// \brief Return a const reference to the arc map storing the |
|
365 /// found flow. |
|
366 /// |
|
367 /// Return a const reference to the arc map storing the found flow. |
|
368 /// |
|
369 /// \pre \ref run() must be called before using this function. |
|
370 const FlowMap& flowMap() const { |
|
371 return *_flow; |
|
372 } |
|
373 |
|
374 /// \brief Return a const reference to the node map storing the |
|
375 /// found potentials (the dual solution). |
|
376 /// |
|
377 /// Return a const reference to the node map storing the found |
|
378 /// potentials (the dual solution). |
|
379 /// |
|
380 /// \pre \ref run() must be called before using this function. |
|
381 const PotentialMap& potentialMap() const { |
|
382 return *_potential; |
|
383 } |
|
384 |
|
385 /// \brief Return the flow on the given arc. |
|
386 /// |
|
387 /// Return the flow on the given arc. |
|
388 /// |
|
389 /// \pre \ref run() must be called before using this function. |
|
390 Capacity flow(const Arc& arc) const { |
|
391 return (*_flow)[arc]; |
|
392 } |
|
393 |
|
394 /// \brief Return the potential of the given node. |
|
395 /// |
|
396 /// Return the potential of the given node. |
|
397 /// |
|
398 /// \pre \ref run() must be called before using this function. |
|
399 Cost potential(const Node& node) const { |
|
400 return (*_potential)[node]; |
|
401 } |
|
402 |
|
403 /// \brief Return the total cost of the found flow. |
|
404 /// |
|
405 /// Return the total cost of the found flow. The complexity of the |
|
406 /// function is \f$ O(e) \f$. |
|
407 /// |
|
408 /// \pre \ref run() must be called before using this function. |
|
409 Cost totalCost() const { |
|
410 Cost c = 0; |
|
411 for (ArcIt e(_graph); e != INVALID; ++e) |
|
412 c += (*_flow)[e] * _cost[e]; |
|
413 return c; |
|
414 } |
|
415 |
|
416 /// @} |
|
417 |
|
418 private: |
|
419 |
|
420 /// Initialize the algorithm. |
|
421 bool init() { |
|
422 if (!_valid_supply) return false; |
|
423 |
|
424 // Initialize flow and potential maps |
|
425 if (!_flow) { |
|
426 _flow = new FlowMap(_graph); |
|
427 _local_flow = true; |
|
428 } |
|
429 if (!_potential) { |
|
430 _potential = new PotentialMap(_graph); |
|
431 _local_potential = true; |
|
432 } |
|
433 |
|
434 _res_graph = new ResDigraph(_graph, _capacity, *_flow); |
|
435 |
|
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(); |
|
442 } |
|
443 |
|
444 bool start() { |
|
445 const double LIMIT_FACTOR = 0.01; |
|
446 const int MIN_LIMIT = 3; |
|
447 |
|
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); |
|
465 |
|
466 Tolerance<double> tol; |
|
467 tol.epsilon(1e-6); |
|
468 |
|
469 Timer t1, t2, t3; |
|
470 t1.reset(); |
|
471 t2.reset(); |
|
472 t3.reset(); |
|
473 |
|
474 // Initialize epsilon and the node potentials |
|
475 double epsilon = 0; |
|
476 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
477 if (_capacity[e] - (*_flow)[e] > 0 && _cost[e] < -epsilon) |
|
478 epsilon = -_cost[e]; |
|
479 else if ((*_flow)[e] > 0 && _cost[e] > epsilon) |
|
480 epsilon = _cost[e]; |
|
481 } |
|
482 for (NodeIt v(_graph); v != INVALID; ++v) { |
|
483 pi[v] = 0; |
|
484 } |
|
485 |
|
486 // Start phases |
|
487 int limit = int(LIMIT_FACTOR * node_num); |
|
488 if (limit < MIN_LIMIT) limit = MIN_LIMIT; |
|
489 int iter = limit; |
|
490 while (epsilon * node_num >= 1) { |
|
491 t1.start(); |
|
492 // Find and cancel cycles in the admissible digraph using DFS |
|
493 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
494 reached[n] = false; |
|
495 processed[n] = false; |
|
496 } |
|
497 int stack_head = -1; |
|
498 int proc_head = -1; |
|
499 |
|
500 for (NodeIt start(_graph); start != INVALID; ++start) { |
|
501 if (reached[start]) continue; |
|
502 |
|
503 // New start node |
|
504 reached[start] = true; |
|
505 pred_arc[start] = INVALID; |
|
506 pred_node[start] = INVALID; |
|
507 |
|
508 // Find the first admissible residual outgoing arc |
|
509 double p = pi[start]; |
|
510 Arc e; |
|
511 _graph.firstOut(e, start); |
|
512 while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 || |
|
513 !tol.negative(_cost[e] + p - pi[_graph.target(e)])) ) |
|
514 _graph.nextOut(e); |
|
515 if (e != INVALID) { |
|
516 stack[++stack_head] = pair(e, true); |
|
517 goto next_step_1; |
|
518 } |
|
519 _graph.firstIn(e, start); |
|
520 while ( e != INVALID && ((*_flow)[e] == 0 || |
|
521 !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) ) |
|
522 _graph.nextIn(e); |
|
523 if (e != INVALID) { |
|
524 stack[++stack_head] = pair(e, false); |
|
525 goto next_step_1; |
|
526 } |
|
527 processed[start] = true; |
|
528 proc_vector[++proc_head] = start; |
|
529 continue; |
|
530 next_step_1: |
|
531 |
|
532 while (stack_head >= 0) { |
|
533 Arc se = stack[stack_head].first; |
|
534 bool sf = stack[stack_head].second; |
|
535 Node u, v; |
|
536 if (sf) { |
|
537 u = _graph.source(se); |
|
538 v = _graph.target(se); |
|
539 } else { |
|
540 u = _graph.target(se); |
|
541 v = _graph.source(se); |
|
542 } |
|
543 |
|
544 if (!reached[v]) { |
|
545 // A new node is reached |
|
546 reached[v] = true; |
|
547 pred_node[v] = u; |
|
548 pred_arc[v] = se; |
|
549 // Find the first admissible residual outgoing arc |
|
550 double p = pi[v]; |
|
551 Arc e; |
|
552 _graph.firstOut(e, v); |
|
553 while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 || |
|
554 !tol.negative(_cost[e] + p - pi[_graph.target(e)])) ) |
|
555 _graph.nextOut(e); |
|
556 if (e != INVALID) { |
|
557 stack[++stack_head] = pair(e, true); |
|
558 goto next_step_2; |
|
559 } |
|
560 _graph.firstIn(e, v); |
|
561 while ( e != INVALID && ((*_flow)[e] == 0 || |
|
562 !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) ) |
|
563 _graph.nextIn(e); |
|
564 stack[++stack_head] = pair(e, false); |
|
565 next_step_2: ; |
|
566 } else { |
|
567 if (!processed[v]) { |
|
568 // A cycle is found |
|
569 Node n, w = u; |
|
570 Capacity d, delta = sf ? _capacity[se] - (*_flow)[se] : |
|
571 (*_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]]; |
|
576 if (d <= delta) { |
|
577 delta = d; |
|
578 w = pred_node[n]; |
|
579 } |
|
580 } |
|
581 |
|
582 /* |
|
583 std::cout << "CYCLE FOUND: "; |
|
584 if (sf) |
|
585 std::cout << _cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)]; |
|
586 else |
|
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])]; |
|
591 else |
|
592 std::cout << " " << -(_cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])]); |
|
593 } |
|
594 std::cout << "\n"; |
|
595 */ |
|
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; |
|
602 else |
|
603 (*_flow)[pred_arc[n]] -= delta; |
|
604 } |
|
605 for (n = u; stack_head > 0 && n != w; n = pred_node[n]) { |
|
606 --stack_head; |
|
607 reached[n] = false; |
|
608 } |
|
609 u = w; |
|
610 } |
|
611 v = u; |
|
612 |
|
613 // Find the next admissible residual outgoing arc |
|
614 double p = pi[v]; |
|
615 Arc e = stack[stack_head].first; |
|
616 if (!stack[stack_head].second) { |
|
617 _graph.nextIn(e); |
|
618 goto in_arc_3; |
|
619 } |
|
620 _graph.nextOut(e); |
|
621 while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 || |
|
622 !tol.negative(_cost[e] + p - pi[_graph.target(e)])) ) |
|
623 _graph.nextOut(e); |
|
624 if (e != INVALID) { |
|
625 stack[stack_head] = pair(e, true); |
|
626 goto next_step_3; |
|
627 } |
|
628 _graph.firstIn(e, v); |
|
629 in_arc_3: |
|
630 while ( e != INVALID && ((*_flow)[e] == 0 || |
|
631 !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) ) |
|
632 _graph.nextIn(e); |
|
633 stack[stack_head] = pair(e, false); |
|
634 next_step_3: ; |
|
635 } |
|
636 |
|
637 while (stack_head >= 0 && stack[stack_head].first == INVALID) { |
|
638 processed[v] = true; |
|
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 |
|
645 double p = pi[v]; |
|
646 Arc e = stack[stack_head].first; |
|
647 if (!stack[stack_head].second) { |
|
648 _graph.nextIn(e); |
|
649 goto in_arc_4; |
|
650 } |
|
651 _graph.nextOut(e); |
|
652 while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 || |
|
653 !tol.negative(_cost[e] + p - pi[_graph.target(e)])) ) |
|
654 _graph.nextOut(e); |
|
655 if (e != INVALID) { |
|
656 stack[stack_head] = pair(e, true); |
|
657 goto next_step_4; |
|
658 } |
|
659 _graph.firstIn(e, v); |
|
660 in_arc_4: |
|
661 while ( e != INVALID && ((*_flow)[e] == 0 || |
|
662 !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) ) |
|
663 _graph.nextIn(e); |
|
664 stack[stack_head] = pair(e, false); |
|
665 next_step_4: ; |
|
666 } |
|
667 } |
|
668 } |
|
669 } |
|
670 t1.stop(); |
|
671 |
|
672 // Tighten potentials and epsilon |
|
673 if (--iter > 0) { |
|
674 // Compute levels |
|
675 t2.start(); |
|
676 for (int i = proc_head; i >= 0; --i) { |
|
677 Node v = proc_vector[i]; |
|
678 double p = pi[v]; |
|
679 int l = 0; |
|
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; |
|
685 } |
|
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; |
|
691 } |
|
692 level[v] = l; |
|
693 } |
|
694 |
|
695 // Modify potentials |
|
696 double p, q = -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; |
|
704 } |
|
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; |
|
709 } |
|
710 } |
|
711 for (NodeIt v(_graph); v != INVALID; ++v) { |
|
712 pi[v] -= q * level[v]; |
|
713 } |
|
714 |
|
715 // Modify epsilon |
|
716 epsilon = 0; |
|
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) |
|
721 epsilon = -curr; |
|
722 else if ((*_flow)[e] > 0 && curr > epsilon) |
|
723 epsilon = curr; |
|
724 } |
|
725 t2.stop(); |
|
726 } else { |
|
727 // Set epsilon to the minimum cycle mean |
|
728 t3.start(); |
|
729 |
|
730 /**/ |
|
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); |
|
737 |
|
738 for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e) |
|
739 static_cost[arc_ref[e]] = _res_cost[e]; |
|
740 |
|
741 Howard<StaticDigraph, typename StaticDigraph::template ArcMap<double> > |
|
742 mmc(static_graph, static_cost); |
|
743 mmc.findMinMean(); |
|
744 epsilon = -mmc.cycleMean(); |
|
745 /**/ |
|
746 |
|
747 /* |
|
748 Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost); |
|
749 mmc.findMinMean(); |
|
750 epsilon = -mmc.cycleMean(); |
|
751 */ |
|
752 |
|
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); |
|
761 bf.start(); |
|
762 for (NodeIt n(_graph); n != INVALID; ++n) |
|
763 pi[n] = static_pi[node_ref[n]]; |
|
764 |
|
765 /* |
|
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); |
|
773 bf.start(); |
|
774 */ |
|
775 |
|
776 iter = limit; |
|
777 t3.stop(); |
|
778 } |
|
779 } |
|
780 |
|
781 // std::cout << t1.realTime() << " " << t2.realTime() << " " << t3.realTime() << "\n"; |
|
782 |
|
783 // Handle non-zero lower bounds |
|
784 if (_lower) { |
|
785 for (ArcIt e(_graph); e != INVALID; ++e) |
|
786 (*_flow)[e] += (*_lower)[e]; |
|
787 } |
|
788 return true; |
|
789 } |
|
790 |
|
791 }; //class CancelAndTighten |
|
792 |
|
793 ///@} |
|
794 |
|
795 } //namespace lemon |
|
796 |
|
797 #endif //LEMON_CANCEL_AND_TIGHTEN_H |
|