1 /* -*- C++ -*- |
|
2 * src/lemon/preflow.h - Part of LEMON, a generic C++ optimization library |
|
3 * |
|
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
|
5 * (Egervary Research Group on Combinatorial Optimization, EGRES). |
|
6 * |
|
7 * Permission to use, modify and distribute this software is granted |
|
8 * provided that this copyright notice appears in all copies. For |
|
9 * precise terms see the accompanying LICENSE file. |
|
10 * |
|
11 * This software is provided "AS IS" with no warranty of any kind, |
|
12 * express or implied, and with no claim as to its suitability for any |
|
13 * purpose. |
|
14 * |
|
15 */ |
|
16 |
|
17 #ifndef LEMON_PREFLOW_H |
|
18 #define LEMON_PREFLOW_H |
|
19 |
|
20 #include <vector> |
|
21 #include <queue> |
|
22 |
|
23 #include <lemon/invalid.h> |
|
24 #include <lemon/maps.h> |
|
25 #include <lemon/graph_utils.h> |
|
26 |
|
27 /// \file |
|
28 /// \ingroup flowalgs |
|
29 /// Implementation of the preflow algorithm. |
|
30 |
|
31 namespace lemon { |
|
32 |
|
33 /// \addtogroup flowalgs |
|
34 /// @{ |
|
35 |
|
36 ///%Preflow algorithms class. |
|
37 |
|
38 ///This class provides an implementation of the \e preflow \e |
|
39 ///algorithm producing a flow of maximum value in a directed |
|
40 ///graph. The preflow algorithms are the fastest known max flow algorithms |
|
41 ///up to now. The \e source node, the \e target node, the \e |
|
42 ///capacity of the edges and the \e starting \e flow value of the |
|
43 ///edges should be passed to the algorithm through the |
|
44 ///constructor. It is possible to change these quantities using the |
|
45 ///functions \ref source, \ref target, \ref capacityMap and \ref |
|
46 ///flowMap. |
|
47 /// |
|
48 ///After running \ref lemon::Preflow::phase1() "phase1()" |
|
49 ///or \ref lemon::Preflow::run() "run()", the maximal flow |
|
50 ///value can be obtained by calling \ref flowValue(). The minimum |
|
51 ///value cut can be written into a <tt>bool</tt> node map by |
|
52 ///calling \ref minCut(). (\ref minMinCut() and \ref maxMinCut() writes |
|
53 ///the inclusionwise minimum and maximum of the minimum value cuts, |
|
54 ///resp.) |
|
55 /// |
|
56 ///\param Graph The directed graph type the algorithm runs on. |
|
57 ///\param Num The number type of the capacities and the flow values. |
|
58 ///\param CapacityMap The capacity map type. |
|
59 ///\param FlowMap The flow map type. |
|
60 /// |
|
61 ///\author Jacint Szabo |
|
62 ///\todo Second template parameter is superfluous |
|
63 template <typename Graph, typename Num, |
|
64 typename CapacityMap=typename Graph::template EdgeMap<Num>, |
|
65 typename FlowMap=typename Graph::template EdgeMap<Num> > |
|
66 class Preflow { |
|
67 protected: |
|
68 typedef typename Graph::Node Node; |
|
69 typedef typename Graph::NodeIt NodeIt; |
|
70 typedef typename Graph::EdgeIt EdgeIt; |
|
71 typedef typename Graph::OutEdgeIt OutEdgeIt; |
|
72 typedef typename Graph::InEdgeIt InEdgeIt; |
|
73 |
|
74 typedef typename Graph::template NodeMap<Node> NNMap; |
|
75 typedef typename std::vector<Node> VecNode; |
|
76 |
|
77 const Graph* _g; |
|
78 Node _source; |
|
79 Node _target; |
|
80 const CapacityMap* _capacity; |
|
81 FlowMap* _flow; |
|
82 int _node_num; //the number of nodes of G |
|
83 |
|
84 typename Graph::template NodeMap<int> level; |
|
85 typename Graph::template NodeMap<Num> excess; |
|
86 |
|
87 // constants used for heuristics |
|
88 static const int H0=20; |
|
89 static const int H1=1; |
|
90 |
|
91 public: |
|
92 |
|
93 ///Indicates the property of the starting flow map. |
|
94 |
|
95 ///Indicates the property of the starting flow map. |
|
96 ///The meanings are as follows: |
|
97 ///- \c ZERO_FLOW: constant zero flow |
|
98 ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to |
|
99 ///the sum of the out-flows in every node except the \e source and |
|
100 ///the \e target. |
|
101 ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at |
|
102 ///least the sum of the out-flows in every node except the \e source. |
|
103 ///- \c NO_FLOW: indicates an unspecified edge map. \c flow will be |
|
104 ///set to the constant zero flow in the beginning of |
|
105 ///the algorithm in this case. |
|
106 /// |
|
107 enum FlowEnum{ |
|
108 NO_FLOW, |
|
109 ZERO_FLOW, |
|
110 GEN_FLOW, |
|
111 PRE_FLOW |
|
112 }; |
|
113 |
|
114 ///Indicates the state of the preflow algorithm. |
|
115 |
|
116 ///Indicates the state of the preflow algorithm. |
|
117 ///The meanings are as follows: |
|
118 ///- \c AFTER_NOTHING: before running the algorithm or |
|
119 /// at an unspecified state. |
|
120 ///- \c AFTER_PREFLOW_PHASE_1: right after running \c phase1 |
|
121 ///- \c AFTER_PREFLOW_PHASE_2: after running \ref phase2() |
|
122 /// |
|
123 enum StatusEnum { |
|
124 AFTER_NOTHING, |
|
125 AFTER_PREFLOW_PHASE_1, |
|
126 AFTER_PREFLOW_PHASE_2 |
|
127 }; |
|
128 |
|
129 protected: |
|
130 FlowEnum flow_prop; |
|
131 StatusEnum status; // Do not needle this flag only if necessary. |
|
132 |
|
133 public: |
|
134 ///The constructor of the class. |
|
135 |
|
136 ///The constructor of the class. |
|
137 ///\param _gr The directed graph the algorithm runs on. |
|
138 ///\param _s The source node. |
|
139 ///\param _t The target node. |
|
140 ///\param _cap The capacity of the edges. |
|
141 ///\param _f The flow of the edges. |
|
142 ///Except the graph, all of these parameters can be reset by |
|
143 ///calling \ref source, \ref target, \ref capacityMap and \ref |
|
144 ///flowMap, resp. |
|
145 Preflow(const Graph& _gr, Node _s, Node _t, |
|
146 const CapacityMap& _cap, FlowMap& _f) : |
|
147 _g(&_gr), _source(_s), _target(_t), _capacity(&_cap), |
|
148 _flow(&_f), _node_num(countNodes(_gr)), level(_gr), excess(_gr,0), |
|
149 flow_prop(NO_FLOW), status(AFTER_NOTHING) { } |
|
150 |
|
151 |
|
152 |
|
153 ///Runs the preflow algorithm. |
|
154 |
|
155 ///Runs the preflow algorithm. |
|
156 /// |
|
157 void run() { |
|
158 phase1(flow_prop); |
|
159 phase2(); |
|
160 } |
|
161 |
|
162 ///Runs the preflow algorithm. |
|
163 |
|
164 ///Runs the preflow algorithm. |
|
165 ///\pre The starting flow map must be |
|
166 /// - a constant zero flow if \c fp is \c ZERO_FLOW, |
|
167 /// - an arbitrary flow if \c fp is \c GEN_FLOW, |
|
168 /// - an arbitrary preflow if \c fp is \c PRE_FLOW, |
|
169 /// - any map if \c fp is NO_FLOW. |
|
170 ///If the starting flow map is a flow or a preflow then |
|
171 ///the algorithm terminates faster. |
|
172 void run(FlowEnum fp) { |
|
173 flow_prop=fp; |
|
174 run(); |
|
175 } |
|
176 |
|
177 ///Runs the first phase of the preflow algorithm. |
|
178 |
|
179 ///The preflow algorithm consists of two phases, this method runs |
|
180 ///the first phase. After the first phase the maximum flow value |
|
181 ///and a minimum value cut can already be computed, although a |
|
182 ///maximum flow is not yet obtained. So after calling this method |
|
183 ///\ref flowValue returns the value of a maximum flow and \ref |
|
184 ///minCut returns a minimum cut. |
|
185 ///\warning \ref minMinCut and \ref maxMinCut do not give minimum |
|
186 ///value cuts unless calling \ref phase2. |
|
187 ///\pre The starting flow must be |
|
188 ///- a constant zero flow if \c fp is \c ZERO_FLOW, |
|
189 ///- an arbitary flow if \c fp is \c GEN_FLOW, |
|
190 ///- an arbitary preflow if \c fp is \c PRE_FLOW, |
|
191 ///- any map if \c fp is NO_FLOW. |
|
192 void phase1(FlowEnum fp) |
|
193 { |
|
194 flow_prop=fp; |
|
195 phase1(); |
|
196 } |
|
197 |
|
198 |
|
199 ///Runs the first phase of the preflow algorithm. |
|
200 |
|
201 ///The preflow algorithm consists of two phases, this method runs |
|
202 ///the first phase. After the first phase the maximum flow value |
|
203 ///and a minimum value cut can already be computed, although a |
|
204 ///maximum flow is not yet obtained. So after calling this method |
|
205 ///\ref flowValue returns the value of a maximum flow and \ref |
|
206 ///minCut returns a minimum cut. |
|
207 ///\warning \ref minCut(), \ref minMinCut() and \ref maxMinCut() do not |
|
208 ///give minimum value cuts unless calling \ref phase2(). |
|
209 void phase1() |
|
210 { |
|
211 int heur0=(int)(H0*_node_num); //time while running 'bound decrease' |
|
212 int heur1=(int)(H1*_node_num); //time while running 'highest label' |
|
213 int heur=heur1; //starting time interval (#of relabels) |
|
214 int numrelabel=0; |
|
215 |
|
216 bool what_heur=1; |
|
217 //It is 0 in case 'bound decrease' and 1 in case 'highest label' |
|
218 |
|
219 bool end=false; |
|
220 //Needed for 'bound decrease', true means no active |
|
221 //nodes are above bound b. |
|
222 |
|
223 int k=_node_num-2; //bound on the highest level under n containing a node |
|
224 int b=k; //bound on the highest level under n of an active node |
|
225 |
|
226 VecNode first(_node_num, INVALID); |
|
227 NNMap next(*_g, INVALID); |
|
228 |
|
229 NNMap left(*_g, INVALID); |
|
230 NNMap right(*_g, INVALID); |
|
231 VecNode level_list(_node_num,INVALID); |
|
232 //List of the nodes in level i<n, set to n. |
|
233 |
|
234 preflowPreproc(first, next, level_list, left, right); |
|
235 |
|
236 //Push/relabel on the highest level active nodes. |
|
237 while ( true ) { |
|
238 if ( b == 0 ) { |
|
239 if ( !what_heur && !end && k > 0 ) { |
|
240 b=k; |
|
241 end=true; |
|
242 } else break; |
|
243 } |
|
244 |
|
245 if ( first[b]==INVALID ) --b; |
|
246 else { |
|
247 end=false; |
|
248 Node w=first[b]; |
|
249 first[b]=next[w]; |
|
250 int newlevel=push(w, next, first); |
|
251 if ( excess[w] > 0 ) relabel(w, newlevel, first, next, level_list, |
|
252 left, right, b, k, what_heur); |
|
253 |
|
254 ++numrelabel; |
|
255 if ( numrelabel >= heur ) { |
|
256 numrelabel=0; |
|
257 if ( what_heur ) { |
|
258 what_heur=0; |
|
259 heur=heur0; |
|
260 end=false; |
|
261 } else { |
|
262 what_heur=1; |
|
263 heur=heur1; |
|
264 b=k; |
|
265 } |
|
266 } |
|
267 } |
|
268 } |
|
269 flow_prop=PRE_FLOW; |
|
270 status=AFTER_PREFLOW_PHASE_1; |
|
271 } |
|
272 // Heuristics: |
|
273 // 2 phase |
|
274 // gap |
|
275 // list 'level_list' on the nodes on level i implemented by hand |
|
276 // stack 'active' on the active nodes on level i |
|
277 // runs heuristic 'highest label' for H1*n relabels |
|
278 // runs heuristic 'bound decrease' for H0*n relabels, |
|
279 // starts with 'highest label' |
|
280 // Parameters H0 and H1 are initialized to 20 and 1. |
|
281 |
|
282 |
|
283 ///Runs the second phase of the preflow algorithm. |
|
284 |
|
285 ///The preflow algorithm consists of two phases, this method runs |
|
286 ///the second phase. After calling \ref phase1 and then \ref |
|
287 ///phase2, \ref flow contains a maximum flow, \ref flowValue |
|
288 ///returns the value of a maximum flow, \ref minCut returns a |
|
289 ///minimum cut, while the methods \ref minMinCut and \ref |
|
290 ///maxMinCut return the inclusionwise minimum and maximum cuts of |
|
291 ///minimum value, resp. \pre \ref phase1 must be called before. |
|
292 void phase2() |
|
293 { |
|
294 |
|
295 int k=_node_num-2; //bound on the highest level under n containing a node |
|
296 int b=k; //bound on the highest level under n of an active node |
|
297 |
|
298 |
|
299 VecNode first(_node_num, INVALID); |
|
300 NNMap next(*_g, INVALID); |
|
301 level.set(_source,0); |
|
302 std::queue<Node> bfs_queue; |
|
303 bfs_queue.push(_source); |
|
304 |
|
305 while ( !bfs_queue.empty() ) { |
|
306 |
|
307 Node v=bfs_queue.front(); |
|
308 bfs_queue.pop(); |
|
309 int l=level[v]+1; |
|
310 |
|
311 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) { |
|
312 if ( (*_capacity)[e] <= (*_flow)[e] ) continue; |
|
313 Node u=_g->source(e); |
|
314 if ( level[u] >= _node_num ) { |
|
315 bfs_queue.push(u); |
|
316 level.set(u, l); |
|
317 if ( excess[u] > 0 ) { |
|
318 next.set(u,first[l]); |
|
319 first[l]=u; |
|
320 } |
|
321 } |
|
322 } |
|
323 |
|
324 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) { |
|
325 if ( 0 >= (*_flow)[e] ) continue; |
|
326 Node u=_g->target(e); |
|
327 if ( level[u] >= _node_num ) { |
|
328 bfs_queue.push(u); |
|
329 level.set(u, l); |
|
330 if ( excess[u] > 0 ) { |
|
331 next.set(u,first[l]); |
|
332 first[l]=u; |
|
333 } |
|
334 } |
|
335 } |
|
336 } |
|
337 b=_node_num-2; |
|
338 |
|
339 while ( true ) { |
|
340 |
|
341 if ( b == 0 ) break; |
|
342 if ( first[b]==INVALID ) --b; |
|
343 else { |
|
344 Node w=first[b]; |
|
345 first[b]=next[w]; |
|
346 int newlevel=push(w,next, first); |
|
347 |
|
348 //relabel |
|
349 if ( excess[w] > 0 ) { |
|
350 level.set(w,++newlevel); |
|
351 next.set(w,first[newlevel]); |
|
352 first[newlevel]=w; |
|
353 b=newlevel; |
|
354 } |
|
355 } |
|
356 } // while(true) |
|
357 flow_prop=GEN_FLOW; |
|
358 status=AFTER_PREFLOW_PHASE_2; |
|
359 } |
|
360 |
|
361 /// Returns the value of the maximum flow. |
|
362 |
|
363 /// Returns the value of the maximum flow by returning the excess |
|
364 /// of the target node \c t. This value equals to the value of |
|
365 /// the maximum flow already after running \ref phase1. |
|
366 Num flowValue() const { |
|
367 return excess[_target]; |
|
368 } |
|
369 |
|
370 |
|
371 ///Returns a minimum value cut. |
|
372 |
|
373 ///Sets \c M to the characteristic vector of a minimum value |
|
374 ///cut. This method can be called both after running \ref |
|
375 ///phase1 and \ref phase2. It is much faster after |
|
376 ///\ref phase1. \pre M should be a bool-valued node-map. \pre |
|
377 ///If \ref minCut() is called after \ref phase2() then M should |
|
378 ///be initialized to false. |
|
379 template<typename _CutMap> |
|
380 void minCut(_CutMap& M) const { |
|
381 switch ( status ) { |
|
382 case AFTER_PREFLOW_PHASE_1: |
|
383 for(NodeIt v(*_g); v!=INVALID; ++v) { |
|
384 if (level[v] < _node_num) { |
|
385 M.set(v, false); |
|
386 } else { |
|
387 M.set(v, true); |
|
388 } |
|
389 } |
|
390 break; |
|
391 case AFTER_PREFLOW_PHASE_2: |
|
392 minMinCut(M); |
|
393 break; |
|
394 case AFTER_NOTHING: |
|
395 break; |
|
396 } |
|
397 } |
|
398 |
|
399 ///Returns the inclusionwise minimum of the minimum value cuts. |
|
400 |
|
401 ///Sets \c M to the characteristic vector of the minimum value cut |
|
402 ///which is inclusionwise minimum. It is computed by processing a |
|
403 ///bfs from the source node \c s in the residual graph. \pre M |
|
404 ///should be a node map of bools initialized to false. \pre \ref |
|
405 ///phase2 should already be run. |
|
406 template<typename _CutMap> |
|
407 void minMinCut(_CutMap& M) const { |
|
408 |
|
409 std::queue<Node> queue; |
|
410 M.set(_source,true); |
|
411 queue.push(_source); |
|
412 |
|
413 while (!queue.empty()) { |
|
414 Node w=queue.front(); |
|
415 queue.pop(); |
|
416 |
|
417 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
|
418 Node v=_g->target(e); |
|
419 if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) { |
|
420 queue.push(v); |
|
421 M.set(v, true); |
|
422 } |
|
423 } |
|
424 |
|
425 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
|
426 Node v=_g->source(e); |
|
427 if (!M[v] && (*_flow)[e] > 0 ) { |
|
428 queue.push(v); |
|
429 M.set(v, true); |
|
430 } |
|
431 } |
|
432 } |
|
433 } |
|
434 |
|
435 ///Returns the inclusionwise maximum of the minimum value cuts. |
|
436 |
|
437 ///Sets \c M to the characteristic vector of the minimum value cut |
|
438 ///which is inclusionwise maximum. It is computed by processing a |
|
439 ///backward bfs from the target node \c t in the residual graph. |
|
440 ///\pre \ref phase2() or run() should already be run. |
|
441 template<typename _CutMap> |
|
442 void maxMinCut(_CutMap& M) const { |
|
443 |
|
444 for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true); |
|
445 |
|
446 std::queue<Node> queue; |
|
447 |
|
448 M.set(_target,false); |
|
449 queue.push(_target); |
|
450 |
|
451 while (!queue.empty()) { |
|
452 Node w=queue.front(); |
|
453 queue.pop(); |
|
454 |
|
455 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
|
456 Node v=_g->source(e); |
|
457 if (M[v] && (*_flow)[e] < (*_capacity)[e] ) { |
|
458 queue.push(v); |
|
459 M.set(v, false); |
|
460 } |
|
461 } |
|
462 |
|
463 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
|
464 Node v=_g->target(e); |
|
465 if (M[v] && (*_flow)[e] > 0 ) { |
|
466 queue.push(v); |
|
467 M.set(v, false); |
|
468 } |
|
469 } |
|
470 } |
|
471 } |
|
472 |
|
473 ///Sets the source node to \c _s. |
|
474 |
|
475 ///Sets the source node to \c _s. |
|
476 /// |
|
477 void source(Node _s) { |
|
478 _source=_s; |
|
479 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW; |
|
480 status=AFTER_NOTHING; |
|
481 } |
|
482 |
|
483 ///Returns the source node. |
|
484 |
|
485 ///Returns the source node. |
|
486 /// |
|
487 Node source() const { |
|
488 return _source; |
|
489 } |
|
490 |
|
491 ///Sets the target node to \c _t. |
|
492 |
|
493 ///Sets the target node to \c _t. |
|
494 /// |
|
495 void target(Node _t) { |
|
496 _target=_t; |
|
497 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW; |
|
498 status=AFTER_NOTHING; |
|
499 } |
|
500 |
|
501 ///Returns the target node. |
|
502 |
|
503 ///Returns the target node. |
|
504 /// |
|
505 Node target() const { |
|
506 return _target; |
|
507 } |
|
508 |
|
509 /// Sets the edge map of the capacities to _cap. |
|
510 |
|
511 /// Sets the edge map of the capacities to _cap. |
|
512 /// |
|
513 void capacityMap(const CapacityMap& _cap) { |
|
514 _capacity=&_cap; |
|
515 status=AFTER_NOTHING; |
|
516 } |
|
517 /// Returns a reference to capacity map. |
|
518 |
|
519 /// Returns a reference to capacity map. |
|
520 /// |
|
521 const CapacityMap &capacityMap() const { |
|
522 return *_capacity; |
|
523 } |
|
524 |
|
525 /// Sets the edge map of the flows to _flow. |
|
526 |
|
527 /// Sets the edge map of the flows to _flow. |
|
528 /// |
|
529 void flowMap(FlowMap& _f) { |
|
530 _flow=&_f; |
|
531 flow_prop=NO_FLOW; |
|
532 status=AFTER_NOTHING; |
|
533 } |
|
534 |
|
535 /// Returns a reference to flow map. |
|
536 |
|
537 /// Returns a reference to flow map. |
|
538 /// |
|
539 const FlowMap &flowMap() const { |
|
540 return *_flow; |
|
541 } |
|
542 |
|
543 private: |
|
544 |
|
545 int push(Node w, NNMap& next, VecNode& first) { |
|
546 |
|
547 int lev=level[w]; |
|
548 Num exc=excess[w]; |
|
549 int newlevel=_node_num; //bound on the next level of w |
|
550 |
|
551 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
|
552 if ( (*_flow)[e] >= (*_capacity)[e] ) continue; |
|
553 Node v=_g->target(e); |
|
554 |
|
555 if( lev > level[v] ) { //Push is allowed now |
|
556 |
|
557 if ( excess[v]<=0 && v!=_target && v!=_source ) { |
|
558 next.set(v,first[level[v]]); |
|
559 first[level[v]]=v; |
|
560 } |
|
561 |
|
562 Num cap=(*_capacity)[e]; |
|
563 Num flo=(*_flow)[e]; |
|
564 Num remcap=cap-flo; |
|
565 |
|
566 if ( remcap >= exc ) { //A nonsaturating push. |
|
567 |
|
568 _flow->set(e, flo+exc); |
|
569 excess.set(v, excess[v]+exc); |
|
570 exc=0; |
|
571 break; |
|
572 |
|
573 } else { //A saturating push. |
|
574 _flow->set(e, cap); |
|
575 excess.set(v, excess[v]+remcap); |
|
576 exc-=remcap; |
|
577 } |
|
578 } else if ( newlevel > level[v] ) newlevel = level[v]; |
|
579 } //for out edges wv |
|
580 |
|
581 if ( exc > 0 ) { |
|
582 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
|
583 |
|
584 if( (*_flow)[e] <= 0 ) continue; |
|
585 Node v=_g->source(e); |
|
586 |
|
587 if( lev > level[v] ) { //Push is allowed now |
|
588 |
|
589 if ( excess[v]<=0 && v!=_target && v!=_source ) { |
|
590 next.set(v,first[level[v]]); |
|
591 first[level[v]]=v; |
|
592 } |
|
593 |
|
594 Num flo=(*_flow)[e]; |
|
595 |
|
596 if ( flo >= exc ) { //A nonsaturating push. |
|
597 |
|
598 _flow->set(e, flo-exc); |
|
599 excess.set(v, excess[v]+exc); |
|
600 exc=0; |
|
601 break; |
|
602 } else { //A saturating push. |
|
603 |
|
604 excess.set(v, excess[v]+flo); |
|
605 exc-=flo; |
|
606 _flow->set(e,0); |
|
607 } |
|
608 } else if ( newlevel > level[v] ) newlevel = level[v]; |
|
609 } //for in edges vw |
|
610 |
|
611 } // if w still has excess after the out edge for cycle |
|
612 |
|
613 excess.set(w, exc); |
|
614 |
|
615 return newlevel; |
|
616 } |
|
617 |
|
618 |
|
619 |
|
620 void preflowPreproc(VecNode& first, NNMap& next, |
|
621 VecNode& level_list, NNMap& left, NNMap& right) |
|
622 { |
|
623 for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num); |
|
624 std::queue<Node> bfs_queue; |
|
625 |
|
626 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) { |
|
627 //Reverse_bfs from t in the residual graph, |
|
628 //to find the starting level. |
|
629 level.set(_target,0); |
|
630 bfs_queue.push(_target); |
|
631 |
|
632 while ( !bfs_queue.empty() ) { |
|
633 |
|
634 Node v=bfs_queue.front(); |
|
635 bfs_queue.pop(); |
|
636 int l=level[v]+1; |
|
637 |
|
638 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
|
639 if ( (*_capacity)[e] <= (*_flow)[e] ) continue; |
|
640 Node w=_g->source(e); |
|
641 if ( level[w] == _node_num && w != _source ) { |
|
642 bfs_queue.push(w); |
|
643 Node z=level_list[l]; |
|
644 if ( z!=INVALID ) left.set(z,w); |
|
645 right.set(w,z); |
|
646 level_list[l]=w; |
|
647 level.set(w, l); |
|
648 } |
|
649 } |
|
650 |
|
651 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
|
652 if ( 0 >= (*_flow)[e] ) continue; |
|
653 Node w=_g->target(e); |
|
654 if ( level[w] == _node_num && w != _source ) { |
|
655 bfs_queue.push(w); |
|
656 Node z=level_list[l]; |
|
657 if ( z!=INVALID ) left.set(z,w); |
|
658 right.set(w,z); |
|
659 level_list[l]=w; |
|
660 level.set(w, l); |
|
661 } |
|
662 } |
|
663 } //while |
|
664 } //if |
|
665 |
|
666 |
|
667 switch (flow_prop) { |
|
668 case NO_FLOW: |
|
669 for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0); |
|
670 case ZERO_FLOW: |
|
671 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0); |
|
672 |
|
673 //Reverse_bfs from t, to find the starting level. |
|
674 level.set(_target,0); |
|
675 bfs_queue.push(_target); |
|
676 |
|
677 while ( !bfs_queue.empty() ) { |
|
678 |
|
679 Node v=bfs_queue.front(); |
|
680 bfs_queue.pop(); |
|
681 int l=level[v]+1; |
|
682 |
|
683 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
|
684 Node w=_g->source(e); |
|
685 if ( level[w] == _node_num && w != _source ) { |
|
686 bfs_queue.push(w); |
|
687 Node z=level_list[l]; |
|
688 if ( z!=INVALID ) left.set(z,w); |
|
689 right.set(w,z); |
|
690 level_list[l]=w; |
|
691 level.set(w, l); |
|
692 } |
|
693 } |
|
694 } |
|
695 |
|
696 //the starting flow |
|
697 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
|
698 Num c=(*_capacity)[e]; |
|
699 if ( c <= 0 ) continue; |
|
700 Node w=_g->target(e); |
|
701 if ( level[w] < _node_num ) { |
|
702 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack |
|
703 next.set(w,first[level[w]]); |
|
704 first[level[w]]=w; |
|
705 } |
|
706 _flow->set(e, c); |
|
707 excess.set(w, excess[w]+c); |
|
708 } |
|
709 } |
|
710 break; |
|
711 |
|
712 case GEN_FLOW: |
|
713 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0); |
|
714 { |
|
715 Num exc=0; |
|
716 for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e]; |
|
717 for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e]; |
|
718 excess.set(_target,exc); |
|
719 } |
|
720 |
|
721 //the starting flow |
|
722 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) { |
|
723 Num rem=(*_capacity)[e]-(*_flow)[e]; |
|
724 if ( rem <= 0 ) continue; |
|
725 Node w=_g->target(e); |
|
726 if ( level[w] < _node_num ) { |
|
727 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack |
|
728 next.set(w,first[level[w]]); |
|
729 first[level[w]]=w; |
|
730 } |
|
731 _flow->set(e, (*_capacity)[e]); |
|
732 excess.set(w, excess[w]+rem); |
|
733 } |
|
734 } |
|
735 |
|
736 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) { |
|
737 if ( (*_flow)[e] <= 0 ) continue; |
|
738 Node w=_g->source(e); |
|
739 if ( level[w] < _node_num ) { |
|
740 if ( excess[w] <= 0 && w!=_target ) { |
|
741 next.set(w,first[level[w]]); |
|
742 first[level[w]]=w; |
|
743 } |
|
744 excess.set(w, excess[w]+(*_flow)[e]); |
|
745 _flow->set(e, 0); |
|
746 } |
|
747 } |
|
748 break; |
|
749 |
|
750 case PRE_FLOW: |
|
751 //the starting flow |
|
752 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
|
753 Num rem=(*_capacity)[e]-(*_flow)[e]; |
|
754 if ( rem <= 0 ) continue; |
|
755 Node w=_g->target(e); |
|
756 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]); |
|
757 } |
|
758 |
|
759 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
|
760 if ( (*_flow)[e] <= 0 ) continue; |
|
761 Node w=_g->source(e); |
|
762 if ( level[w] < _node_num ) _flow->set(e, 0); |
|
763 } |
|
764 |
|
765 //computing the excess |
|
766 for(NodeIt w(*_g); w!=INVALID; ++w) { |
|
767 Num exc=0; |
|
768 for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e]; |
|
769 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e]; |
|
770 excess.set(w,exc); |
|
771 |
|
772 //putting the active nodes into the stack |
|
773 int lev=level[w]; |
|
774 if ( exc > 0 && lev < _node_num && Node(w) != _target ) { |
|
775 next.set(w,first[lev]); |
|
776 first[lev]=w; |
|
777 } |
|
778 } |
|
779 break; |
|
780 } //switch |
|
781 } //preflowPreproc |
|
782 |
|
783 |
|
784 void relabel(Node w, int newlevel, VecNode& first, NNMap& next, |
|
785 VecNode& level_list, NNMap& left, |
|
786 NNMap& right, int& b, int& k, bool what_heur ) |
|
787 { |
|
788 |
|
789 int lev=level[w]; |
|
790 |
|
791 Node right_n=right[w]; |
|
792 Node left_n=left[w]; |
|
793 |
|
794 //unlacing starts |
|
795 if ( right_n!=INVALID ) { |
|
796 if ( left_n!=INVALID ) { |
|
797 right.set(left_n, right_n); |
|
798 left.set(right_n, left_n); |
|
799 } else { |
|
800 level_list[lev]=right_n; |
|
801 left.set(right_n, INVALID); |
|
802 } |
|
803 } else { |
|
804 if ( left_n!=INVALID ) { |
|
805 right.set(left_n, INVALID); |
|
806 } else { |
|
807 level_list[lev]=INVALID; |
|
808 } |
|
809 } |
|
810 //unlacing ends |
|
811 |
|
812 if ( level_list[lev]==INVALID ) { |
|
813 |
|
814 //gapping starts |
|
815 for (int i=lev; i!=k ; ) { |
|
816 Node v=level_list[++i]; |
|
817 while ( v!=INVALID ) { |
|
818 level.set(v,_node_num); |
|
819 v=right[v]; |
|
820 } |
|
821 level_list[i]=INVALID; |
|
822 if ( !what_heur ) first[i]=INVALID; |
|
823 } |
|
824 |
|
825 level.set(w,_node_num); |
|
826 b=lev-1; |
|
827 k=b; |
|
828 //gapping ends |
|
829 |
|
830 } else { |
|
831 |
|
832 if ( newlevel == _node_num ) level.set(w,_node_num); |
|
833 else { |
|
834 level.set(w,++newlevel); |
|
835 next.set(w,first[newlevel]); |
|
836 first[newlevel]=w; |
|
837 if ( what_heur ) b=newlevel; |
|
838 if ( k < newlevel ) ++k; //now k=newlevel |
|
839 Node z=level_list[newlevel]; |
|
840 if ( z!=INVALID ) left.set(z,w); |
|
841 right.set(w,z); |
|
842 left.set(w,INVALID); |
|
843 level_list[newlevel]=w; |
|
844 } |
|
845 } |
|
846 } //relabel |
|
847 |
|
848 }; |
|
849 |
|
850 ///Function type interface for Preflow algorithm. |
|
851 |
|
852 /// \ingroup flowalgs |
|
853 ///Function type interface for Preflow algorithm. |
|
854 ///\sa Preflow |
|
855 template<class GR, class CM, class FM> |
|
856 Preflow<GR,typename CM::Value,CM,FM> preflow(const GR &g, |
|
857 typename GR::Node source, |
|
858 typename GR::Node target, |
|
859 const CM &cap, |
|
860 FM &flow |
|
861 ) |
|
862 { |
|
863 return Preflow<GR,typename CM::Value,CM,FM>(g,source,target,cap,flow); |
|
864 } |
|
865 |
|
866 } //namespace lemon |
|
867 |
|
868 #endif //LEMON_PREFLOW_H |
|