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