|
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_CAPACITY_SCALING_H |
|
20 #define LEMON_CAPACITY_SCALING_H |
|
21 |
|
22 /// \ingroup min_cost_flow |
|
23 /// |
|
24 /// \file |
|
25 /// \brief Capacity scaling algorithm for finding a minimum cost flow. |
|
26 |
|
27 #include <vector> |
|
28 #include <lemon/bin_heap.h> |
|
29 |
|
30 namespace lemon { |
|
31 |
|
32 /// \addtogroup min_cost_flow |
|
33 /// @{ |
|
34 |
|
35 /// \brief Implementation of the capacity scaling algorithm for |
|
36 /// finding a minimum cost flow. |
|
37 /// |
|
38 /// \ref CapacityScaling implements the capacity scaling version |
|
39 /// of the successive shortest path algorithm for finding a minimum |
|
40 /// cost flow. |
|
41 /// |
|
42 /// \tparam Digraph The digraph type the algorithm runs on. |
|
43 /// \tparam LowerMap The type of the lower bound map. |
|
44 /// \tparam CapacityMap The type of the capacity (upper bound) map. |
|
45 /// \tparam CostMap The type of the cost (length) map. |
|
46 /// \tparam SupplyMap The type of the supply map. |
|
47 /// |
|
48 /// \warning |
|
49 /// - Arc capacities and costs should be \e non-negative \e integers. |
|
50 /// - Supply values should be \e signed \e integers. |
|
51 /// - The value types of the maps should be convertible to each other. |
|
52 /// - \c CostMap::Value must be signed type. |
|
53 /// |
|
54 /// \author Peter Kovacs |
|
55 template < typename Digraph, |
|
56 typename LowerMap = typename Digraph::template ArcMap<int>, |
|
57 typename CapacityMap = typename Digraph::template ArcMap<int>, |
|
58 typename CostMap = typename Digraph::template ArcMap<int>, |
|
59 typename SupplyMap = typename Digraph::template NodeMap<int> > |
|
60 class CapacityScaling |
|
61 { |
|
62 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
|
63 |
|
64 typedef typename CapacityMap::Value Capacity; |
|
65 typedef typename CostMap::Value Cost; |
|
66 typedef typename SupplyMap::Value Supply; |
|
67 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap; |
|
68 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap; |
|
69 typedef typename Digraph::template NodeMap<Arc> PredMap; |
|
70 |
|
71 public: |
|
72 |
|
73 /// The type of the flow map. |
|
74 typedef typename Digraph::template ArcMap<Capacity> FlowMap; |
|
75 /// The type of the potential map. |
|
76 typedef typename Digraph::template NodeMap<Cost> PotentialMap; |
|
77 |
|
78 private: |
|
79 |
|
80 /// \brief Special implementation of the \ref Dijkstra algorithm |
|
81 /// for finding shortest paths in the residual network. |
|
82 /// |
|
83 /// \ref ResidualDijkstra is a special implementation of the |
|
84 /// \ref Dijkstra algorithm for finding shortest paths in the |
|
85 /// residual network of the digraph with respect to the reduced arc |
|
86 /// costs and modifying the node potentials according to the |
|
87 /// distance of the nodes. |
|
88 class ResidualDijkstra |
|
89 { |
|
90 typedef typename Digraph::template NodeMap<int> HeapCrossRef; |
|
91 typedef BinHeap<Cost, HeapCrossRef> Heap; |
|
92 |
|
93 private: |
|
94 |
|
95 // The digraph the algorithm runs on |
|
96 const Digraph &_graph; |
|
97 |
|
98 // The main maps |
|
99 const FlowMap &_flow; |
|
100 const CapacityArcMap &_res_cap; |
|
101 const CostMap &_cost; |
|
102 const SupplyNodeMap &_excess; |
|
103 PotentialMap &_potential; |
|
104 |
|
105 // The distance map |
|
106 PotentialMap _dist; |
|
107 // The pred arc map |
|
108 PredMap &_pred; |
|
109 // The processed (i.e. permanently labeled) nodes |
|
110 std::vector<Node> _proc_nodes; |
|
111 |
|
112 public: |
|
113 |
|
114 /// Constructor. |
|
115 ResidualDijkstra( const Digraph &digraph, |
|
116 const FlowMap &flow, |
|
117 const CapacityArcMap &res_cap, |
|
118 const CostMap &cost, |
|
119 const SupplyMap &excess, |
|
120 PotentialMap &potential, |
|
121 PredMap &pred ) : |
|
122 _graph(digraph), _flow(flow), _res_cap(res_cap), _cost(cost), |
|
123 _excess(excess), _potential(potential), _dist(digraph), |
|
124 _pred(pred) |
|
125 {} |
|
126 |
|
127 /// Run the algorithm from the given source node. |
|
128 Node run(Node s, Capacity delta = 1) { |
|
129 HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); |
|
130 Heap heap(heap_cross_ref); |
|
131 heap.push(s, 0); |
|
132 _pred[s] = INVALID; |
|
133 _proc_nodes.clear(); |
|
134 |
|
135 // Processing nodes |
|
136 while (!heap.empty() && _excess[heap.top()] > -delta) { |
|
137 Node u = heap.top(), v; |
|
138 Cost d = heap.prio() + _potential[u], nd; |
|
139 _dist[u] = heap.prio(); |
|
140 heap.pop(); |
|
141 _proc_nodes.push_back(u); |
|
142 |
|
143 // Traversing outgoing arcs |
|
144 for (OutArcIt e(_graph, u); e != INVALID; ++e) { |
|
145 if (_res_cap[e] >= delta) { |
|
146 v = _graph.target(e); |
|
147 switch(heap.state(v)) { |
|
148 case Heap::PRE_HEAP: |
|
149 heap.push(v, d + _cost[e] - _potential[v]); |
|
150 _pred[v] = e; |
|
151 break; |
|
152 case Heap::IN_HEAP: |
|
153 nd = d + _cost[e] - _potential[v]; |
|
154 if (nd < heap[v]) { |
|
155 heap.decrease(v, nd); |
|
156 _pred[v] = e; |
|
157 } |
|
158 break; |
|
159 case Heap::POST_HEAP: |
|
160 break; |
|
161 } |
|
162 } |
|
163 } |
|
164 |
|
165 // Traversing incoming arcs |
|
166 for (InArcIt e(_graph, u); e != INVALID; ++e) { |
|
167 if (_flow[e] >= delta) { |
|
168 v = _graph.source(e); |
|
169 switch(heap.state(v)) { |
|
170 case Heap::PRE_HEAP: |
|
171 heap.push(v, d - _cost[e] - _potential[v]); |
|
172 _pred[v] = e; |
|
173 break; |
|
174 case Heap::IN_HEAP: |
|
175 nd = d - _cost[e] - _potential[v]; |
|
176 if (nd < heap[v]) { |
|
177 heap.decrease(v, nd); |
|
178 _pred[v] = e; |
|
179 } |
|
180 break; |
|
181 case Heap::POST_HEAP: |
|
182 break; |
|
183 } |
|
184 } |
|
185 } |
|
186 } |
|
187 if (heap.empty()) return INVALID; |
|
188 |
|
189 // Updating potentials of processed nodes |
|
190 Node t = heap.top(); |
|
191 Cost t_dist = heap.prio(); |
|
192 for (int i = 0; i < int(_proc_nodes.size()); ++i) |
|
193 _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; |
|
194 |
|
195 return t; |
|
196 } |
|
197 |
|
198 }; //class ResidualDijkstra |
|
199 |
|
200 private: |
|
201 |
|
202 // The digraph the algorithm runs on |
|
203 const Digraph &_graph; |
|
204 // The original lower bound map |
|
205 const LowerMap *_lower; |
|
206 // The modified capacity map |
|
207 CapacityArcMap _capacity; |
|
208 // The original cost map |
|
209 const CostMap &_cost; |
|
210 // The modified supply map |
|
211 SupplyNodeMap _supply; |
|
212 bool _valid_supply; |
|
213 |
|
214 // Arc map of the current flow |
|
215 FlowMap *_flow; |
|
216 bool _local_flow; |
|
217 // Node map of the current potentials |
|
218 PotentialMap *_potential; |
|
219 bool _local_potential; |
|
220 |
|
221 // The residual capacity map |
|
222 CapacityArcMap _res_cap; |
|
223 // The excess map |
|
224 SupplyNodeMap _excess; |
|
225 // The excess nodes (i.e. nodes with positive excess) |
|
226 std::vector<Node> _excess_nodes; |
|
227 // The deficit nodes (i.e. nodes with negative excess) |
|
228 std::vector<Node> _deficit_nodes; |
|
229 |
|
230 // The delta parameter used for capacity scaling |
|
231 Capacity _delta; |
|
232 // The maximum number of phases |
|
233 int _phase_num; |
|
234 |
|
235 // The pred arc map |
|
236 PredMap _pred; |
|
237 // Implementation of the Dijkstra algorithm for finding augmenting |
|
238 // shortest paths in the residual network |
|
239 ResidualDijkstra *_dijkstra; |
|
240 |
|
241 public: |
|
242 |
|
243 /// \brief General constructor (with lower bounds). |
|
244 /// |
|
245 /// General constructor (with lower bounds). |
|
246 /// |
|
247 /// \param digraph The digraph the algorithm runs on. |
|
248 /// \param lower The lower bounds of the arcs. |
|
249 /// \param capacity The capacities (upper bounds) of the arcs. |
|
250 /// \param cost The cost (length) values of the arcs. |
|
251 /// \param supply The supply values of the nodes (signed). |
|
252 CapacityScaling( const Digraph &digraph, |
|
253 const LowerMap &lower, |
|
254 const CapacityMap &capacity, |
|
255 const CostMap &cost, |
|
256 const SupplyMap &supply ) : |
|
257 _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost), |
|
258 _supply(digraph), _flow(NULL), _local_flow(false), |
|
259 _potential(NULL), _local_potential(false), |
|
260 _res_cap(digraph), _excess(digraph), _pred(digraph), _dijkstra(NULL) |
|
261 { |
|
262 Supply sum = 0; |
|
263 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
264 _supply[n] = supply[n]; |
|
265 _excess[n] = supply[n]; |
|
266 sum += supply[n]; |
|
267 } |
|
268 _valid_supply = sum == 0; |
|
269 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
270 _capacity[a] = capacity[a]; |
|
271 _res_cap[a] = capacity[a]; |
|
272 } |
|
273 |
|
274 // Remove non-zero lower bounds |
|
275 typename LowerMap::Value lcap; |
|
276 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
277 if ((lcap = lower[e]) != 0) { |
|
278 _capacity[e] -= lcap; |
|
279 _res_cap[e] -= lcap; |
|
280 _supply[_graph.source(e)] -= lcap; |
|
281 _supply[_graph.target(e)] += lcap; |
|
282 _excess[_graph.source(e)] -= lcap; |
|
283 _excess[_graph.target(e)] += lcap; |
|
284 } |
|
285 } |
|
286 } |
|
287 /* |
|
288 /// \brief General constructor (without lower bounds). |
|
289 /// |
|
290 /// General constructor (without lower bounds). |
|
291 /// |
|
292 /// \param digraph The digraph the algorithm runs on. |
|
293 /// \param capacity The capacities (upper bounds) of the arcs. |
|
294 /// \param cost The cost (length) values of the arcs. |
|
295 /// \param supply The supply values of the nodes (signed). |
|
296 CapacityScaling( const Digraph &digraph, |
|
297 const CapacityMap &capacity, |
|
298 const CostMap &cost, |
|
299 const SupplyMap &supply ) : |
|
300 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), |
|
301 _supply(supply), _flow(NULL), _local_flow(false), |
|
302 _potential(NULL), _local_potential(false), |
|
303 _res_cap(capacity), _excess(supply), _pred(digraph), _dijkstra(NULL) |
|
304 { |
|
305 // Check the sum of supply values |
|
306 Supply sum = 0; |
|
307 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; |
|
308 _valid_supply = sum == 0; |
|
309 } |
|
310 |
|
311 /// \brief Simple constructor (with lower bounds). |
|
312 /// |
|
313 /// Simple constructor (with lower bounds). |
|
314 /// |
|
315 /// \param digraph The digraph the algorithm runs on. |
|
316 /// \param lower The lower bounds of the arcs. |
|
317 /// \param capacity The capacities (upper bounds) of the arcs. |
|
318 /// \param cost The cost (length) values of the arcs. |
|
319 /// \param s The source node. |
|
320 /// \param t The target node. |
|
321 /// \param flow_value The required amount of flow from node \c s |
|
322 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
|
323 CapacityScaling( const Digraph &digraph, |
|
324 const LowerMap &lower, |
|
325 const CapacityMap &capacity, |
|
326 const CostMap &cost, |
|
327 Node s, Node t, |
|
328 Supply flow_value ) : |
|
329 _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost), |
|
330 _supply(digraph, 0), _flow(NULL), _local_flow(false), |
|
331 _potential(NULL), _local_potential(false), |
|
332 _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL) |
|
333 { |
|
334 // Remove non-zero lower bounds |
|
335 _supply[s] = _excess[s] = flow_value; |
|
336 _supply[t] = _excess[t] = -flow_value; |
|
337 typename LowerMap::Value lcap; |
|
338 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
339 if ((lcap = lower[e]) != 0) { |
|
340 _capacity[e] -= lcap; |
|
341 _res_cap[e] -= lcap; |
|
342 _supply[_graph.source(e)] -= lcap; |
|
343 _supply[_graph.target(e)] += lcap; |
|
344 _excess[_graph.source(e)] -= lcap; |
|
345 _excess[_graph.target(e)] += lcap; |
|
346 } |
|
347 } |
|
348 _valid_supply = true; |
|
349 } |
|
350 |
|
351 /// \brief Simple constructor (without lower bounds). |
|
352 /// |
|
353 /// Simple constructor (without lower bounds). |
|
354 /// |
|
355 /// \param digraph The digraph the algorithm runs on. |
|
356 /// \param capacity The capacities (upper bounds) of the arcs. |
|
357 /// \param cost The cost (length) values of the arcs. |
|
358 /// \param s The source node. |
|
359 /// \param t The target node. |
|
360 /// \param flow_value The required amount of flow from node \c s |
|
361 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
|
362 CapacityScaling( const Digraph &digraph, |
|
363 const CapacityMap &capacity, |
|
364 const CostMap &cost, |
|
365 Node s, Node t, |
|
366 Supply flow_value ) : |
|
367 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), |
|
368 _supply(digraph, 0), _flow(NULL), _local_flow(false), |
|
369 _potential(NULL), _local_potential(false), |
|
370 _res_cap(capacity), _excess(digraph, 0), _pred(digraph), _dijkstra(NULL) |
|
371 { |
|
372 _supply[s] = _excess[s] = flow_value; |
|
373 _supply[t] = _excess[t] = -flow_value; |
|
374 _valid_supply = true; |
|
375 } |
|
376 */ |
|
377 /// Destructor. |
|
378 ~CapacityScaling() { |
|
379 if (_local_flow) delete _flow; |
|
380 if (_local_potential) delete _potential; |
|
381 delete _dijkstra; |
|
382 } |
|
383 |
|
384 /// \brief Set the flow map. |
|
385 /// |
|
386 /// Set the flow map. |
|
387 /// |
|
388 /// \return \c (*this) |
|
389 CapacityScaling& flowMap(FlowMap &map) { |
|
390 if (_local_flow) { |
|
391 delete _flow; |
|
392 _local_flow = false; |
|
393 } |
|
394 _flow = ↦ |
|
395 return *this; |
|
396 } |
|
397 |
|
398 /// \brief Set the potential map. |
|
399 /// |
|
400 /// Set the potential map. |
|
401 /// |
|
402 /// \return \c (*this) |
|
403 CapacityScaling& potentialMap(PotentialMap &map) { |
|
404 if (_local_potential) { |
|
405 delete _potential; |
|
406 _local_potential = false; |
|
407 } |
|
408 _potential = ↦ |
|
409 return *this; |
|
410 } |
|
411 |
|
412 /// \name Execution control |
|
413 |
|
414 /// @{ |
|
415 |
|
416 /// \brief Run the algorithm. |
|
417 /// |
|
418 /// This function runs the algorithm. |
|
419 /// |
|
420 /// \param scaling Enable or disable capacity scaling. |
|
421 /// If the maximum arc capacity and/or the amount of total supply |
|
422 /// is rather small, the algorithm could be slightly faster without |
|
423 /// scaling. |
|
424 /// |
|
425 /// \return \c true if a feasible flow can be found. |
|
426 bool run(bool scaling = true) { |
|
427 return init(scaling) && start(); |
|
428 } |
|
429 |
|
430 /// @} |
|
431 |
|
432 /// \name Query Functions |
|
433 /// The results of the algorithm can be obtained using these |
|
434 /// functions.\n |
|
435 /// \ref lemon::CapacityScaling::run() "run()" must be called before |
|
436 /// using them. |
|
437 |
|
438 /// @{ |
|
439 |
|
440 /// \brief Return a const reference to the arc map storing the |
|
441 /// found flow. |
|
442 /// |
|
443 /// Return a const reference to the arc map storing the found flow. |
|
444 /// |
|
445 /// \pre \ref run() must be called before using this function. |
|
446 const FlowMap& flowMap() const { |
|
447 return *_flow; |
|
448 } |
|
449 |
|
450 /// \brief Return a const reference to the node map storing the |
|
451 /// found potentials (the dual solution). |
|
452 /// |
|
453 /// Return a const reference to the node map storing the found |
|
454 /// potentials (the dual solution). |
|
455 /// |
|
456 /// \pre \ref run() must be called before using this function. |
|
457 const PotentialMap& potentialMap() const { |
|
458 return *_potential; |
|
459 } |
|
460 |
|
461 /// \brief Return the flow on the given arc. |
|
462 /// |
|
463 /// Return the flow on the given arc. |
|
464 /// |
|
465 /// \pre \ref run() must be called before using this function. |
|
466 Capacity flow(const Arc& arc) const { |
|
467 return (*_flow)[arc]; |
|
468 } |
|
469 |
|
470 /// \brief Return the potential of the given node. |
|
471 /// |
|
472 /// Return the potential of the given node. |
|
473 /// |
|
474 /// \pre \ref run() must be called before using this function. |
|
475 Cost potential(const Node& node) const { |
|
476 return (*_potential)[node]; |
|
477 } |
|
478 |
|
479 /// \brief Return the total cost of the found flow. |
|
480 /// |
|
481 /// Return the total cost of the found flow. The complexity of the |
|
482 /// function is \f$ O(e) \f$. |
|
483 /// |
|
484 /// \pre \ref run() must be called before using this function. |
|
485 Cost totalCost() const { |
|
486 Cost c = 0; |
|
487 for (ArcIt e(_graph); e != INVALID; ++e) |
|
488 c += (*_flow)[e] * _cost[e]; |
|
489 return c; |
|
490 } |
|
491 |
|
492 /// @} |
|
493 |
|
494 private: |
|
495 |
|
496 /// Initialize the algorithm. |
|
497 bool init(bool scaling) { |
|
498 if (!_valid_supply) return false; |
|
499 |
|
500 // Initializing maps |
|
501 if (!_flow) { |
|
502 _flow = new FlowMap(_graph); |
|
503 _local_flow = true; |
|
504 } |
|
505 if (!_potential) { |
|
506 _potential = new PotentialMap(_graph); |
|
507 _local_potential = true; |
|
508 } |
|
509 for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0; |
|
510 for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0; |
|
511 |
|
512 _dijkstra = new ResidualDijkstra( _graph, *_flow, _res_cap, _cost, |
|
513 _excess, *_potential, _pred ); |
|
514 |
|
515 // Initializing delta value |
|
516 if (scaling) { |
|
517 // With scaling |
|
518 Supply max_sup = 0, max_dem = 0; |
|
519 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
520 if ( _supply[n] > max_sup) max_sup = _supply[n]; |
|
521 if (-_supply[n] > max_dem) max_dem = -_supply[n]; |
|
522 } |
|
523 Capacity max_cap = 0; |
|
524 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
525 if (_capacity[e] > max_cap) max_cap = _capacity[e]; |
|
526 } |
|
527 max_sup = std::min(std::min(max_sup, max_dem), max_cap); |
|
528 _phase_num = 0; |
|
529 for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) |
|
530 ++_phase_num; |
|
531 } else { |
|
532 // Without scaling |
|
533 _delta = 1; |
|
534 } |
|
535 |
|
536 return true; |
|
537 } |
|
538 |
|
539 bool start() { |
|
540 if (_delta > 1) |
|
541 return startWithScaling(); |
|
542 else |
|
543 return startWithoutScaling(); |
|
544 } |
|
545 |
|
546 /// Execute the capacity scaling algorithm. |
|
547 bool startWithScaling() { |
|
548 // Processing capacity scaling phases |
|
549 Node s, t; |
|
550 int phase_cnt = 0; |
|
551 int factor = 4; |
|
552 while (true) { |
|
553 // Saturating all arcs not satisfying the optimality condition |
|
554 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
555 Node u = _graph.source(e), v = _graph.target(e); |
|
556 Cost c = _cost[e] + (*_potential)[u] - (*_potential)[v]; |
|
557 if (c < 0 && _res_cap[e] >= _delta) { |
|
558 _excess[u] -= _res_cap[e]; |
|
559 _excess[v] += _res_cap[e]; |
|
560 (*_flow)[e] = _capacity[e]; |
|
561 _res_cap[e] = 0; |
|
562 } |
|
563 else if (c > 0 && (*_flow)[e] >= _delta) { |
|
564 _excess[u] += (*_flow)[e]; |
|
565 _excess[v] -= (*_flow)[e]; |
|
566 (*_flow)[e] = 0; |
|
567 _res_cap[e] = _capacity[e]; |
|
568 } |
|
569 } |
|
570 |
|
571 // Finding excess nodes and deficit nodes |
|
572 _excess_nodes.clear(); |
|
573 _deficit_nodes.clear(); |
|
574 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
575 if (_excess[n] >= _delta) _excess_nodes.push_back(n); |
|
576 if (_excess[n] <= -_delta) _deficit_nodes.push_back(n); |
|
577 } |
|
578 int next_node = 0, next_def_node = 0; |
|
579 |
|
580 // Finding augmenting shortest paths |
|
581 while (next_node < int(_excess_nodes.size())) { |
|
582 // Checking deficit nodes |
|
583 if (_delta > 1) { |
|
584 bool delta_deficit = false; |
|
585 for ( ; next_def_node < int(_deficit_nodes.size()); |
|
586 ++next_def_node ) { |
|
587 if (_excess[_deficit_nodes[next_def_node]] <= -_delta) { |
|
588 delta_deficit = true; |
|
589 break; |
|
590 } |
|
591 } |
|
592 if (!delta_deficit) break; |
|
593 } |
|
594 |
|
595 // Running Dijkstra |
|
596 s = _excess_nodes[next_node]; |
|
597 if ((t = _dijkstra->run(s, _delta)) == INVALID) { |
|
598 if (_delta > 1) { |
|
599 ++next_node; |
|
600 continue; |
|
601 } |
|
602 return false; |
|
603 } |
|
604 |
|
605 // Augmenting along a shortest path from s to t. |
|
606 Capacity d = std::min(_excess[s], -_excess[t]); |
|
607 Node u = t; |
|
608 Arc e; |
|
609 if (d > _delta) { |
|
610 while ((e = _pred[u]) != INVALID) { |
|
611 Capacity rc; |
|
612 if (u == _graph.target(e)) { |
|
613 rc = _res_cap[e]; |
|
614 u = _graph.source(e); |
|
615 } else { |
|
616 rc = (*_flow)[e]; |
|
617 u = _graph.target(e); |
|
618 } |
|
619 if (rc < d) d = rc; |
|
620 } |
|
621 } |
|
622 u = t; |
|
623 while ((e = _pred[u]) != INVALID) { |
|
624 if (u == _graph.target(e)) { |
|
625 (*_flow)[e] += d; |
|
626 _res_cap[e] -= d; |
|
627 u = _graph.source(e); |
|
628 } else { |
|
629 (*_flow)[e] -= d; |
|
630 _res_cap[e] += d; |
|
631 u = _graph.target(e); |
|
632 } |
|
633 } |
|
634 _excess[s] -= d; |
|
635 _excess[t] += d; |
|
636 |
|
637 if (_excess[s] < _delta) ++next_node; |
|
638 } |
|
639 |
|
640 if (_delta == 1) break; |
|
641 if (++phase_cnt > _phase_num / 4) factor = 2; |
|
642 _delta = _delta <= factor ? 1 : _delta / factor; |
|
643 } |
|
644 |
|
645 // Handling non-zero lower bounds |
|
646 if (_lower) { |
|
647 for (ArcIt e(_graph); e != INVALID; ++e) |
|
648 (*_flow)[e] += (*_lower)[e]; |
|
649 } |
|
650 return true; |
|
651 } |
|
652 |
|
653 /// Execute the successive shortest path algorithm. |
|
654 bool startWithoutScaling() { |
|
655 // Finding excess nodes |
|
656 for (NodeIt n(_graph); n != INVALID; ++n) |
|
657 if (_excess[n] > 0) _excess_nodes.push_back(n); |
|
658 if (_excess_nodes.size() == 0) return true; |
|
659 int next_node = 0; |
|
660 |
|
661 // Finding shortest paths |
|
662 Node s, t; |
|
663 while ( _excess[_excess_nodes[next_node]] > 0 || |
|
664 ++next_node < int(_excess_nodes.size()) ) |
|
665 { |
|
666 // Running Dijkstra |
|
667 s = _excess_nodes[next_node]; |
|
668 if ((t = _dijkstra->run(s)) == INVALID) return false; |
|
669 |
|
670 // Augmenting along a shortest path from s to t |
|
671 Capacity d = std::min(_excess[s], -_excess[t]); |
|
672 Node u = t; |
|
673 Arc e; |
|
674 if (d > 1) { |
|
675 while ((e = _pred[u]) != INVALID) { |
|
676 Capacity rc; |
|
677 if (u == _graph.target(e)) { |
|
678 rc = _res_cap[e]; |
|
679 u = _graph.source(e); |
|
680 } else { |
|
681 rc = (*_flow)[e]; |
|
682 u = _graph.target(e); |
|
683 } |
|
684 if (rc < d) d = rc; |
|
685 } |
|
686 } |
|
687 u = t; |
|
688 while ((e = _pred[u]) != INVALID) { |
|
689 if (u == _graph.target(e)) { |
|
690 (*_flow)[e] += d; |
|
691 _res_cap[e] -= d; |
|
692 u = _graph.source(e); |
|
693 } else { |
|
694 (*_flow)[e] -= d; |
|
695 _res_cap[e] += d; |
|
696 u = _graph.target(e); |
|
697 } |
|
698 } |
|
699 _excess[s] -= d; |
|
700 _excess[t] += d; |
|
701 } |
|
702 |
|
703 // Handling non-zero lower bounds |
|
704 if (_lower) { |
|
705 for (ArcIt e(_graph); e != INVALID; ++e) |
|
706 (*_flow)[e] += (*_lower)[e]; |
|
707 } |
|
708 return true; |
|
709 } |
|
710 |
|
711 }; //class CapacityScaling |
|
712 |
|
713 ///@} |
|
714 |
|
715 } //namespace lemon |
|
716 |
|
717 #endif //LEMON_CAPACITY_SCALING_H |