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