|
1 // -*- C++ -*- |
|
2 |
|
3 //hogy kell azt megcsinalni hogy a konstruktorban a FlowMap-nek legyen default erteke (a 0 map) |
|
4 |
|
5 /* |
|
6 preflow_aug.h |
|
7 by jacint. |
|
8 |
|
9 The same as preflow.h, the only difference is that here |
|
10 we start from a given flow. |
|
11 |
|
12 |
|
13 The constructor runs the algorithm: |
|
14 |
|
15 template <typename Graph, typename T, |
|
16 typename Graph_EdgeMap_T_1=typename Graph::EdgeMap<T>, |
|
17 typename Graph_EdgeMap_T_2=typename Graph::EdgeMap<T> > |
|
18 PreflowAug(Graph G, G_NodeIt s, G_NodeIt t, G_EdgeMap_T_1 capacity, |
|
19 G_EdgeMap_T_2 flow, bool flow_type) |
|
20 |
|
21 'capacity' must be non-negative, if flow_type=0 then |
|
22 'flow' must be a flow, otherwise it must be a preflow. |
|
23 |
|
24 |
|
25 Members: |
|
26 |
|
27 T maxFlow() : returns the value of a maximum flow |
|
28 |
|
29 T flow(G_EdgeIt e) : for a fixed maximum flow x it returns x(e) |
|
30 |
|
31 G_EdgeMap_T_2 flow() : returns the fixed maximum flow x |
|
32 |
|
33 void minMinCut(G_NodeMap_bool& M) : |
|
34 sets M to the characteristic vector of the |
|
35 minimum min cut. M must be initialized to false. |
|
36 |
|
37 void maxMinCut(G_NodeMap_bool& M) : |
|
38 sets M to the characteristic vector of the |
|
39 maximum min cut. M must be initialized to false. |
|
40 |
|
41 void minCut(G_NodeMap_bool& M) : |
|
42 sets M to the characteristic vector of a |
|
43 min cut. M must be initialized to false. |
|
44 */ |
|
45 |
|
46 #ifndef PREFLOW_H |
|
47 #define PREFLOW_H |
|
48 |
|
49 #define H0 20 |
|
50 #define H1 1 |
|
51 |
|
52 #include <vector> |
|
53 #include <queue> |
|
54 |
|
55 #include <iostream> //for error handling |
|
56 |
|
57 #include <time_measure.h> |
|
58 |
|
59 namespace hugo { |
|
60 |
|
61 template <typename Graph, typename T, |
|
62 typename CapMap=typename Graph::EdgeMap<T>, |
|
63 typename FlowMap=typename Graph::EdgeMap<T> > |
|
64 class PreflowAug { |
|
65 |
|
66 typedef typename Graph::NodeIt NodeIt; |
|
67 typedef typename Graph::EdgeIt EdgeIt; |
|
68 typedef typename Graph::EachNodeIt EachNodeIt; |
|
69 typedef typename Graph::EachEdgeIt EachEdgeIt; |
|
70 typedef typename Graph::OutEdgeIt OutEdgeIt; |
|
71 typedef typename Graph::InEdgeIt InEdgeIt; |
|
72 |
|
73 Graph& G; |
|
74 NodeIt s; |
|
75 NodeIt t; |
|
76 CapMap& capacity; |
|
77 FlowMap& _flow; |
|
78 bool flow_type; |
|
79 T value; |
|
80 |
|
81 public: |
|
82 double time; |
|
83 PreflowAug(Graph& _G, NodeIt _s, NodeIt _t, |
|
84 CapMap& _capacity, FlowMap& __flow, bool _flow_type ) : |
|
85 G(_G), s(_s), t(_t), capacity(_capacity), _flow(__flow), |
|
86 flow_type(_flow_type) |
|
87 { |
|
88 |
|
89 |
|
90 bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd |
|
91 int n=G.nodeNum(); |
|
92 int heur0=(int)(H0*n); //time while running 'bound decrease' |
|
93 int heur1=(int)(H1*n); //time while running 'highest label' |
|
94 int heur=heur1; //starting time interval (#of relabels) |
|
95 bool what_heur=1; |
|
96 /* |
|
97 what_heur is 0 in case 'bound decrease' |
|
98 and 1 in case 'highest label' |
|
99 */ |
|
100 bool end=false; |
|
101 /* |
|
102 Needed for 'bound decrease', 'true' |
|
103 means no active nodes are above bound b. |
|
104 */ |
|
105 int relabel=0; |
|
106 int k=n-2; //bound on the highest level under n containing a node |
|
107 int b=k; //bound on the highest level under n of an active node |
|
108 |
|
109 typename Graph::NodeMap<int> level(G,n); |
|
110 typename Graph::NodeMap<T> excess(G); |
|
111 |
|
112 std::vector<NodeIt> active(n); |
|
113 typename Graph::NodeMap<NodeIt> next(G); |
|
114 //Stack of the active nodes in level i < n. |
|
115 //We use it in both phases. |
|
116 |
|
117 typename Graph::NodeMap<NodeIt> left(G); |
|
118 typename Graph::NodeMap<NodeIt> right(G); |
|
119 std::vector<NodeIt> level_list(n); |
|
120 /* |
|
121 List of the nodes in level i<n. |
|
122 */ |
|
123 |
|
124 /*Reverse_bfs from t, to find the starting level.*/ |
|
125 level.set(t,0); |
|
126 std::queue<NodeIt> bfs_queue; |
|
127 bfs_queue.push(t); |
|
128 |
|
129 |
|
130 while (!bfs_queue.empty()) { |
|
131 |
|
132 NodeIt v=bfs_queue.front(); |
|
133 bfs_queue.pop(); |
|
134 int l=level.get(v)+1; |
|
135 |
|
136 for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) { |
|
137 if ( capacity.get(e) == _flow.get(e) ) continue; |
|
138 NodeIt w=G.tail(e); |
|
139 if ( level.get(w) == n && w != s ) { |
|
140 bfs_queue.push(w); |
|
141 NodeIt first=level_list[l]; |
|
142 if ( first.valid() ) left.set(first,w); |
|
143 right.set(w,first); |
|
144 level_list[l]=w; |
|
145 level.set(w, l); |
|
146 } |
|
147 } |
|
148 |
|
149 for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) { |
|
150 if ( 0 == _flow.get(e) ) continue; |
|
151 NodeIt w=G.head(e); |
|
152 if ( level.get(w) == n && w != s ) { |
|
153 bfs_queue.push(w); |
|
154 NodeIt first=level_list[l]; |
|
155 if ( first != 0 ) left.set(first,w); |
|
156 right.set(w,first); |
|
157 level_list[l]=w; |
|
158 level.set(w, l); |
|
159 } |
|
160 } |
|
161 } |
|
162 |
|
163 level.set(s,n); |
|
164 |
|
165 |
|
166 |
|
167 /*The starting excess.*/ |
|
168 if ( flow_type ) { |
|
169 for(EachEdgeIt e=G.template first<EachEdgeIt>(); e.valid(); ++e) { |
|
170 if ( _flow.get(e) > 0 ) { |
|
171 T flo=_flow.get(e); |
|
172 NodeIt u=G.tail(e); |
|
173 NodeIt v=G.head(e); |
|
174 excess.set(u, excess.get(u)-flo); |
|
175 excess.set(v, excess.get(v)+flo); |
|
176 } |
|
177 } |
|
178 |
|
179 for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) { |
|
180 if ( excess.get(v) < 0 ) { |
|
181 std::cerr<<"It is not a pre_flow."<<std::endl; |
|
182 exit(1); |
|
183 } else { |
|
184 next.set(v,active[level.get(v)]); |
|
185 active[level.get(v)]=v;} |
|
186 } |
|
187 } |
|
188 |
|
189 |
|
190 /* Starting flow. It is everywhere 0 at the moment. */ |
|
191 for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) |
|
192 { |
|
193 T c=capacity.get(e)-_flow.get(e); |
|
194 if ( c == 0 ) continue; |
|
195 NodeIt w=G.head(e); |
|
196 if ( level.get(w) < n ) { |
|
197 if ( excess.get(w) == 0 && w!=t ) { |
|
198 next.set(w,active[level.get(w)]); |
|
199 active[level.get(w)]=w; |
|
200 } |
|
201 _flow.set(e, capacity.get(e)); |
|
202 excess.set(w, excess.get(w)+c); |
|
203 } |
|
204 } |
|
205 |
|
206 /* |
|
207 End of preprocessing |
|
208 */ |
|
209 |
|
210 |
|
211 |
|
212 /* |
|
213 Push/relabel on the highest level active nodes. |
|
214 */ |
|
215 while ( true ) { |
|
216 |
|
217 if ( b == 0 ) { |
|
218 if ( phase ) break; |
|
219 |
|
220 if ( !what_heur && !end && k > 0 ) { |
|
221 b=k; |
|
222 end=true; |
|
223 } else { |
|
224 phase=1; |
|
225 time=currTime(); |
|
226 level.set(s,0); |
|
227 std::queue<NodeIt> bfs_queue; |
|
228 bfs_queue.push(s); |
|
229 |
|
230 while (!bfs_queue.empty()) { |
|
231 |
|
232 NodeIt v=bfs_queue.front(); |
|
233 bfs_queue.pop(); |
|
234 int l=level.get(v)+1; |
|
235 |
|
236 for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) { |
|
237 if ( capacity.get(e) == _flow.get(e) ) continue; |
|
238 NodeIt u=G.tail(e); |
|
239 if ( level.get(u) >= n ) { |
|
240 bfs_queue.push(u); |
|
241 level.set(u, l); |
|
242 if ( excess.get(u) > 0 ) { |
|
243 next.set(u,active[l]); |
|
244 active[l]=u; |
|
245 } |
|
246 } |
|
247 } |
|
248 |
|
249 for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) { |
|
250 if ( 0 == _flow.get(e) ) continue; |
|
251 NodeIt u=G.head(e); |
|
252 if ( level.get(u) >= n ) { |
|
253 bfs_queue.push(u); |
|
254 level.set(u, l); |
|
255 if ( excess.get(u) > 0 ) { |
|
256 next.set(u,active[l]); |
|
257 active[l]=u; |
|
258 } |
|
259 } |
|
260 } |
|
261 } |
|
262 b=n-2; |
|
263 } |
|
264 |
|
265 } |
|
266 |
|
267 |
|
268 if ( active[b] == 0 ) --b; |
|
269 else { |
|
270 end=false; |
|
271 |
|
272 NodeIt w=active[b]; |
|
273 active[b]=next.get(w); |
|
274 int lev=level.get(w); |
|
275 T exc=excess.get(w); |
|
276 int newlevel=n; //bound on the next level of w |
|
277 |
|
278 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) { |
|
279 |
|
280 if ( _flow.get(e) == capacity.get(e) ) continue; |
|
281 NodeIt v=G.head(e); |
|
282 //e=wv |
|
283 |
|
284 if( lev > level.get(v) ) { |
|
285 /*Push is allowed now*/ |
|
286 |
|
287 if ( excess.get(v)==0 && v!=t && v!=s ) { |
|
288 int lev_v=level.get(v); |
|
289 next.set(v,active[lev_v]); |
|
290 active[lev_v]=v; |
|
291 } |
|
292 |
|
293 T cap=capacity.get(e); |
|
294 T flo=_flow.get(e); |
|
295 T remcap=cap-flo; |
|
296 |
|
297 if ( remcap >= exc ) { |
|
298 /*A nonsaturating push.*/ |
|
299 |
|
300 _flow.set(e, flo+exc); |
|
301 excess.set(v, excess.get(v)+exc); |
|
302 exc=0; |
|
303 break; |
|
304 |
|
305 } else { |
|
306 /*A saturating push.*/ |
|
307 |
|
308 _flow.set(e, cap); |
|
309 excess.set(v, excess.get(v)+remcap); |
|
310 exc-=remcap; |
|
311 } |
|
312 } else if ( newlevel > level.get(v) ){ |
|
313 newlevel = level.get(v); |
|
314 } |
|
315 |
|
316 } //for out edges wv |
|
317 |
|
318 |
|
319 if ( exc > 0 ) { |
|
320 for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) { |
|
321 |
|
322 if( _flow.get(e) == 0 ) continue; |
|
323 NodeIt v=G.tail(e); |
|
324 //e=vw |
|
325 |
|
326 if( lev > level.get(v) ) { |
|
327 /*Push is allowed now*/ |
|
328 |
|
329 if ( excess.get(v)==0 && v!=t && v!=s ) { |
|
330 int lev_v=level.get(v); |
|
331 next.set(v,active[lev_v]); |
|
332 active[lev_v]=v; |
|
333 } |
|
334 |
|
335 T flo=_flow.get(e); |
|
336 |
|
337 if ( flo >= exc ) { |
|
338 /*A nonsaturating push.*/ |
|
339 |
|
340 _flow.set(e, flo-exc); |
|
341 excess.set(v, excess.get(v)+exc); |
|
342 exc=0; |
|
343 break; |
|
344 } else { |
|
345 /*A saturating push.*/ |
|
346 |
|
347 excess.set(v, excess.get(v)+flo); |
|
348 exc-=flo; |
|
349 _flow.set(e,0); |
|
350 } |
|
351 } else if ( newlevel > level.get(v) ) { |
|
352 newlevel = level.get(v); |
|
353 } |
|
354 } //for in edges vw |
|
355 |
|
356 } // if w still has excess after the out edge for cycle |
|
357 |
|
358 excess.set(w, exc); |
|
359 |
|
360 /* |
|
361 Relabel |
|
362 */ |
|
363 |
|
364 |
|
365 if ( exc > 0 ) { |
|
366 //now 'lev' is the old level of w |
|
367 |
|
368 if ( phase ) { |
|
369 level.set(w,++newlevel); |
|
370 next.set(w,active[newlevel]); |
|
371 active[newlevel]=w; |
|
372 b=newlevel; |
|
373 } else { |
|
374 //unlacing starts |
|
375 NodeIt right_n=right.get(w); |
|
376 NodeIt left_n=left.get(w); |
|
377 |
|
378 if ( right_n != 0 ) { |
|
379 if ( left_n != 0 ) { |
|
380 right.set(left_n, right_n); |
|
381 left.set(right_n, left_n); |
|
382 } else { |
|
383 level_list[lev]=right_n; |
|
384 left.set(right_n, 0); |
|
385 } |
|
386 } else { |
|
387 if ( left_n != 0 ) { |
|
388 right.set(left_n, 0); |
|
389 } else { |
|
390 level_list[lev]=0; |
|
391 |
|
392 } |
|
393 } |
|
394 //unlacing ends |
|
395 |
|
396 //gapping starts |
|
397 if ( level_list[lev]==0 ) { |
|
398 |
|
399 for (int i=lev; i!=k ; ) { |
|
400 NodeIt v=level_list[++i]; |
|
401 while ( v != 0 ) { |
|
402 level.set(v,n); |
|
403 v=right.get(v); |
|
404 } |
|
405 level_list[i]=0; |
|
406 if ( !what_heur ) active[i]=0; |
|
407 } |
|
408 |
|
409 level.set(w,n); |
|
410 b=lev-1; |
|
411 k=b; |
|
412 //gapping ends |
|
413 } else { |
|
414 |
|
415 if ( newlevel == n ) level.set(w,n); |
|
416 else { |
|
417 level.set(w,++newlevel); |
|
418 next.set(w,active[newlevel]); |
|
419 active[newlevel]=w; |
|
420 if ( what_heur ) b=newlevel; |
|
421 if ( k < newlevel ) ++k; |
|
422 NodeIt first=level_list[newlevel]; |
|
423 if ( first != 0 ) left.set(first,w); |
|
424 right.set(w,first); |
|
425 left.set(w,0); |
|
426 level_list[newlevel]=w; |
|
427 } |
|
428 } |
|
429 |
|
430 |
|
431 ++relabel; |
|
432 if ( relabel >= heur ) { |
|
433 relabel=0; |
|
434 if ( what_heur ) { |
|
435 what_heur=0; |
|
436 heur=heur0; |
|
437 end=false; |
|
438 } else { |
|
439 what_heur=1; |
|
440 heur=heur1; |
|
441 b=k; |
|
442 } |
|
443 } |
|
444 } //phase 0 |
|
445 |
|
446 |
|
447 } // if ( exc > 0 ) |
|
448 |
|
449 |
|
450 } // if stack[b] is nonempty |
|
451 |
|
452 } // while(true) |
|
453 |
|
454 |
|
455 value = excess.get(t); |
|
456 /*Max flow value.*/ |
|
457 |
|
458 } //void run() |
|
459 |
|
460 |
|
461 |
|
462 |
|
463 |
|
464 /* |
|
465 Returns the maximum value of a flow. |
|
466 */ |
|
467 |
|
468 T maxFlow() { |
|
469 return value; |
|
470 } |
|
471 |
|
472 |
|
473 |
|
474 /* |
|
475 For the maximum flow x found by the algorithm, |
|
476 it returns the flow value on edge e, i.e. x(e). |
|
477 */ |
|
478 |
|
479 T flow(EdgeIt e) { |
|
480 return _flow.get(e); |
|
481 } |
|
482 |
|
483 |
|
484 |
|
485 FlowMap flow() { |
|
486 return _flow; |
|
487 } |
|
488 |
|
489 |
|
490 |
|
491 void flow(FlowMap& __flow ) { |
|
492 for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) |
|
493 __flow.set(v,_flow.get(v)); |
|
494 } |
|
495 |
|
496 |
|
497 |
|
498 /* |
|
499 Returns the minimum min cut, by a bfs from s in the residual graph. |
|
500 */ |
|
501 |
|
502 template<typename _CutMap> |
|
503 void minMinCut(_CutMap& M) { |
|
504 |
|
505 std::queue<NodeIt> queue; |
|
506 |
|
507 M.set(s,true); |
|
508 queue.push(s); |
|
509 |
|
510 while (!queue.empty()) { |
|
511 NodeIt w=queue.front(); |
|
512 queue.pop(); |
|
513 |
|
514 for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) { |
|
515 NodeIt v=G.head(e); |
|
516 if (!M.get(v) && _flow.get(e) < capacity.get(e) ) { |
|
517 queue.push(v); |
|
518 M.set(v, true); |
|
519 } |
|
520 } |
|
521 |
|
522 for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) { |
|
523 NodeIt v=G.tail(e); |
|
524 if (!M.get(v) && _flow.get(e) > 0 ) { |
|
525 queue.push(v); |
|
526 M.set(v, true); |
|
527 } |
|
528 } |
|
529 } |
|
530 } |
|
531 |
|
532 |
|
533 |
|
534 /* |
|
535 Returns the maximum min cut, by a reverse bfs |
|
536 from t in the residual graph. |
|
537 */ |
|
538 |
|
539 template<typename _CutMap> |
|
540 void maxMinCut(_CutMap& M) { |
|
541 |
|
542 std::queue<NodeIt> queue; |
|
543 |
|
544 M.set(t,true); |
|
545 queue.push(t); |
|
546 |
|
547 while (!queue.empty()) { |
|
548 NodeIt w=queue.front(); |
|
549 queue.pop(); |
|
550 |
|
551 for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) { |
|
552 NodeIt v=G.tail(e); |
|
553 if (!M.get(v) && _flow.get(e) < capacity.get(e) ) { |
|
554 queue.push(v); |
|
555 M.set(v, true); |
|
556 } |
|
557 } |
|
558 |
|
559 for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) { |
|
560 NodeIt v=G.head(e); |
|
561 if (!M.get(v) && _flow.get(e) > 0 ) { |
|
562 queue.push(v); |
|
563 M.set(v, true); |
|
564 } |
|
565 } |
|
566 } |
|
567 |
|
568 for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) { |
|
569 M.set(v, !M.get(v)); |
|
570 } |
|
571 |
|
572 } |
|
573 |
|
574 |
|
575 |
|
576 template<typename CutMap> |
|
577 void minCut(CutMap& M) { |
|
578 minMinCut(M); |
|
579 } |
|
580 |
|
581 |
|
582 }; |
|
583 }//namespace marci |
|
584 #endif |
|
585 |
|
586 |
|
587 |
|
588 |