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