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