|
1 /* -*- mode: C++; indent-tabs-mode: nil; -*- |
|
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_PREFLOW_H |
|
20 #define LEMON_PREFLOW_H |
|
21 |
|
22 #include <lemon/error.h> |
|
23 #include <lemon/tolerance.h> |
|
24 #include <lemon/elevator.h> |
|
25 |
|
26 /// \file |
|
27 /// \ingroup max_flow |
|
28 /// \brief Implementation of the preflow algorithm. |
|
29 |
|
30 namespace lemon { |
|
31 |
|
32 /// \brief Default traits class of Preflow class. |
|
33 /// |
|
34 /// Default traits class of Preflow class. |
|
35 /// \param _Graph Digraph type. |
|
36 /// \param _CapacityMap Type of capacity map. |
|
37 template <typename _Graph, typename _CapacityMap> |
|
38 struct PreflowDefaultTraits { |
|
39 |
|
40 /// \brief The digraph type the algorithm runs on. |
|
41 typedef _Graph Digraph; |
|
42 |
|
43 /// \brief The type of the map that stores the arc capacities. |
|
44 /// |
|
45 /// The type of the map that stores the arc capacities. |
|
46 /// It must meet the \ref concepts::ReadMap "ReadMap" concept. |
|
47 typedef _CapacityMap CapacityMap; |
|
48 |
|
49 /// \brief The type of the length of the arcs. |
|
50 typedef typename CapacityMap::Value Value; |
|
51 |
|
52 /// \brief The map type that stores the flow values. |
|
53 /// |
|
54 /// The map type that stores the flow values. |
|
55 /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. |
|
56 typedef typename Digraph::template ArcMap<Value> FlowMap; |
|
57 |
|
58 /// \brief Instantiates a FlowMap. |
|
59 /// |
|
60 /// This function instantiates a \ref FlowMap. |
|
61 /// \param digraph The digraph, to which we would like to define |
|
62 /// the flow map. |
|
63 static FlowMap* createFlowMap(const Digraph& digraph) { |
|
64 return new FlowMap(digraph); |
|
65 } |
|
66 |
|
67 /// \brief The eleavator type used by Preflow algorithm. |
|
68 /// |
|
69 /// The elevator type used by Preflow algorithm. |
|
70 /// |
|
71 /// \sa Elevator |
|
72 /// \sa LinkedElevator |
|
73 typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator; |
|
74 |
|
75 /// \brief Instantiates an Elevator. |
|
76 /// |
|
77 /// This function instantiates a \ref Elevator. |
|
78 /// \param digraph The digraph, to which we would like to define |
|
79 /// the elevator. |
|
80 /// \param max_level The maximum level of the elevator. |
|
81 static Elevator* createElevator(const Digraph& digraph, int max_level) { |
|
82 return new Elevator(digraph, max_level); |
|
83 } |
|
84 |
|
85 /// \brief The tolerance used by the algorithm |
|
86 /// |
|
87 /// The tolerance used by the algorithm to handle inexact computation. |
|
88 typedef lemon::Tolerance<Value> Tolerance; |
|
89 |
|
90 }; |
|
91 |
|
92 |
|
93 /// \ingroup max_flow |
|
94 /// |
|
95 /// \brief %Preflow algorithms class. |
|
96 /// |
|
97 /// This class provides an implementation of the Goldberg's \e |
|
98 /// preflow \e algorithm producing a flow of maximum value in a |
|
99 /// digraph. The preflow algorithms are the fastest known max |
|
100 /// flow algorithms. The current implementation use a mixture of the |
|
101 /// \e "highest label" and the \e "bound decrease" heuristics. |
|
102 /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$. |
|
103 /// |
|
104 /// The algorithm consists from two phases. After the first phase |
|
105 /// the maximal flow value and the minimum cut can be obtained. The |
|
106 /// second phase constructs the feasible maximum flow on each arc. |
|
107 /// |
|
108 /// \param _Graph The digraph type the algorithm runs on. |
|
109 /// \param _CapacityMap The flow map type. |
|
110 /// \param _Traits Traits class to set various data types used by |
|
111 /// the algorithm. The default traits class is \ref |
|
112 /// PreflowDefaultTraits. See \ref PreflowDefaultTraits for the |
|
113 /// documentation of a %Preflow traits class. |
|
114 /// |
|
115 ///\author Jacint Szabo and Balazs Dezso |
|
116 #ifdef DOXYGEN |
|
117 template <typename _Graph, typename _CapacityMap, typename _Traits> |
|
118 #else |
|
119 template <typename _Graph, |
|
120 typename _CapacityMap = typename _Graph::template ArcMap<int>, |
|
121 typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> > |
|
122 #endif |
|
123 class Preflow { |
|
124 public: |
|
125 |
|
126 typedef _Traits Traits; |
|
127 typedef typename Traits::Digraph Digraph; |
|
128 typedef typename Traits::CapacityMap CapacityMap; |
|
129 typedef typename Traits::Value Value; |
|
130 |
|
131 typedef typename Traits::FlowMap FlowMap; |
|
132 typedef typename Traits::Elevator Elevator; |
|
133 typedef typename Traits::Tolerance Tolerance; |
|
134 |
|
135 /// \brief \ref Exception for uninitialized parameters. |
|
136 /// |
|
137 /// This error represents problems in the initialization |
|
138 /// of the parameters of the algorithms. |
|
139 class UninitializedParameter : public lemon::Exception { |
|
140 public: |
|
141 virtual const char* what() const throw() { |
|
142 return "lemon::Preflow::UninitializedParameter"; |
|
143 } |
|
144 }; |
|
145 |
|
146 private: |
|
147 |
|
148 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
|
149 |
|
150 const Digraph& _graph; |
|
151 const CapacityMap* _capacity; |
|
152 |
|
153 int _node_num; |
|
154 |
|
155 Node _source, _target; |
|
156 |
|
157 FlowMap* _flow; |
|
158 bool _local_flow; |
|
159 |
|
160 Elevator* _level; |
|
161 bool _local_level; |
|
162 |
|
163 typedef typename Digraph::template NodeMap<Value> ExcessMap; |
|
164 ExcessMap* _excess; |
|
165 |
|
166 Tolerance _tolerance; |
|
167 |
|
168 bool _phase; |
|
169 |
|
170 |
|
171 void createStructures() { |
|
172 _node_num = countNodes(_graph); |
|
173 |
|
174 if (!_flow) { |
|
175 _flow = Traits::createFlowMap(_graph); |
|
176 _local_flow = true; |
|
177 } |
|
178 if (!_level) { |
|
179 _level = Traits::createElevator(_graph, _node_num); |
|
180 _local_level = true; |
|
181 } |
|
182 if (!_excess) { |
|
183 _excess = new ExcessMap(_graph); |
|
184 } |
|
185 } |
|
186 |
|
187 void destroyStructures() { |
|
188 if (_local_flow) { |
|
189 delete _flow; |
|
190 } |
|
191 if (_local_level) { |
|
192 delete _level; |
|
193 } |
|
194 if (_excess) { |
|
195 delete _excess; |
|
196 } |
|
197 } |
|
198 |
|
199 public: |
|
200 |
|
201 typedef Preflow Create; |
|
202 |
|
203 ///\name Named template parameters |
|
204 |
|
205 ///@{ |
|
206 |
|
207 template <typename _FlowMap> |
|
208 struct DefFlowMapTraits : public Traits { |
|
209 typedef _FlowMap FlowMap; |
|
210 static FlowMap *createFlowMap(const Digraph&) { |
|
211 throw UninitializedParameter(); |
|
212 } |
|
213 }; |
|
214 |
|
215 /// \brief \ref named-templ-param "Named parameter" for setting |
|
216 /// FlowMap type |
|
217 /// |
|
218 /// \ref named-templ-param "Named parameter" for setting FlowMap |
|
219 /// type |
|
220 template <typename _FlowMap> |
|
221 struct DefFlowMap |
|
222 : public Preflow<Digraph, CapacityMap, DefFlowMapTraits<_FlowMap> > { |
|
223 typedef Preflow<Digraph, CapacityMap, |
|
224 DefFlowMapTraits<_FlowMap> > Create; |
|
225 }; |
|
226 |
|
227 template <typename _Elevator> |
|
228 struct DefElevatorTraits : public Traits { |
|
229 typedef _Elevator Elevator; |
|
230 static Elevator *createElevator(const Digraph&, int) { |
|
231 throw UninitializedParameter(); |
|
232 } |
|
233 }; |
|
234 |
|
235 /// \brief \ref named-templ-param "Named parameter" for setting |
|
236 /// Elevator type |
|
237 /// |
|
238 /// \ref named-templ-param "Named parameter" for setting Elevator |
|
239 /// type |
|
240 template <typename _Elevator> |
|
241 struct DefElevator |
|
242 : public Preflow<Digraph, CapacityMap, DefElevatorTraits<_Elevator> > { |
|
243 typedef Preflow<Digraph, CapacityMap, |
|
244 DefElevatorTraits<_Elevator> > Create; |
|
245 }; |
|
246 |
|
247 template <typename _Elevator> |
|
248 struct DefStandardElevatorTraits : public Traits { |
|
249 typedef _Elevator Elevator; |
|
250 static Elevator *createElevator(const Digraph& digraph, int max_level) { |
|
251 return new Elevator(digraph, max_level); |
|
252 } |
|
253 }; |
|
254 |
|
255 /// \brief \ref named-templ-param "Named parameter" for setting |
|
256 /// Elevator type |
|
257 /// |
|
258 /// \ref named-templ-param "Named parameter" for setting Elevator |
|
259 /// type. The Elevator should be standard constructor interface, ie. |
|
260 /// the digraph and the maximum level should be passed to it. |
|
261 template <typename _Elevator> |
|
262 struct DefStandardElevator |
|
263 : public Preflow<Digraph, CapacityMap, |
|
264 DefStandardElevatorTraits<_Elevator> > { |
|
265 typedef Preflow<Digraph, CapacityMap, |
|
266 DefStandardElevatorTraits<_Elevator> > Create; |
|
267 }; |
|
268 |
|
269 /// @} |
|
270 |
|
271 protected: |
|
272 |
|
273 Preflow() {} |
|
274 |
|
275 public: |
|
276 |
|
277 |
|
278 /// \brief The constructor of the class. |
|
279 /// |
|
280 /// The constructor of the class. |
|
281 /// \param digraph The digraph the algorithm runs on. |
|
282 /// \param capacity The capacity of the arcs. |
|
283 /// \param source The source node. |
|
284 /// \param target The target node. |
|
285 Preflow(const Digraph& digraph, const CapacityMap& capacity, |
|
286 Node source, Node target) |
|
287 : _graph(digraph), _capacity(&capacity), |
|
288 _node_num(0), _source(source), _target(target), |
|
289 _flow(0), _local_flow(false), |
|
290 _level(0), _local_level(false), |
|
291 _excess(0), _tolerance(), _phase() {} |
|
292 |
|
293 /// \brief Destrcutor. |
|
294 /// |
|
295 /// Destructor. |
|
296 ~Preflow() { |
|
297 destroyStructures(); |
|
298 } |
|
299 |
|
300 /// \brief Sets the capacity map. |
|
301 /// |
|
302 /// Sets the capacity map. |
|
303 /// \return \c (*this) |
|
304 Preflow& capacityMap(const CapacityMap& map) { |
|
305 _capacity = ↦ |
|
306 return *this; |
|
307 } |
|
308 |
|
309 /// \brief Sets the flow map. |
|
310 /// |
|
311 /// Sets the flow map. |
|
312 /// \return \c (*this) |
|
313 Preflow& flowMap(FlowMap& map) { |
|
314 if (_local_flow) { |
|
315 delete _flow; |
|
316 _local_flow = false; |
|
317 } |
|
318 _flow = ↦ |
|
319 return *this; |
|
320 } |
|
321 |
|
322 /// \brief Returns the flow map. |
|
323 /// |
|
324 /// \return The flow map. |
|
325 const FlowMap& flowMap() { |
|
326 return *_flow; |
|
327 } |
|
328 |
|
329 /// \brief Sets the elevator. |
|
330 /// |
|
331 /// Sets the elevator. |
|
332 /// \return \c (*this) |
|
333 Preflow& elevator(Elevator& elevator) { |
|
334 if (_local_level) { |
|
335 delete _level; |
|
336 _local_level = false; |
|
337 } |
|
338 _level = &elevator; |
|
339 return *this; |
|
340 } |
|
341 |
|
342 /// \brief Returns the elevator. |
|
343 /// |
|
344 /// \return The elevator. |
|
345 const Elevator& elevator() { |
|
346 return *_level; |
|
347 } |
|
348 |
|
349 /// \brief Sets the source node. |
|
350 /// |
|
351 /// Sets the source node. |
|
352 /// \return \c (*this) |
|
353 Preflow& source(const Node& node) { |
|
354 _source = node; |
|
355 return *this; |
|
356 } |
|
357 |
|
358 /// \brief Sets the target node. |
|
359 /// |
|
360 /// Sets the target node. |
|
361 /// \return \c (*this) |
|
362 Preflow& target(const Node& node) { |
|
363 _target = node; |
|
364 return *this; |
|
365 } |
|
366 |
|
367 /// \brief Sets the tolerance used by algorithm. |
|
368 /// |
|
369 /// Sets the tolerance used by algorithm. |
|
370 Preflow& tolerance(const Tolerance& tolerance) const { |
|
371 _tolerance = tolerance; |
|
372 return *this; |
|
373 } |
|
374 |
|
375 /// \brief Returns the tolerance used by algorithm. |
|
376 /// |
|
377 /// Returns the tolerance used by algorithm. |
|
378 const Tolerance& tolerance() const { |
|
379 return tolerance; |
|
380 } |
|
381 |
|
382 /// \name Execution control The simplest way to execute the |
|
383 /// algorithm is to use one of the member functions called \c |
|
384 /// run(). |
|
385 /// \n |
|
386 /// If you need more control on initial solution or |
|
387 /// execution then you have to call one \ref init() function and then |
|
388 /// the startFirstPhase() and if you need the startSecondPhase(). |
|
389 |
|
390 ///@{ |
|
391 |
|
392 /// \brief Initializes the internal data structures. |
|
393 /// |
|
394 /// Initializes the internal data structures. |
|
395 /// |
|
396 void init() { |
|
397 createStructures(); |
|
398 |
|
399 _phase = true; |
|
400 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
401 _excess->set(n, 0); |
|
402 } |
|
403 |
|
404 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
405 _flow->set(e, 0); |
|
406 } |
|
407 |
|
408 typename Digraph::template NodeMap<bool> reached(_graph, false); |
|
409 |
|
410 _level->initStart(); |
|
411 _level->initAddItem(_target); |
|
412 |
|
413 std::vector<Node> queue; |
|
414 reached.set(_source, true); |
|
415 |
|
416 queue.push_back(_target); |
|
417 reached.set(_target, true); |
|
418 while (!queue.empty()) { |
|
419 _level->initNewLevel(); |
|
420 std::vector<Node> nqueue; |
|
421 for (int i = 0; i < int(queue.size()); ++i) { |
|
422 Node n = queue[i]; |
|
423 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
424 Node u = _graph.source(e); |
|
425 if (!reached[u] && _tolerance.positive((*_capacity)[e])) { |
|
426 reached.set(u, true); |
|
427 _level->initAddItem(u); |
|
428 nqueue.push_back(u); |
|
429 } |
|
430 } |
|
431 } |
|
432 queue.swap(nqueue); |
|
433 } |
|
434 _level->initFinish(); |
|
435 |
|
436 for (OutArcIt e(_graph, _source); e != INVALID; ++e) { |
|
437 if (_tolerance.positive((*_capacity)[e])) { |
|
438 Node u = _graph.target(e); |
|
439 if ((*_level)[u] == _level->maxLevel()) continue; |
|
440 _flow->set(e, (*_capacity)[e]); |
|
441 _excess->set(u, (*_excess)[u] + (*_capacity)[e]); |
|
442 if (u != _target && !_level->active(u)) { |
|
443 _level->activate(u); |
|
444 } |
|
445 } |
|
446 } |
|
447 } |
|
448 |
|
449 /// \brief Initializes the internal data structures. |
|
450 /// |
|
451 /// Initializes the internal data structures and sets the initial |
|
452 /// flow to the given \c flowMap. The \c flowMap should contain a |
|
453 /// flow or at least a preflow, ie. in each node excluding the |
|
454 /// target the incoming flow should greater or equal to the |
|
455 /// outgoing flow. |
|
456 /// \return %False when the given \c flowMap is not a preflow. |
|
457 template <typename FlowMap> |
|
458 bool flowInit(const FlowMap& flowMap) { |
|
459 createStructures(); |
|
460 |
|
461 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
462 _flow->set(e, flowMap[e]); |
|
463 } |
|
464 |
|
465 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
466 Value excess = 0; |
|
467 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
468 excess += (*_flow)[e]; |
|
469 } |
|
470 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
471 excess -= (*_flow)[e]; |
|
472 } |
|
473 if (excess < 0 && n != _source) return false; |
|
474 _excess->set(n, excess); |
|
475 } |
|
476 |
|
477 typename Digraph::template NodeMap<bool> reached(_graph, false); |
|
478 |
|
479 _level->initStart(); |
|
480 _level->initAddItem(_target); |
|
481 |
|
482 std::vector<Node> queue; |
|
483 reached.set(_source, true); |
|
484 |
|
485 queue.push_back(_target); |
|
486 reached.set(_target, true); |
|
487 while (!queue.empty()) { |
|
488 _level->initNewLevel(); |
|
489 std::vector<Node> nqueue; |
|
490 for (int i = 0; i < int(queue.size()); ++i) { |
|
491 Node n = queue[i]; |
|
492 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
493 Node u = _graph.source(e); |
|
494 if (!reached[u] && |
|
495 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) { |
|
496 reached.set(u, true); |
|
497 _level->initAddItem(u); |
|
498 nqueue.push_back(u); |
|
499 } |
|
500 } |
|
501 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
502 Node v = _graph.target(e); |
|
503 if (!reached[v] && _tolerance.positive((*_flow)[e])) { |
|
504 reached.set(v, true); |
|
505 _level->initAddItem(v); |
|
506 nqueue.push_back(v); |
|
507 } |
|
508 } |
|
509 } |
|
510 queue.swap(nqueue); |
|
511 } |
|
512 _level->initFinish(); |
|
513 |
|
514 for (OutArcIt e(_graph, _source); e != INVALID; ++e) { |
|
515 Value rem = (*_capacity)[e] - (*_flow)[e]; |
|
516 if (_tolerance.positive(rem)) { |
|
517 Node u = _graph.target(e); |
|
518 if ((*_level)[u] == _level->maxLevel()) continue; |
|
519 _flow->set(e, (*_capacity)[e]); |
|
520 _excess->set(u, (*_excess)[u] + rem); |
|
521 if (u != _target && !_level->active(u)) { |
|
522 _level->activate(u); |
|
523 } |
|
524 } |
|
525 } |
|
526 for (InArcIt e(_graph, _source); e != INVALID; ++e) { |
|
527 Value rem = (*_flow)[e]; |
|
528 if (_tolerance.positive(rem)) { |
|
529 Node v = _graph.source(e); |
|
530 if ((*_level)[v] == _level->maxLevel()) continue; |
|
531 _flow->set(e, 0); |
|
532 _excess->set(v, (*_excess)[v] + rem); |
|
533 if (v != _target && !_level->active(v)) { |
|
534 _level->activate(v); |
|
535 } |
|
536 } |
|
537 } |
|
538 return true; |
|
539 } |
|
540 |
|
541 /// \brief Starts the first phase of the preflow algorithm. |
|
542 /// |
|
543 /// The preflow algorithm consists of two phases, this method runs |
|
544 /// the first phase. After the first phase the maximum flow value |
|
545 /// and a minimum value cut can already be computed, although a |
|
546 /// maximum flow is not yet obtained. So after calling this method |
|
547 /// \ref flowValue() returns the value of a maximum flow and \ref |
|
548 /// minCut() returns a minimum cut. |
|
549 /// \pre One of the \ref init() functions should be called. |
|
550 void startFirstPhase() { |
|
551 _phase = true; |
|
552 |
|
553 Node n = _level->highestActive(); |
|
554 int level = _level->highestActiveLevel(); |
|
555 while (n != INVALID) { |
|
556 int num = _node_num; |
|
557 |
|
558 while (num > 0 && n != INVALID) { |
|
559 Value excess = (*_excess)[n]; |
|
560 int new_level = _level->maxLevel(); |
|
561 |
|
562 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
563 Value rem = (*_capacity)[e] - (*_flow)[e]; |
|
564 if (!_tolerance.positive(rem)) continue; |
|
565 Node v = _graph.target(e); |
|
566 if ((*_level)[v] < level) { |
|
567 if (!_level->active(v) && v != _target) { |
|
568 _level->activate(v); |
|
569 } |
|
570 if (!_tolerance.less(rem, excess)) { |
|
571 _flow->set(e, (*_flow)[e] + excess); |
|
572 _excess->set(v, (*_excess)[v] + excess); |
|
573 excess = 0; |
|
574 goto no_more_push_1; |
|
575 } else { |
|
576 excess -= rem; |
|
577 _excess->set(v, (*_excess)[v] + rem); |
|
578 _flow->set(e, (*_capacity)[e]); |
|
579 } |
|
580 } else if (new_level > (*_level)[v]) { |
|
581 new_level = (*_level)[v]; |
|
582 } |
|
583 } |
|
584 |
|
585 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
586 Value rem = (*_flow)[e]; |
|
587 if (!_tolerance.positive(rem)) continue; |
|
588 Node v = _graph.source(e); |
|
589 if ((*_level)[v] < level) { |
|
590 if (!_level->active(v) && v != _target) { |
|
591 _level->activate(v); |
|
592 } |
|
593 if (!_tolerance.less(rem, excess)) { |
|
594 _flow->set(e, (*_flow)[e] - excess); |
|
595 _excess->set(v, (*_excess)[v] + excess); |
|
596 excess = 0; |
|
597 goto no_more_push_1; |
|
598 } else { |
|
599 excess -= rem; |
|
600 _excess->set(v, (*_excess)[v] + rem); |
|
601 _flow->set(e, 0); |
|
602 } |
|
603 } else if (new_level > (*_level)[v]) { |
|
604 new_level = (*_level)[v]; |
|
605 } |
|
606 } |
|
607 |
|
608 no_more_push_1: |
|
609 |
|
610 _excess->set(n, excess); |
|
611 |
|
612 if (excess != 0) { |
|
613 if (new_level + 1 < _level->maxLevel()) { |
|
614 _level->liftHighestActive(new_level + 1); |
|
615 } else { |
|
616 _level->liftHighestActiveToTop(); |
|
617 } |
|
618 if (_level->emptyLevel(level)) { |
|
619 _level->liftToTop(level); |
|
620 } |
|
621 } else { |
|
622 _level->deactivate(n); |
|
623 } |
|
624 |
|
625 n = _level->highestActive(); |
|
626 level = _level->highestActiveLevel(); |
|
627 --num; |
|
628 } |
|
629 |
|
630 num = _node_num * 20; |
|
631 while (num > 0 && n != INVALID) { |
|
632 Value excess = (*_excess)[n]; |
|
633 int new_level = _level->maxLevel(); |
|
634 |
|
635 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
636 Value rem = (*_capacity)[e] - (*_flow)[e]; |
|
637 if (!_tolerance.positive(rem)) continue; |
|
638 Node v = _graph.target(e); |
|
639 if ((*_level)[v] < level) { |
|
640 if (!_level->active(v) && v != _target) { |
|
641 _level->activate(v); |
|
642 } |
|
643 if (!_tolerance.less(rem, excess)) { |
|
644 _flow->set(e, (*_flow)[e] + excess); |
|
645 _excess->set(v, (*_excess)[v] + excess); |
|
646 excess = 0; |
|
647 goto no_more_push_2; |
|
648 } else { |
|
649 excess -= rem; |
|
650 _excess->set(v, (*_excess)[v] + rem); |
|
651 _flow->set(e, (*_capacity)[e]); |
|
652 } |
|
653 } else if (new_level > (*_level)[v]) { |
|
654 new_level = (*_level)[v]; |
|
655 } |
|
656 } |
|
657 |
|
658 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
659 Value rem = (*_flow)[e]; |
|
660 if (!_tolerance.positive(rem)) continue; |
|
661 Node v = _graph.source(e); |
|
662 if ((*_level)[v] < level) { |
|
663 if (!_level->active(v) && v != _target) { |
|
664 _level->activate(v); |
|
665 } |
|
666 if (!_tolerance.less(rem, excess)) { |
|
667 _flow->set(e, (*_flow)[e] - excess); |
|
668 _excess->set(v, (*_excess)[v] + excess); |
|
669 excess = 0; |
|
670 goto no_more_push_2; |
|
671 } else { |
|
672 excess -= rem; |
|
673 _excess->set(v, (*_excess)[v] + rem); |
|
674 _flow->set(e, 0); |
|
675 } |
|
676 } else if (new_level > (*_level)[v]) { |
|
677 new_level = (*_level)[v]; |
|
678 } |
|
679 } |
|
680 |
|
681 no_more_push_2: |
|
682 |
|
683 _excess->set(n, excess); |
|
684 |
|
685 if (excess != 0) { |
|
686 if (new_level + 1 < _level->maxLevel()) { |
|
687 _level->liftActiveOn(level, new_level + 1); |
|
688 } else { |
|
689 _level->liftActiveToTop(level); |
|
690 } |
|
691 if (_level->emptyLevel(level)) { |
|
692 _level->liftToTop(level); |
|
693 } |
|
694 } else { |
|
695 _level->deactivate(n); |
|
696 } |
|
697 |
|
698 while (level >= 0 && _level->activeFree(level)) { |
|
699 --level; |
|
700 } |
|
701 if (level == -1) { |
|
702 n = _level->highestActive(); |
|
703 level = _level->highestActiveLevel(); |
|
704 } else { |
|
705 n = _level->activeOn(level); |
|
706 } |
|
707 --num; |
|
708 } |
|
709 } |
|
710 } |
|
711 |
|
712 /// \brief Starts the second phase of the preflow algorithm. |
|
713 /// |
|
714 /// The preflow algorithm consists of two phases, this method runs |
|
715 /// the second phase. After calling \ref init() and \ref |
|
716 /// startFirstPhase() and then \ref startSecondPhase(), \ref |
|
717 /// flowMap() return a maximum flow, \ref flowValue() returns the |
|
718 /// value of a maximum flow, \ref minCut() returns a minimum cut |
|
719 /// \pre The \ref init() and startFirstPhase() functions should be |
|
720 /// called before. |
|
721 void startSecondPhase() { |
|
722 _phase = false; |
|
723 |
|
724 typename Digraph::template NodeMap<bool> reached(_graph); |
|
725 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
726 reached.set(n, (*_level)[n] < _level->maxLevel()); |
|
727 } |
|
728 |
|
729 _level->initStart(); |
|
730 _level->initAddItem(_source); |
|
731 |
|
732 std::vector<Node> queue; |
|
733 queue.push_back(_source); |
|
734 reached.set(_source, true); |
|
735 |
|
736 while (!queue.empty()) { |
|
737 _level->initNewLevel(); |
|
738 std::vector<Node> nqueue; |
|
739 for (int i = 0; i < int(queue.size()); ++i) { |
|
740 Node n = queue[i]; |
|
741 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
742 Node v = _graph.target(e); |
|
743 if (!reached[v] && _tolerance.positive((*_flow)[e])) { |
|
744 reached.set(v, true); |
|
745 _level->initAddItem(v); |
|
746 nqueue.push_back(v); |
|
747 } |
|
748 } |
|
749 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
750 Node u = _graph.source(e); |
|
751 if (!reached[u] && |
|
752 _tolerance.positive((*_capacity)[e] - (*_flow)[e])) { |
|
753 reached.set(u, true); |
|
754 _level->initAddItem(u); |
|
755 nqueue.push_back(u); |
|
756 } |
|
757 } |
|
758 } |
|
759 queue.swap(nqueue); |
|
760 } |
|
761 _level->initFinish(); |
|
762 |
|
763 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
764 if (!reached[n]) { |
|
765 _level->dirtyTopButOne(n); |
|
766 } else if ((*_excess)[n] > 0 && _target != n) { |
|
767 _level->activate(n); |
|
768 } |
|
769 } |
|
770 |
|
771 Node n; |
|
772 while ((n = _level->highestActive()) != INVALID) { |
|
773 Value excess = (*_excess)[n]; |
|
774 int level = _level->highestActiveLevel(); |
|
775 int new_level = _level->maxLevel(); |
|
776 |
|
777 for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
|
778 Value rem = (*_capacity)[e] - (*_flow)[e]; |
|
779 if (!_tolerance.positive(rem)) continue; |
|
780 Node v = _graph.target(e); |
|
781 if ((*_level)[v] < level) { |
|
782 if (!_level->active(v) && v != _source) { |
|
783 _level->activate(v); |
|
784 } |
|
785 if (!_tolerance.less(rem, excess)) { |
|
786 _flow->set(e, (*_flow)[e] + excess); |
|
787 _excess->set(v, (*_excess)[v] + excess); |
|
788 excess = 0; |
|
789 goto no_more_push; |
|
790 } else { |
|
791 excess -= rem; |
|
792 _excess->set(v, (*_excess)[v] + rem); |
|
793 _flow->set(e, (*_capacity)[e]); |
|
794 } |
|
795 } else if (new_level > (*_level)[v]) { |
|
796 new_level = (*_level)[v]; |
|
797 } |
|
798 } |
|
799 |
|
800 for (InArcIt e(_graph, n); e != INVALID; ++e) { |
|
801 Value rem = (*_flow)[e]; |
|
802 if (!_tolerance.positive(rem)) continue; |
|
803 Node v = _graph.source(e); |
|
804 if ((*_level)[v] < level) { |
|
805 if (!_level->active(v) && v != _source) { |
|
806 _level->activate(v); |
|
807 } |
|
808 if (!_tolerance.less(rem, excess)) { |
|
809 _flow->set(e, (*_flow)[e] - excess); |
|
810 _excess->set(v, (*_excess)[v] + excess); |
|
811 excess = 0; |
|
812 goto no_more_push; |
|
813 } else { |
|
814 excess -= rem; |
|
815 _excess->set(v, (*_excess)[v] + rem); |
|
816 _flow->set(e, 0); |
|
817 } |
|
818 } else if (new_level > (*_level)[v]) { |
|
819 new_level = (*_level)[v]; |
|
820 } |
|
821 } |
|
822 |
|
823 no_more_push: |
|
824 |
|
825 _excess->set(n, excess); |
|
826 |
|
827 if (excess != 0) { |
|
828 if (new_level + 1 < _level->maxLevel()) { |
|
829 _level->liftHighestActive(new_level + 1); |
|
830 } else { |
|
831 // Calculation error |
|
832 _level->liftHighestActiveToTop(); |
|
833 } |
|
834 if (_level->emptyLevel(level)) { |
|
835 // Calculation error |
|
836 _level->liftToTop(level); |
|
837 } |
|
838 } else { |
|
839 _level->deactivate(n); |
|
840 } |
|
841 |
|
842 } |
|
843 } |
|
844 |
|
845 /// \brief Runs the preflow algorithm. |
|
846 /// |
|
847 /// Runs the preflow algorithm. |
|
848 /// \note pf.run() is just a shortcut of the following code. |
|
849 /// \code |
|
850 /// pf.init(); |
|
851 /// pf.startFirstPhase(); |
|
852 /// pf.startSecondPhase(); |
|
853 /// \endcode |
|
854 void run() { |
|
855 init(); |
|
856 startFirstPhase(); |
|
857 startSecondPhase(); |
|
858 } |
|
859 |
|
860 /// \brief Runs the preflow algorithm to compute the minimum cut. |
|
861 /// |
|
862 /// Runs the preflow algorithm to compute the minimum cut. |
|
863 /// \note pf.runMinCut() is just a shortcut of the following code. |
|
864 /// \code |
|
865 /// pf.init(); |
|
866 /// pf.startFirstPhase(); |
|
867 /// \endcode |
|
868 void runMinCut() { |
|
869 init(); |
|
870 startFirstPhase(); |
|
871 } |
|
872 |
|
873 /// @} |
|
874 |
|
875 /// \name Query Functions |
|
876 /// The result of the %Preflow algorithm can be obtained using these |
|
877 /// functions.\n |
|
878 /// Before the use of these functions, |
|
879 /// either run() or start() must be called. |
|
880 |
|
881 ///@{ |
|
882 |
|
883 /// \brief Returns the value of the maximum flow. |
|
884 /// |
|
885 /// Returns the value of the maximum flow by returning the excess |
|
886 /// of the target node \c t. This value equals to the value of |
|
887 /// the maximum flow already after the first phase. |
|
888 Value flowValue() const { |
|
889 return (*_excess)[_target]; |
|
890 } |
|
891 |
|
892 /// \brief Returns true when the node is on the source side of minimum cut. |
|
893 /// |
|
894 /// Returns true when the node is on the source side of minimum |
|
895 /// cut. This method can be called both after running \ref |
|
896 /// startFirstPhase() and \ref startSecondPhase(). |
|
897 bool minCut(const Node& node) const { |
|
898 return ((*_level)[node] == _level->maxLevel()) == _phase; |
|
899 } |
|
900 |
|
901 /// \brief Returns a minimum value cut. |
|
902 /// |
|
903 /// Sets the \c cutMap to the characteristic vector of a minimum value |
|
904 /// cut. This method can be called both after running \ref |
|
905 /// startFirstPhase() and \ref startSecondPhase(). The result after second |
|
906 /// phase could be changed slightly if inexact computation is used. |
|
907 /// \pre The \c cutMap should be a bool-valued node-map. |
|
908 template <typename CutMap> |
|
909 void minCutMap(CutMap& cutMap) const { |
|
910 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
911 cutMap.set(n, minCut(n)); |
|
912 } |
|
913 } |
|
914 |
|
915 /// \brief Returns the flow on the arc. |
|
916 /// |
|
917 /// Sets the \c flowMap to the flow on the arcs. This method can |
|
918 /// be called after the second phase of algorithm. |
|
919 Value flow(const Arc& arc) const { |
|
920 return (*_flow)[arc]; |
|
921 } |
|
922 |
|
923 /// @} |
|
924 }; |
|
925 } |
|
926 |
|
927 #endif |