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