13 |
15 |
14 T maxflow() : returns the value of a maximum flow |
16 T maxflow() : returns the value of a maximum flow |
15 |
17 |
16 T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e) |
18 T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e) |
17 |
19 |
18 Graph::EdgeMap<T> allflow() : returns the fixed maximum flow x |
20 FlowMap allflow() : returns the fixed maximum flow x |
19 |
21 |
20 Graph::NodeMap<bool> mincut() : returns a |
22 void mincut(CutMap& M) : sets M to the characteristic vector of a |
21 characteristic vector of a minimum cut. (An empty level |
23 minimum cut. M should be a map of bools initialized to false. |
22 in the algorithm gives a minimum cut.) |
24 |
|
25 void min_mincut(CutMap& M) : sets M to the characteristic vector of the |
|
26 minimum min cut. M should be a map of bools initialized to false. |
|
27 |
|
28 void max_mincut(CutMap& M) : sets M to the characteristic vector of the |
|
29 maximum min cut. M should be a map of bools initialized to false. |
|
30 |
23 */ |
31 */ |
24 |
32 |
25 #ifndef PREFLOW_PUSH_HL_H |
33 #ifndef PREFLOW_PUSH_HL_H |
26 #define PREFLOW_PUSH_HL_H |
34 #define PREFLOW_PUSH_HL_H |
27 |
35 |
28 #define A 1 |
36 #define A 1 |
29 |
37 |
30 #include <vector> |
38 #include <vector> |
31 #include <stack> |
39 #include <stack> |
32 |
40 #include <queue> |
33 #include <reverse_bfs.h> |
|
34 |
41 |
35 namespace marci { |
42 namespace marci { |
36 |
43 |
37 template <typename Graph, typename T> |
44 template <typename Graph, typename T, |
|
45 typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>, |
|
46 typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> > |
38 class preflow_push_hl { |
47 class preflow_push_hl { |
39 |
48 |
40 typedef typename Graph::NodeIt NodeIt; |
49 typedef typename Graph::NodeIt NodeIt; |
41 typedef typename Graph::EdgeIt EdgeIt; |
50 typedef typename Graph::EdgeIt EdgeIt; |
42 typedef typename Graph::EachNodeIt EachNodeIt; |
51 typedef typename Graph::EachNodeIt EachNodeIt; |
44 typedef typename Graph::InEdgeIt InEdgeIt; |
53 typedef typename Graph::InEdgeIt InEdgeIt; |
45 |
54 |
46 Graph& G; |
55 Graph& G; |
47 NodeIt s; |
56 NodeIt s; |
48 NodeIt t; |
57 NodeIt t; |
49 typename Graph::EdgeMap<T> flow; |
58 FlowMap flow; |
50 typename Graph::EdgeMap<T> capacity; |
59 CapMap& capacity; |
51 T value; |
60 T value; |
52 typename Graph::NodeMap<bool> mincutvector; |
61 |
53 |
|
54 public: |
62 public: |
55 |
63 |
56 preflow_push_hl(Graph& _G, NodeIt _s, NodeIt _t, |
64 preflow_push_hl(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) : |
57 typename Graph::EdgeMap<T>& _capacity) : |
65 G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { } |
58 G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), |
66 |
59 mincutvector(_G, true) { } |
67 |
60 |
68 |
61 |
69 |
62 /* |
|
63 The run() function runs the highest label preflow-push, |
|
64 running time: O(n^2\sqrt(m)) |
|
65 */ |
|
66 void run() { |
70 void run() { |
67 |
71 |
68 std::cout<<"A is "<<A<<" "; |
|
69 |
|
70 typename Graph::NodeMap<int> level(G); |
|
71 typename Graph::NodeMap<T> excess(G); |
|
72 |
|
73 int n=G.nodeNum(); |
72 int n=G.nodeNum(); |
74 int b=n-2; |
73 int b=n-2; |
75 /* |
74 /* |
76 b is a bound on the highest level of an active node. |
75 b is a bound on the highest level of an active node. |
77 In the beginning it is at most n-2. |
|
78 */ |
76 */ |
79 |
77 |
80 std::vector<int> numb(n); //The number of nodes on level i < n. |
78 IntMap level(G,n); |
|
79 TMap excess(G); |
|
80 |
|
81 std::vector<int> numb(n); |
|
82 /* |
|
83 The number of nodes on level i < n. It is |
|
84 initialized to n+1, because of the reverse_bfs-part. |
|
85 */ |
|
86 |
81 std::vector<std::stack<NodeIt> > stack(2*n-1); |
87 std::vector<std::stack<NodeIt> > stack(2*n-1); |
82 //Stack of the active nodes in level i. |
88 //Stack of the active nodes in level i. |
83 |
89 |
84 |
90 |
85 /*Reverse_bfs from t, to find the starting level.*/ |
91 /*Reverse_bfs from t, to find the starting level.*/ |
86 reverse_bfs<Graph> bfs(G, t); |
92 level.set(t,0); |
87 bfs.run(); |
93 std::queue<NodeIt> bfs_queue; |
88 for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) |
94 bfs_queue.push(t); |
89 { |
95 |
90 int dist=bfs.dist(v); |
96 while (!bfs_queue.empty()) { |
91 level.set(v, dist); |
97 |
92 ++numb[dist]; |
98 NodeIt v=bfs_queue.front(); |
93 } |
99 bfs_queue.pop(); |
94 |
100 int l=level.get(v)+1; |
|
101 |
|
102 for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) { |
|
103 NodeIt w=G.tail(e); |
|
104 if ( level.get(w) == n ) { |
|
105 bfs_queue.push(w); |
|
106 ++numb[l]; |
|
107 level.set(w, l); |
|
108 } |
|
109 } |
|
110 } |
|
111 |
95 level.set(s,n); |
112 level.set(s,n); |
|
113 |
96 |
114 |
97 |
115 |
98 /* Starting flow. It is everywhere 0 at the moment. */ |
116 /* Starting flow. It is everywhere 0 at the moment. */ |
99 for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) |
117 for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) |
100 { |
118 { |
101 if ( capacity.get(e) > 0 ) { |
119 if ( capacity.get(e) == 0 ) continue; |
102 NodeIt w=G.head(e); |
120 NodeIt w=G.head(e); |
103 if ( w!=s ) { |
121 if ( w!=s ) { |
104 if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); |
122 if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); |
105 flow.set(e, capacity.get(e)); |
123 flow.set(e, capacity.get(e)); |
106 excess.set(w, excess.get(w)+capacity.get(e)); |
124 excess.set(w, excess.get(w)+capacity.get(e)); |
107 } |
125 } |
108 } |
126 } |
109 } |
127 |
110 |
|
111 /* |
128 /* |
112 End of preprocessing |
129 End of preprocessing |
113 */ |
130 */ |
114 |
131 |
115 |
132 |
116 |
|
117 /* |
133 /* |
118 Push/relabel on the highest level active nodes. |
134 Push/relabel on the highest level active nodes. |
119 */ |
135 */ |
120 |
|
121 /*While there exists an active node.*/ |
136 /*While there exists an active node.*/ |
122 while (b) { |
137 while (b) { |
123 |
138 if ( stack[b].empty() ) { |
124 /*We decrease the bound if there is no active node of level b.*/ |
|
125 if (stack[b].empty()) { |
|
126 --b; |
139 --b; |
127 } else { |
140 continue; |
128 |
141 } |
129 NodeIt w=stack[b].top(); //w is a highest label active node. |
142 |
130 stack[b].pop(); |
143 NodeIt w=stack[b].top(); //w is a highest label active node. |
131 |
144 stack[b].pop(); |
132 int newlevel=2*n-2; //In newlevel we bound the next level of w. |
145 int lev=level.get(w); |
133 |
146 int exc=excess.get(w); |
|
147 int newlevel=2*n-2; //In newlevel we bound the next level of w. |
|
148 |
|
149 // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e |
134 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) { |
150 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) { |
135 |
151 |
136 if ( flow.get(e) < capacity.get(e) ) { |
152 if ( flow.get(e) == capacity.get(e) ) continue; |
137 /*e is an edge of the residual graph */ |
153 NodeIt v=G.head(e); |
138 |
154 //e=wv |
139 NodeIt v=G.head(e); /*e is the edge wv.*/ |
155 |
140 |
156 if( lev > level.get(v) ) { |
141 if( level.get(w) == level.get(v)+1 ) { |
157 /*Push is allowed now*/ |
142 /*Push is allowed now*/ |
158 |
143 |
159 if ( excess.get(v)==0 && v != s && v !=t ) |
144 if ( excess.get(v)==0 && v != s && v !=t ) stack[level.get(v)].push(v); |
160 stack[level.get(v)].push(v); |
145 /*v becomes active.*/ |
161 /*v becomes active.*/ |
146 |
162 |
147 if ( capacity.get(e)-flow.get(e) > excess.get(w) ) { |
163 int cap=capacity.get(e); |
148 /*A nonsaturating push.*/ |
164 int flo=flow.get(e); |
149 |
165 int remcap=cap-flo; |
150 flow.set(e, flow.get(e)+excess.get(w)); |
166 |
151 excess.set(v, excess.get(v)+excess.get(w)); |
167 if ( remcap >= exc ) { |
152 excess.set(w,0); |
168 /*A nonsaturating push.*/ |
153 break; |
169 |
154 |
170 flow.set(e, flo+exc); |
155 } else { |
171 excess.set(v, excess.get(v)+exc); |
156 /*A saturating push.*/ |
172 exc=0; |
157 |
173 break; |
158 excess.set(v, excess.get(v)+capacity.get(e)-flow.get(e)); |
174 |
159 excess.set(w, excess.get(w)-capacity.get(e)+flow.get(e)); |
175 } else { |
160 flow.set(e, capacity.get(e)); |
176 /*A saturating push.*/ |
161 if ( excess.get(w)==0 ) break; |
177 |
162 /*If w is not active any more, then we go on to the next node.*/ |
178 flow.set(e, cap ); |
163 |
179 excess.set(v, excess.get(v)+remcap); |
164 } |
180 exc-=remcap; |
165 } else { |
|
166 newlevel = newlevel < level.get(v) ? newlevel : level.get(v); |
|
167 } |
181 } |
168 |
182 } else if ( newlevel > level.get(v) ) newlevel = level.get(v); |
169 } //if the out edge wv is in the res graph |
183 |
170 |
|
171 } //for out edges wv |
184 } //for out edges wv |
|
185 |
|
186 |
|
187 if ( exc > 0 ) { |
|
188 for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) { |
|
189 |
|
190 if( flow.get(e) == 0 ) continue; |
|
191 NodeIt v=G.tail(e); |
|
192 //e=vw |
|
193 |
|
194 if( lev > level.get(v) ) { |
|
195 /*Push is allowed now*/ |
|
196 |
|
197 if ( excess.get(v)==0 && v != s && v !=t) |
|
198 stack[level.get(v)].push(v); |
|
199 /*v becomes active.*/ |
|
200 |
|
201 int flo=flow.get(e); |
|
202 |
|
203 if ( flo >= exc ) { |
|
204 /*A nonsaturating push.*/ |
|
205 |
|
206 flow.set(e, flo-exc); |
|
207 excess.set(v, excess.get(v)+exc); |
|
208 exc=0; |
|
209 break; |
|
210 } else { |
|
211 /*A saturating push.*/ |
|
212 |
|
213 excess.set(v, excess.get(v)+flo); |
|
214 exc-=flo; |
|
215 flow.set(e,0); |
|
216 } |
|
217 } else if ( newlevel > level.get(v) ) newlevel = level.get(v); |
|
218 |
|
219 } //for in edges vw |
172 |
220 |
173 |
221 } // if w still has excess after the out edge for cycle |
174 if ( excess.get(w) > 0 ) { |
222 |
175 |
223 excess.set(w, exc); |
176 for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) { |
224 |
177 NodeIt v=G.tail(e); /*e is the edge vw.*/ |
225 |
178 |
226 |
179 if( flow.get(e) > 0 ) { |
227 |
180 /*e is an edge of the residual graph */ |
228 /* |
181 |
229 Relabel |
182 if( level.get(w)==level.get(v)+1 ) { |
230 */ |
183 /*Push is allowed now*/ |
231 |
184 |
232 if ( exc > 0 ) { |
185 if ( excess.get(v)==0 && v != s && v !=t) stack[level.get(v)].push(v); |
233 //now 'lev' is the old level of w |
186 /*v becomes active.*/ |
234 level.set(w,++newlevel); |
187 |
|
188 if ( flow.get(e) > excess.get(w) ) { |
|
189 /*A nonsaturating push.*/ |
|
190 |
|
191 flow.set(e, flow.get(e)-excess.get(w)); |
|
192 excess.set(v, excess.get(v)+excess.get(w)); |
|
193 excess.set(w,0); |
|
194 break; |
|
195 } else { |
|
196 /*A saturating push.*/ |
|
197 |
|
198 excess.set(v, excess.get(v)+flow.get(e)); |
|
199 excess.set(w, excess.get(w)-flow.get(e)); |
|
200 flow.set(e,0); |
|
201 if ( excess.get(w)==0 ) break; |
|
202 } |
|
203 } else { |
|
204 newlevel = newlevel < level.get(v) ? newlevel : level.get(v); |
|
205 } |
|
206 |
|
207 } //if in edge vw is in the res graph |
|
208 |
|
209 } //for in edges vw |
|
210 |
|
211 } // if w still has excess after the out edge for cycle |
|
212 |
|
213 |
|
214 /* |
|
215 Relabel |
|
216 */ |
|
217 |
235 |
218 if ( excess.get(w) > 0 ) { |
236 if ( lev < n ) { |
219 |
237 --numb[lev]; |
220 int oldlevel=level.get(w); |
238 |
221 level.set(w,++newlevel); |
239 if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty. |
222 |
240 |
223 if ( oldlevel < n ) { |
241 for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) { |
224 --numb[oldlevel]; |
242 if (level.get(v) > lev && level.get(v) < n ) level.set(v,n); |
225 |
|
226 if ( !numb[oldlevel] && oldlevel < A*n ) { //If the level of w gets empty. |
|
227 |
|
228 for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) { |
|
229 if (level.get(v) > oldlevel && level.get(v) < n ) level.set(v,n); |
|
230 } |
|
231 for (int i=oldlevel+1 ; i!=n ; ++i) numb[i]=0; |
|
232 if ( newlevel < n ) newlevel=n; |
|
233 } else { |
|
234 if ( newlevel < n ) ++numb[newlevel]; |
|
235 } |
243 } |
|
244 for (int i=lev+1 ; i!=n ; ++i) numb[i]=0; |
|
245 if ( newlevel < n ) newlevel=n; |
236 } else { |
246 } else { |
237 if ( newlevel < n ) ++numb[newlevel]; |
247 if ( newlevel < n ) ++numb[newlevel]; |
238 } |
248 } |
239 |
249 } |
240 stack[newlevel].push(w); |
250 |
241 b=newlevel; |
251 stack[newlevel].push(w); |
242 |
252 b=newlevel; |
243 } |
253 |
244 |
254 } |
245 } // if stack[b] is nonempty |
255 |
246 |
|
247 } // while(b) |
256 } // while(b) |
248 |
257 |
249 |
258 |
250 value = excess.get(t); |
259 value = excess.get(t); |
251 /*Max flow value.*/ |
260 /*Max flow value.*/ |
252 |
261 |
253 |
262 |
254 } //void run() |
263 } //void run() |
269 |
278 |
270 /* |
279 /* |
271 For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). |
280 For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). |
272 */ |
281 */ |
273 |
282 |
274 T flowonedge(EdgeIt e) { |
283 T flowonedge(const EdgeIt e) { |
275 return flow.get(e); |
284 return flow.get(e); |
276 } |
285 } |
277 |
286 |
278 |
287 |
279 |
288 |
280 /* |
289 /* |
281 Returns the maximum flow x found by the algorithm. |
290 Returns the maximum flow x found by the algorithm. |
282 */ |
291 */ |
283 |
292 |
284 typename Graph::EdgeMap<T> allflow() { |
293 FlowMap allflow() { |
285 return flow; |
294 return flow; |
286 } |
295 } |
287 |
296 |
288 |
297 |
289 |
298 |
290 /* |
299 |
291 Returns a minimum cut by using a reverse bfs from t in the residual graph. |
300 /* |
|
301 Returns the minimum min cut, by a bfs from s in the residual graph. |
292 */ |
302 */ |
293 |
303 |
294 typename Graph::NodeMap<bool> mincut() { |
304 template<typename CutMap> |
|
305 void mincut(CutMap& M) { |
295 |
306 |
296 std::queue<NodeIt> queue; |
307 std::queue<NodeIt> queue; |
297 |
308 |
298 mincutvector.set(t,false); |
309 M.set(s,true); |
299 queue.push(t); |
310 queue.push(s); |
300 |
311 |
301 while (!queue.empty()) { |
312 while (!queue.empty()) { |
302 NodeIt w=queue.front(); |
313 NodeIt w=queue.front(); |
303 queue.pop(); |
314 queue.pop(); |
|
315 |
|
316 for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) { |
|
317 NodeIt v=G.head(e); |
|
318 if (!M.get(v) && flow.get(e) < capacity.get(e) ) { |
|
319 queue.push(v); |
|
320 M.set(v, true); |
|
321 } |
|
322 } |
304 |
323 |
305 for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) { |
324 for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) { |
306 NodeIt v=G.tail(e); |
325 NodeIt v=G.tail(e); |
307 if (mincutvector.get(v) && flow.get(e) < capacity.get(e) ) { |
326 if (!M.get(v) && flow.get(e) > 0 ) { |
308 queue.push(v); |
327 queue.push(v); |
309 mincutvector.set(v, false); |
328 M.set(v, true); |
310 } |
329 } |
311 } // for |
330 } |
|
331 |
|
332 } |
|
333 } |
|
334 |
|
335 |
|
336 |
|
337 /* |
|
338 Returns the maximum min cut, by a reverse bfs |
|
339 from t in the residual graph. |
|
340 */ |
|
341 |
|
342 template<typename CutMap> |
|
343 void max_mincut(CutMap& M) { |
|
344 |
|
345 std::queue<NodeIt> queue; |
|
346 |
|
347 M.set(t,true); |
|
348 queue.push(t); |
|
349 |
|
350 while (!queue.empty()) { |
|
351 NodeIt w=queue.front(); |
|
352 queue.pop(); |
|
353 |
|
354 for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) { |
|
355 NodeIt v=G.tail(e); |
|
356 if (!M.get(v) && flow.get(e) < capacity.get(e) ) { |
|
357 queue.push(v); |
|
358 M.set(v, true); |
|
359 } |
|
360 } |
312 |
361 |
313 for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) { |
362 for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) { |
314 NodeIt v=G.head(e); |
363 NodeIt v=G.head(e); |
315 if (mincutvector.get(v) && flow.get(e) > 0 ) { |
364 if (!M.get(v) && flow.get(e) > 0 ) { |
316 queue.push(v); |
365 queue.push(v); |
317 mincutvector.set(v, false); |
366 M.set(v, true); |
318 } |
367 } |
319 } // for |
368 } |
320 |
|
321 } |
369 } |
322 |
370 |
323 return mincutvector; |
371 for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) { |
324 |
372 M.set(v, !M.get(v)); |
325 } |
373 } |
|
374 |
|
375 } |
|
376 |
|
377 |
|
378 |
|
379 template<typename CutMap> |
|
380 void min_mincut(CutMap& M) { |
|
381 mincut(M); |
|
382 } |
|
383 |
|
384 |
|
385 |
326 }; |
386 }; |
327 }//namespace marci |
387 }//namespace marci |
328 #endif |
388 #endif |
329 |
389 |
330 |
390 |