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