1 // -*- C++ -*- |
|
2 |
|
3 //run gyorsan tudna adni a minmincutot a 2 fazis elejen , ne vegyuk be konstruktorba egy cutmapet? |
|
4 //constzero jo igy? |
|
5 |
|
6 //majd marci megmondja betegyem-e bfs-t meg resgraphot |
|
7 |
|
8 /* |
|
9 Heuristics: |
|
10 2 phase |
|
11 gap |
|
12 list 'level_list' on the nodes on level i implemented by hand |
|
13 stack 'active' on the active nodes on level i implemented by hand |
|
14 runs heuristic 'highest label' for H1*n relabels |
|
15 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label' |
|
16 |
|
17 Parameters H0 and H1 are initialized to 20 and 10. |
|
18 |
|
19 Constructors: |
|
20 |
|
21 Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if |
|
22 FlowMap is not constant zero, and should be true if it is |
|
23 |
|
24 Members: |
|
25 |
|
26 void run() |
|
27 |
|
28 T flowValue() : returns the value of a maximum flow |
|
29 |
|
30 void minMinCut(CutMap& M) : sets M to the characteristic vector of the |
|
31 minimum min cut. M should be a map of bools initialized to false. |
|
32 |
|
33 void maxMinCut(CutMap& M) : sets M to the characteristic vector of the |
|
34 maximum min cut. M should be a map of bools initialized to false. |
|
35 |
|
36 void minCut(CutMap& M) : sets M to the characteristic vector of |
|
37 a min cut. M should be a map of bools initialized to false. |
|
38 |
|
39 FIXME reset |
|
40 |
|
41 */ |
|
42 |
|
43 #ifndef HUGO_PREFLOW_PROBA_H |
|
44 #define HUGO_PREFLOW_PROBA_H |
|
45 |
|
46 #define H0 20 |
|
47 #define H1 1 |
|
48 |
|
49 #include <vector> |
|
50 #include <queue> |
|
51 #include<graph_wrapper.h> |
|
52 |
|
53 namespace hugo { |
|
54 |
|
55 template <typename Graph, typename T, |
|
56 typename CapMap=typename Graph::EdgeMap<T>, |
|
57 typename FlowMap=typename Graph::EdgeMap<T> > |
|
58 class PreflowProba { |
|
59 |
|
60 typedef typename Graph::Node Node; |
|
61 typedef typename Graph::Edge Edge; |
|
62 typedef typename Graph::NodeIt NodeIt; |
|
63 typedef typename Graph::OutEdgeIt OutEdgeIt; |
|
64 typedef typename Graph::InEdgeIt InEdgeIt; |
|
65 |
|
66 const Graph& G; |
|
67 Node s; |
|
68 Node t; |
|
69 const CapMap& capacity; |
|
70 FlowMap& flow; |
|
71 T value; |
|
72 bool constzero; |
|
73 bool res; |
|
74 |
|
75 typedef ResGraphWrapper<const Graph, T, CapMap, FlowMap> ResGW; |
|
76 typedef typename ResGW::OutEdgeIt ResOutEdgeIt; |
|
77 typedef typename ResGW::InEdgeIt ResInEdgeIt; |
|
78 typedef typename ResGW::Edge ResEdge; |
|
79 |
|
80 public: |
|
81 PreflowProba(Graph& _G, Node _s, Node _t, CapMap& _capacity, |
|
82 FlowMap& _flow, bool _constzero, bool _res ) : |
|
83 G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero), res(_res) {} |
|
84 |
|
85 |
|
86 void run() { |
|
87 |
|
88 ResGW res_graph(G, capacity, flow); |
|
89 |
|
90 value=0; //for the subsequent runs |
|
91 |
|
92 bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd |
|
93 int n=G.nodeNum(); |
|
94 int heur0=(int)(H0*n); //time while running 'bound decrease' |
|
95 int heur1=(int)(H1*n); //time while running 'highest label' |
|
96 int heur=heur1; //starting time interval (#of relabels) |
|
97 bool what_heur=1; |
|
98 /* |
|
99 what_heur is 0 in case 'bound decrease' |
|
100 and 1 in case 'highest label' |
|
101 */ |
|
102 bool end=false; |
|
103 /* |
|
104 Needed for 'bound decrease', 'true' |
|
105 means no active nodes are above bound b. |
|
106 */ |
|
107 int relabel=0; |
|
108 int k=n-2; //bound on the highest level under n containing a node |
|
109 int b=k; //bound on the highest level under n of an active node |
|
110 |
|
111 typename Graph::NodeMap<int> level(G,n); |
|
112 typename Graph::NodeMap<T> excess(G); |
|
113 |
|
114 std::vector<Node> active(n-1,INVALID); |
|
115 typename Graph::NodeMap<Node> next(G,INVALID); |
|
116 //Stack of the active nodes in level i < n. |
|
117 //We use it in both phases. |
|
118 |
|
119 typename Graph::NodeMap<Node> left(G,INVALID); |
|
120 typename Graph::NodeMap<Node> right(G,INVALID); |
|
121 std::vector<Node> level_list(n,INVALID); |
|
122 /* |
|
123 List of the nodes in level i<n. |
|
124 */ |
|
125 |
|
126 |
|
127 if ( constzero ) { |
|
128 |
|
129 /*Reverse_bfs from t, to find the starting level.*/ |
|
130 level.set(t,0); |
|
131 std::queue<Node> bfs_queue; |
|
132 bfs_queue.push(t); |
|
133 |
|
134 while (!bfs_queue.empty()) { |
|
135 |
|
136 Node v=bfs_queue.front(); |
|
137 bfs_queue.pop(); |
|
138 int l=level[v]+1; |
|
139 |
|
140 InEdgeIt e; |
|
141 for(G.first(e,v); G.valid(e); G.next(e)) { |
|
142 Node w=G.tail(e); |
|
143 if ( level[w] == n && w != s ) { |
|
144 bfs_queue.push(w); |
|
145 Node first=level_list[l]; |
|
146 if ( G.valid(first) ) left.set(first,w); |
|
147 right.set(w,first); |
|
148 level_list[l]=w; |
|
149 level.set(w, l); |
|
150 } |
|
151 } |
|
152 } |
|
153 |
|
154 //the starting flow |
|
155 OutEdgeIt e; |
|
156 for(G.first(e,s); G.valid(e); G.next(e)) |
|
157 { |
|
158 T c=capacity[e]; |
|
159 if ( c == 0 ) continue; |
|
160 Node w=G.head(e); |
|
161 if ( level[w] < n ) { |
|
162 if ( excess[w] == 0 && w!=t ) { |
|
163 next.set(w,active[level[w]]); |
|
164 active[level[w]]=w; |
|
165 } |
|
166 flow.set(e, c); |
|
167 excess.set(w, excess[w]+c); |
|
168 } |
|
169 } |
|
170 } |
|
171 else |
|
172 { |
|
173 |
|
174 /* |
|
175 Reverse_bfs from t in the residual graph, |
|
176 to find the starting level. |
|
177 */ |
|
178 level.set(t,0); |
|
179 std::queue<Node> bfs_queue; |
|
180 bfs_queue.push(t); |
|
181 |
|
182 while (!bfs_queue.empty()) { |
|
183 |
|
184 Node v=bfs_queue.front(); |
|
185 bfs_queue.pop(); |
|
186 int l=level[v]+1; |
|
187 |
|
188 InEdgeIt e; |
|
189 for(G.first(e,v); G.valid(e); G.next(e)) { |
|
190 if ( capacity[e] == flow[e] ) continue; |
|
191 Node w=G.tail(e); |
|
192 if ( level[w] == n && w != s ) { |
|
193 bfs_queue.push(w); |
|
194 Node first=level_list[l]; |
|
195 if ( G.valid(first) ) left.set(first,w); |
|
196 right.set(w,first); |
|
197 level_list[l]=w; |
|
198 level.set(w, l); |
|
199 } |
|
200 } |
|
201 |
|
202 OutEdgeIt f; |
|
203 for(G.first(f,v); G.valid(f); G.next(f)) { |
|
204 if ( 0 == flow[f] ) continue; |
|
205 Node w=G.head(f); |
|
206 if ( level[w] == n && w != s ) { |
|
207 bfs_queue.push(w); |
|
208 Node first=level_list[l]; |
|
209 if ( G.valid(first) ) left.set(first,w); |
|
210 right.set(w,first); |
|
211 level_list[l]=w; |
|
212 level.set(w, l); |
|
213 } |
|
214 } |
|
215 } |
|
216 |
|
217 |
|
218 /* |
|
219 Counting the excess |
|
220 */ |
|
221 NodeIt v; |
|
222 for(G.first(v); G.valid(v); G.next(v)) { |
|
223 T exc=0; |
|
224 |
|
225 InEdgeIt e; |
|
226 for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e]; |
|
227 OutEdgeIt f; |
|
228 for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[e]; |
|
229 |
|
230 excess.set(v,exc); |
|
231 |
|
232 //putting the active nodes into the stack |
|
233 int lev=level[v]; |
|
234 if ( exc > 0 && lev < n ) { |
|
235 next.set(v,active[lev]); |
|
236 active[lev]=v; |
|
237 } |
|
238 } |
|
239 |
|
240 |
|
241 //the starting flow |
|
242 OutEdgeIt e; |
|
243 for(G.first(e,s); G.valid(e); G.next(e)) |
|
244 { |
|
245 T rem=capacity[e]-flow[e]; |
|
246 if ( rem == 0 ) continue; |
|
247 Node w=G.head(e); |
|
248 if ( level[w] < n ) { |
|
249 if ( excess[w] == 0 && w!=t ) { |
|
250 next.set(w,active[level[w]]); |
|
251 active[level[w]]=w; |
|
252 } |
|
253 flow.set(e, capacity[e]); |
|
254 excess.set(w, excess[w]+rem); |
|
255 } |
|
256 } |
|
257 |
|
258 InEdgeIt f; |
|
259 for(G.first(f,s); G.valid(f); G.next(f)) |
|
260 { |
|
261 if ( flow[f] == 0 ) continue; |
|
262 Node w=G.head(f); |
|
263 if ( level[w] < n ) { |
|
264 if ( excess[w] == 0 && w!=t ) { |
|
265 next.set(w,active[level[w]]); |
|
266 active[level[w]]=w; |
|
267 } |
|
268 excess.set(w, excess[w]+flow[f]); |
|
269 flow.set(f, 0); |
|
270 } |
|
271 } |
|
272 } |
|
273 |
|
274 |
|
275 |
|
276 |
|
277 /* |
|
278 End of preprocessing |
|
279 */ |
|
280 |
|
281 |
|
282 |
|
283 /* |
|
284 Push/relabel on the highest level active nodes. |
|
285 */ |
|
286 while ( true ) { |
|
287 |
|
288 if ( b == 0 ) { |
|
289 if ( phase ) break; |
|
290 |
|
291 if ( !what_heur && !end && k > 0 ) { |
|
292 b=k; |
|
293 end=true; |
|
294 } else { |
|
295 phase=1; |
|
296 level.set(s,0); |
|
297 std::queue<Node> bfs_queue; |
|
298 bfs_queue.push(s); |
|
299 |
|
300 while (!bfs_queue.empty()) { |
|
301 |
|
302 Node v=bfs_queue.front(); |
|
303 bfs_queue.pop(); |
|
304 int l=level[v]+1; |
|
305 |
|
306 if (res){ |
|
307 ResInEdgeIt e; |
|
308 for(res_graph.first(e,v); res_graph.valid(e); |
|
309 res_graph.next(e)) { |
|
310 Node u=res_graph.tail(e); |
|
311 if ( level[u] >= n ) { |
|
312 bfs_queue.push(u); |
|
313 level.set(u, l); |
|
314 if ( excess[u] > 0 ) { |
|
315 next.set(u,active[l]); |
|
316 active[l]=u; |
|
317 } |
|
318 } |
|
319 } |
|
320 } else { |
|
321 InEdgeIt e; |
|
322 for(G.first(e,v); G.valid(e); G.next(e)) { |
|
323 if ( capacity[e] == flow[e] ) continue; |
|
324 Node u=G.tail(e); |
|
325 if ( level[u] >= n ) { |
|
326 bfs_queue.push(u); |
|
327 level.set(u, l); |
|
328 if ( excess[u] > 0 ) { |
|
329 next.set(u,active[l]); |
|
330 active[l]=u; |
|
331 } |
|
332 } |
|
333 } |
|
334 |
|
335 OutEdgeIt f; |
|
336 for(G.first(f,v); G.valid(f); G.next(f)) { |
|
337 if ( 0 == flow[f] ) continue; |
|
338 Node u=G.head(f); |
|
339 if ( level[u] >= n ) { |
|
340 bfs_queue.push(u); |
|
341 level.set(u, l); |
|
342 if ( excess[u] > 0 ) { |
|
343 next.set(u,active[l]); |
|
344 active[l]=u; |
|
345 } |
|
346 } |
|
347 } |
|
348 } |
|
349 } |
|
350 b=n-2; |
|
351 } |
|
352 |
|
353 } |
|
354 |
|
355 |
|
356 if ( !G.valid(active[b]) ) --b; |
|
357 else { |
|
358 end=false; |
|
359 |
|
360 Node w=active[b]; |
|
361 active[b]=next[w]; |
|
362 int lev=level[w]; |
|
363 T exc=excess[w]; |
|
364 int newlevel=n; //bound on the next level of w |
|
365 |
|
366 OutEdgeIt e; |
|
367 for(G.first(e,w); G.valid(e); G.next(e)) { |
|
368 |
|
369 if ( flow[e] == capacity[e] ) continue; |
|
370 Node v=G.head(e); |
|
371 //e=wv |
|
372 |
|
373 if( lev > level[v] ) { |
|
374 /*Push is allowed now*/ |
|
375 |
|
376 if ( excess[v]==0 && v!=t && v!=s ) { |
|
377 int lev_v=level[v]; |
|
378 next.set(v,active[lev_v]); |
|
379 active[lev_v]=v; |
|
380 } |
|
381 |
|
382 T cap=capacity[e]; |
|
383 T flo=flow[e]; |
|
384 T remcap=cap-flo; |
|
385 |
|
386 if ( remcap >= exc ) { |
|
387 /*A nonsaturating push.*/ |
|
388 |
|
389 flow.set(e, flo+exc); |
|
390 excess.set(v, excess[v]+exc); |
|
391 exc=0; |
|
392 break; |
|
393 |
|
394 } else { |
|
395 /*A saturating push.*/ |
|
396 |
|
397 flow.set(e, cap); |
|
398 excess.set(v, excess[v]+remcap); |
|
399 exc-=remcap; |
|
400 } |
|
401 } else if ( newlevel > level[v] ){ |
|
402 newlevel = level[v]; |
|
403 } |
|
404 |
|
405 } //for out edges wv |
|
406 |
|
407 |
|
408 if ( exc > 0 ) { |
|
409 InEdgeIt e; |
|
410 for(G.first(e,w); G.valid(e); G.next(e)) { |
|
411 |
|
412 if( flow[e] == 0 ) continue; |
|
413 Node v=G.tail(e); |
|
414 //e=vw |
|
415 |
|
416 if( lev > level[v] ) { |
|
417 /*Push is allowed now*/ |
|
418 |
|
419 if ( excess[v]==0 && v!=t && v!=s ) { |
|
420 int lev_v=level[v]; |
|
421 next.set(v,active[lev_v]); |
|
422 active[lev_v]=v; |
|
423 } |
|
424 |
|
425 T flo=flow[e]; |
|
426 |
|
427 if ( flo >= exc ) { |
|
428 /*A nonsaturating push.*/ |
|
429 |
|
430 flow.set(e, flo-exc); |
|
431 excess.set(v, excess[v]+exc); |
|
432 exc=0; |
|
433 break; |
|
434 } else { |
|
435 /*A saturating push.*/ |
|
436 |
|
437 excess.set(v, excess[v]+flo); |
|
438 exc-=flo; |
|
439 flow.set(e,0); |
|
440 } |
|
441 } else if ( newlevel > level[v] ) { |
|
442 newlevel = level[v]; |
|
443 } |
|
444 } //for in edges vw |
|
445 |
|
446 } // if w still has excess after the out edge for cycle |
|
447 |
|
448 excess.set(w, exc); |
|
449 |
|
450 /* |
|
451 Relabel |
|
452 */ |
|
453 |
|
454 |
|
455 if ( exc > 0 ) { |
|
456 //now 'lev' is the old level of w |
|
457 |
|
458 if ( phase ) { |
|
459 level.set(w,++newlevel); |
|
460 next.set(w,active[newlevel]); |
|
461 active[newlevel]=w; |
|
462 b=newlevel; |
|
463 } else { |
|
464 //unlacing starts |
|
465 Node right_n=right[w]; |
|
466 Node left_n=left[w]; |
|
467 |
|
468 if ( G.valid(right_n) ) { |
|
469 if ( G.valid(left_n) ) { |
|
470 right.set(left_n, right_n); |
|
471 left.set(right_n, left_n); |
|
472 } else { |
|
473 level_list[lev]=right_n; |
|
474 left.set(right_n, INVALID); |
|
475 } |
|
476 } else { |
|
477 if ( G.valid(left_n) ) { |
|
478 right.set(left_n, INVALID); |
|
479 } else { |
|
480 level_list[lev]=INVALID; |
|
481 } |
|
482 } |
|
483 //unlacing ends |
|
484 |
|
485 if ( !G.valid(level_list[lev]) ) { |
|
486 |
|
487 //gapping starts |
|
488 for (int i=lev; i!=k ; ) { |
|
489 Node v=level_list[++i]; |
|
490 while ( G.valid(v) ) { |
|
491 level.set(v,n); |
|
492 v=right[v]; |
|
493 } |
|
494 level_list[i]=INVALID; |
|
495 if ( !what_heur ) active[i]=INVALID; |
|
496 } |
|
497 |
|
498 level.set(w,n); |
|
499 b=lev-1; |
|
500 k=b; |
|
501 //gapping ends |
|
502 |
|
503 } else { |
|
504 |
|
505 if ( newlevel == n ) level.set(w,n); |
|
506 else { |
|
507 level.set(w,++newlevel); |
|
508 next.set(w,active[newlevel]); |
|
509 active[newlevel]=w; |
|
510 if ( what_heur ) b=newlevel; |
|
511 if ( k < newlevel ) ++k; //now k=newlevel |
|
512 Node first=level_list[newlevel]; |
|
513 if ( G.valid(first) ) left.set(first,w); |
|
514 right.set(w,first); |
|
515 left.set(w,INVALID); |
|
516 level_list[newlevel]=w; |
|
517 } |
|
518 } |
|
519 |
|
520 |
|
521 ++relabel; |
|
522 if ( relabel >= heur ) { |
|
523 relabel=0; |
|
524 if ( what_heur ) { |
|
525 what_heur=0; |
|
526 heur=heur0; |
|
527 end=false; |
|
528 } else { |
|
529 what_heur=1; |
|
530 heur=heur1; |
|
531 b=k; |
|
532 } |
|
533 } |
|
534 } //phase 0 |
|
535 |
|
536 |
|
537 } // if ( exc > 0 ) |
|
538 |
|
539 |
|
540 } // if stack[b] is nonempty |
|
541 |
|
542 } // while(true) |
|
543 |
|
544 |
|
545 value = excess[t]; |
|
546 /*Max flow value.*/ |
|
547 |
|
548 } //void run() |
|
549 |
|
550 |
|
551 |
|
552 |
|
553 |
|
554 /* |
|
555 Returns the maximum value of a flow. |
|
556 */ |
|
557 |
|
558 T flowValue() { |
|
559 return value; |
|
560 } |
|
561 |
|
562 |
|
563 FlowMap Flow() { |
|
564 return flow; |
|
565 } |
|
566 |
|
567 |
|
568 |
|
569 void Flow(FlowMap& _flow ) { |
|
570 NodeIt v; |
|
571 for(G.first(v) ; G.valid(v); G.next(v)) |
|
572 _flow.set(v,flow[v]); |
|
573 } |
|
574 |
|
575 |
|
576 |
|
577 /* |
|
578 Returns the minimum min cut, by a bfs from s in the residual graph. |
|
579 */ |
|
580 |
|
581 template<typename _CutMap> |
|
582 void minMinCut(_CutMap& M) { |
|
583 |
|
584 std::queue<Node> queue; |
|
585 |
|
586 M.set(s,true); |
|
587 queue.push(s); |
|
588 |
|
589 while (!queue.empty()) { |
|
590 Node w=queue.front(); |
|
591 queue.pop(); |
|
592 |
|
593 OutEdgeIt e; |
|
594 for(G.first(e,w) ; G.valid(e); G.next(e)) { |
|
595 Node v=G.head(e); |
|
596 if (!M[v] && flow[e] < capacity[e] ) { |
|
597 queue.push(v); |
|
598 M.set(v, true); |
|
599 } |
|
600 } |
|
601 |
|
602 InEdgeIt f; |
|
603 for(G.first(f,w) ; G.valid(f); G.next(f)) { |
|
604 Node v=G.tail(f); |
|
605 if (!M[v] && flow[f] > 0 ) { |
|
606 queue.push(v); |
|
607 M.set(v, true); |
|
608 } |
|
609 } |
|
610 } |
|
611 } |
|
612 |
|
613 |
|
614 |
|
615 /* |
|
616 Returns the maximum min cut, by a reverse bfs |
|
617 from t in the residual graph. |
|
618 */ |
|
619 |
|
620 template<typename _CutMap> |
|
621 void maxMinCut(_CutMap& M) { |
|
622 |
|
623 std::queue<Node> queue; |
|
624 |
|
625 M.set(t,true); |
|
626 queue.push(t); |
|
627 |
|
628 while (!queue.empty()) { |
|
629 Node w=queue.front(); |
|
630 queue.pop(); |
|
631 |
|
632 |
|
633 InEdgeIt e; |
|
634 for(G.first(e,w) ; G.valid(e); G.next(e)) { |
|
635 Node v=G.tail(e); |
|
636 if (!M[v] && flow[e] < capacity[e] ) { |
|
637 queue.push(v); |
|
638 M.set(v, true); |
|
639 } |
|
640 } |
|
641 |
|
642 OutEdgeIt f; |
|
643 for(G.first(f,w) ; G.valid(f); G.next(f)) { |
|
644 Node v=G.head(f); |
|
645 if (!M[v] && flow[f] > 0 ) { |
|
646 queue.push(v); |
|
647 M.set(v, true); |
|
648 } |
|
649 } |
|
650 } |
|
651 |
|
652 NodeIt v; |
|
653 for(G.first(v) ; G.valid(v); G.next(v)) { |
|
654 M.set(v, !M[v]); |
|
655 } |
|
656 |
|
657 } |
|
658 |
|
659 |
|
660 |
|
661 template<typename CutMap> |
|
662 void minCut(CutMap& M) { |
|
663 minMinCut(M); |
|
664 } |
|
665 |
|
666 |
|
667 void reset_target (Node _t) {t=_t;} |
|
668 void reset_source (Node _s) {s=_s;} |
|
669 |
|
670 template<typename _CapMap> |
|
671 void reset_cap (_CapMap _cap) {capacity=_cap;} |
|
672 |
|
673 template<typename _FlowMap> |
|
674 void reset_cap (_FlowMap _flow, bool _constzero) { |
|
675 flow=_flow; |
|
676 constzero=_constzero; |
|
677 } |
|
678 |
|
679 |
|
680 |
|
681 }; |
|
682 |
|
683 } //namespace hugo |
|
684 |
|
685 #endif //PREFLOW_PROBA_H |
|
686 |
|
687 |
|
688 |
|
689 |
|