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