13 |
14 |
14 The following functions should be used after run() was already run. |
15 The following functions should be used after run() was already run. |
15 |
16 |
16 T maxflow() : returns the value of a maximum flow |
17 T maxflow() : returns the value of a maximum flow |
17 |
18 |
18 NodeMap<bool> mincut(): returns a |
19 void mincut(CutMap& M) : sets M to the characteristic vector of a |
19 characteristic vector of a minimum cut. |
20 minimum cut. M should be a map of bools initialized to false. |
|
21 |
20 */ |
22 */ |
21 |
23 |
22 #ifndef PREFLOW_PUSH_MAX_FLOW_H |
24 #ifndef PREFLOW_PUSH_MAX_FLOW_H |
23 #define PREFLOW_PUSH_MAX_FLOW_H |
25 #define PREFLOW_PUSH_MAX_FLOW_H |
|
26 |
|
27 #define A 1 |
24 |
28 |
25 #include <algorithm> |
29 #include <algorithm> |
26 #include <vector> |
30 #include <vector> |
27 #include <stack> |
31 #include <stack> |
28 |
32 |
29 #include <reverse_bfs.h> |
33 #include <reverse_bfs.h> |
30 |
34 |
31 |
35 |
32 namespace marci { |
36 namespace marci { |
33 |
37 |
34 template <typename Graph, typename T> |
38 template <typename Graph, typename T, |
|
39 typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>, |
|
40 typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> > |
35 class preflow_push_max_flow { |
41 class preflow_push_max_flow { |
36 |
42 |
37 typedef typename Graph::NodeIt NodeIt; |
43 typedef typename Graph::NodeIt NodeIt; |
38 typedef typename Graph::EachNodeIt EachNodeIt; |
44 typedef typename Graph::EachNodeIt EachNodeIt; |
39 typedef typename Graph::OutEdgeIt OutEdgeIt; |
45 typedef typename Graph::OutEdgeIt OutEdgeIt; |
40 typedef typename Graph::InEdgeIt InEdgeIt; |
46 typedef typename Graph::InEdgeIt InEdgeIt; |
41 |
47 |
42 Graph& G; |
48 Graph& G; |
43 NodeIt s; |
49 NodeIt s; |
44 NodeIt t; |
50 NodeIt t; |
45 typename Graph::EdgeMap<T>& capacity; |
51 IntMap level; |
46 T value; |
52 CapMap& capacity; |
47 typename Graph::NodeMap<bool> mincutvector; |
53 int empty_level; //an empty level in the end of run() |
48 |
54 T value; |
49 |
55 |
50 |
|
51 public: |
56 public: |
52 |
57 |
53 preflow_push_max_flow ( Graph& _G, NodeIt _s, NodeIt _t, |
58 preflow_push_max_flow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) : |
54 typename Graph::EdgeMap<T>& _capacity) : |
59 G(_G), s(_s), t(_t), level(_G), capacity(_capacity) { } |
55 G(_G), s(_s), t(_t), capacity(_capacity), mincutvector(_G, false) { } |
|
56 |
60 |
57 |
61 |
58 /* |
62 /* |
59 The run() function runs a modified version of the |
63 The run() function runs a modified version of the |
60 highest label preflow-push, which only |
64 highest label preflow-push, which only |
61 finds a maximum preflow, hence giving the value of a maximum flow. |
65 finds a maximum preflow, hence giving the value of a maximum flow. |
62 */ |
66 */ |
63 void run() { |
67 void run() { |
64 |
68 |
65 typename Graph::EdgeMap<T> flow(G, 0); |
69 int n=G.nodeNum(); |
66 typename Graph::NodeMap<int> level(G); |
|
67 typename Graph::NodeMap<T> excess(G); |
|
68 |
|
69 int n=G.nodeNum(); |
|
70 int b=n-2; |
70 int b=n-2; |
71 /* |
71 /* |
72 b is a bound on the highest level of an active Node. |
72 b is a bound on the highest level of an active node. |
73 In the beginning it is at most n-2. |
73 */ |
74 */ |
74 |
75 |
75 IntMap level(G,n); |
76 std::vector<int> numb(n); //The number of Nodes on level i < n. |
76 TMap excess(G); |
77 std::vector<std::stack<NodeIt> > stack(2*n-1); |
77 FlowMap flow(G,0); |
78 //Stack of the active Nodes in level i. |
78 |
|
79 std::vector<int> numb(n); |
|
80 /* |
|
81 The number of nodes on level i < n. It is |
|
82 initialized to n+1, because of the reverse_bfs-part. |
|
83 */ |
|
84 |
|
85 std::vector<std::stack<NodeIt> > stack(n); |
|
86 //Stack of the active nodes in level i. |
|
87 |
79 |
88 |
80 /*Reverse_bfs from t, to find the starting level.*/ |
89 /*Reverse_bfs from t, to find the starting level.*/ |
81 reverse_bfs<Graph> bfs(G, t); |
90 level.set(t,0); |
82 bfs.run(); |
91 std::queue<NodeIt> bfs_queue; |
83 for(EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) |
92 bfs_queue.push(t); |
84 { |
93 |
85 int dist=bfs.dist(v); |
94 while (!bfs_queue.empty()) { |
86 level.set(v, dist); |
95 |
87 ++numb[dist]; |
96 NodeIt v=bfs_queue.front(); |
|
97 bfs_queue.pop(); |
|
98 int l=level.get(v)+1; |
|
99 |
|
100 for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) { |
|
101 NodeIt w=G.tail(e); |
|
102 if ( level.get(w) == n ) { |
|
103 bfs_queue.push(w); |
|
104 ++numb[l]; |
|
105 level.set(w, l); |
|
106 } |
88 } |
107 } |
89 |
108 } |
|
109 |
90 level.set(s,n); |
110 level.set(s,n); |
91 |
111 |
92 /* Starting flow. It is everywhere 0 at the moment. */ |
112 |
|
113 |
|
114 /* Starting flow. It is everywhere 0 at the moment. */ |
93 for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) |
115 for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) |
94 { |
116 { |
95 if ( capacity.get(e) > 0 ) { |
117 if ( capacity.get(e) == 0 ) continue; |
96 NodeIt w=G.head(e); |
118 NodeIt w=G.head(e); |
|
119 if ( level.get(w) < n ) { |
|
120 if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); |
97 flow.set(e, capacity.get(e)); |
121 flow.set(e, capacity.get(e)); |
98 stack[level.get(w)].push(w); |
|
99 excess.set(w, excess.get(w)+capacity.get(e)); |
122 excess.set(w, excess.get(w)+capacity.get(e)); |
100 } |
123 } |
101 } |
124 } |
102 |
125 |
103 /* |
126 /* |
104 End of preprocessing |
127 End of preprocessing |
105 */ |
128 */ |
106 |
129 |
107 |
130 |
108 |
|
109 /* |
131 /* |
110 Push/relabel on the highest level active Nodes. |
132 Push/relabel on the highest level active nodes. |
111 */ |
133 */ |
112 |
134 /*While there exists an active node.*/ |
113 /*While there exists an active Node.*/ |
|
114 while (b) { |
135 while (b) { |
115 |
136 if ( stack[b].empty() ) { |
116 /*We decrease the bound if there is no active node of level b.*/ |
|
117 if (stack[b].empty()) { |
|
118 --b; |
137 --b; |
119 } else { |
138 continue; |
120 |
139 } |
121 NodeIt w=stack[b].top(); //w is the highest label active Node. |
140 |
122 stack[b].pop(); //We delete w from the stack. |
141 NodeIt w=stack[b].top(); //w is a highest label active node. |
123 |
142 stack[b].pop(); |
124 int newlevel=2*n-2; //In newlevel we maintain the next level of w. |
143 int lev=level.get(w); |
125 |
144 int exc=excess.get(w); |
|
145 int newlevel=2*n-2; //In newlevel we bound the next level of w. |
|
146 |
|
147 // if ( level.get(w) < n ) { //Nem tudom ez mukodik-e |
126 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) { |
148 for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) { |
127 NodeIt v=G.head(e); |
149 |
128 /*e is the Edge wv.*/ |
150 if ( flow.get(e) == capacity.get(e) ) continue; |
129 |
151 NodeIt v=G.head(e); |
130 if (flow.get(e)<capacity.get(e)) { |
152 //e=wv |
131 /*e is an Edge of the residual graph */ |
153 |
132 |
154 if( lev > level.get(v) ) { |
133 if(level.get(w)==level.get(v)+1) { |
155 /*Push is allowed now*/ |
134 /*Push is allowed now*/ |
156 |
135 |
157 if ( excess.get(v)==0 && v != s && v !=t ) |
136 if (capacity.get(e)-flow.get(e) > excess.get(w)) { |
158 stack[level.get(v)].push(v); |
137 /*A nonsaturating push.*/ |
159 /*v becomes active.*/ |
138 |
160 |
139 if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v); |
161 int cap=capacity.get(e); |
140 /*v becomes active.*/ |
162 int flo=flow.get(e); |
141 |
163 int remcap=cap-flo; |
142 flow.set(e, flow.get(e)+excess.get(w)); |
164 |
143 excess.set(v, excess.get(v)+excess.get(w)); |
165 if ( remcap >= exc ) { |
144 excess.set(w,0); |
166 /*A nonsaturating push.*/ |
145 //std::cout << w << " " << v <<" elore elen nonsat pump " << std::endl; |
167 |
146 break; |
168 flow.set(e, flo+exc); |
147 } else { |
169 excess.set(v, excess.get(v)+exc); |
148 /*A saturating push.*/ |
170 exc=0; |
149 |
171 break; |
150 if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v); |
172 |
151 /*v becomes active.*/ |
173 } else { |
152 |
174 /*A saturating push.*/ |
153 excess.set(v, excess.get(v)+capacity.get(e)-flow.get(e)); |
175 |
154 excess.set(w, excess.get(w)-capacity.get(e)+flow.get(e)); |
176 flow.set(e, cap ); |
155 flow.set(e, capacity.get(e)); |
177 excess.set(v, excess.get(v)+remcap); |
156 //std::cout << w <<" " << v <<" elore elen sat pump " << std::endl; |
178 exc-=remcap; |
157 if (excess.get(w)==0) break; |
179 } |
158 /*If w is not active any more, then we go on to the next Node.*/ |
180 } else if ( newlevel > level.get(v) ) newlevel = level.get(v); |
159 |
181 |
160 } // if (capacity.get(e)-flow.get(e) > excess.get(w)) |
182 } //for out edges wv |
161 } // if (level.get(w)==level.get(v)+1) |
183 |
162 |
184 |
163 else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);} |
185 if ( exc > 0 ) { |
164 |
186 for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) { |
165 } //if (flow.get(e)<capacity.get(e)) |
187 |
166 |
188 if( flow.get(e) == 0 ) continue; |
167 } //for(OutEdgeIt e=G.first_OutEdge(w); e.valid(); ++e) |
189 NodeIt v=G.tail(e); |
|
190 //e=vw |
|
191 |
|
192 if( lev > level.get(v) ) { |
|
193 /*Push is allowed now*/ |
|
194 |
|
195 if ( excess.get(v)==0 && v != s && v !=t) |
|
196 stack[level.get(v)].push(v); |
|
197 /*v becomes active.*/ |
|
198 |
|
199 int flo=flow.get(e); |
|
200 |
|
201 if ( flo >= exc ) { |
|
202 /*A nonsaturating push.*/ |
|
203 |
|
204 flow.set(e, flo-exc); |
|
205 excess.set(v, excess.get(v)+exc); |
|
206 exc=0; |
|
207 break; |
|
208 } else { |
|
209 /*A saturating push.*/ |
|
210 |
|
211 excess.set(v, excess.get(v)+flo); |
|
212 exc-=flo; |
|
213 flow.set(e,0); |
|
214 } |
|
215 } else if ( newlevel > level.get(v) ) newlevel = level.get(v); |
|
216 |
|
217 } //for in edges vw |
168 |
218 |
169 |
219 } // if w still has excess after the out edge for cycle |
170 |
220 |
171 for(InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) { |
221 excess.set(w, exc); |
172 NodeIt v=G.tail(e); |
222 |
173 /*e is the Edge vw.*/ |
223 |
174 |
224 /* |
175 if (excess.get(w)==0) break; |
225 Relabel |
176 /*It may happen, that w became inactive in the first 'for' cycle.*/ |
226 */ |
177 |
227 |
178 if(flow.get(e)>0) { |
228 if ( exc > 0 ) { |
179 /*e is an Edge of the residual graph */ |
229 //now 'lev' is the old level of w |
180 |
230 level.set(w,++newlevel); |
181 if(level.get(w)==level.get(v)+1) { |
231 --numb[lev]; |
182 /*Push is allowed now*/ |
232 |
183 |
233 if ( !numb[lev] && lev < A*n ) { //If the level of w gets empty. |
184 if (flow.get(e) > excess.get(w)) { |
234 |
185 /*A nonsaturating push.*/ |
235 for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) { |
186 |
236 if (level.get(v) > lev ) level.set(v,n); |
187 if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v); |
237 } |
188 /*v becomes active.*/ |
238 for (int i=lev+1 ; i!=n ; ++i) numb[i]=0; |
189 |
239 if ( newlevel < n ) newlevel=n; |
190 flow.set(e, flow.get(e)-excess.get(w)); |
240 } else if ( newlevel < n ) { |
191 excess.set(v, excess.get(v)+excess.get(w)); |
241 ++numb[newlevel]; |
192 excess.set(w,0); |
242 stack[newlevel].push(w); |
193 //std::cout << v << " " << w << " vissza elen nonsat pump " << std::endl; |
243 b=newlevel; |
194 break; |
244 } |
195 } else { |
245 } |
196 /*A saturating push.*/ |
246 |
197 |
247 |
198 if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v); |
248 |
199 /*v becomes active.*/ |
|
200 |
|
201 flow.set(e,0); |
|
202 excess.set(v, excess.get(v)+flow.get(e)); |
|
203 excess.set(w, excess.get(w)-flow.get(e)); |
|
204 //std::cout << v <<" " << w << " vissza elen sat pump " << std::endl; |
|
205 if (excess.get(w)==0) { break;} |
|
206 } //if (flow.get(e) > excess.get(v)) |
|
207 } //if(level.get(w)==level.get(v)+1) |
|
208 |
|
209 else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);} |
|
210 //std::cout << "Leveldecrease of Node " << w << " to " << newlevel << std::endl; |
|
211 |
|
212 } //if (flow.get(e)>0) |
|
213 |
|
214 } //for in-Edge |
|
215 |
|
216 |
|
217 |
|
218 |
|
219 /* |
|
220 Relabel |
|
221 */ |
|
222 if (excess.get(w)>0) { |
|
223 /*Now newlevel <= n*/ |
|
224 |
|
225 int l=level.get(w); //l is the old level of w. |
|
226 --numb[l]; |
|
227 |
|
228 if (newlevel == n) { |
|
229 level.set(w,n); |
|
230 |
|
231 } else { |
|
232 |
|
233 if (numb[l]) { |
|
234 /*If the level of w remains nonempty.*/ |
|
235 |
|
236 level.set(w,++newlevel); |
|
237 ++numb[newlevel]; |
|
238 stack[newlevel].push(w); |
|
239 b=newlevel; |
|
240 } else { |
|
241 /*If the level of w gets empty.*/ |
|
242 |
|
243 for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) { |
|
244 if (level.get(v) >= l ) { |
|
245 level.set(v,n); |
|
246 } |
|
247 } |
|
248 |
|
249 for (int i=l+1 ; i!=n ; ++i) numb[i]=0; |
|
250 } //if (numb[l]) |
|
251 |
|
252 } // if (newlevel = n) |
|
253 |
|
254 } // if (excess.get(w)>0) |
|
255 |
|
256 |
|
257 } //else |
|
258 |
|
259 } //while(b) |
249 } //while(b) |
260 |
250 |
261 value=excess.get(t); |
251 value=excess.get(t); |
262 /*Max flow value.*/ |
252 /*Max flow value.*/ |
263 |
253 |
264 |
254 |
265 |
255 /* |
266 /* |
256 We count empty_level. The nodes above this level is a mincut. |
267 We find an empty level, e. The Nodes above this level give |
257 */ |
268 a minimum cut. |
258 while(true) { |
269 */ |
259 if(numb[empty_level]) ++empty_level; |
270 |
|
271 int e=1; |
|
272 |
|
273 while(e) { |
|
274 if(numb[e]) ++e; |
|
275 else break; |
260 else break; |
276 } |
261 } |
277 for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid(); ++v) { |
262 |
278 if (level.get(v) > e) mincutvector.set(v, true); |
|
279 } |
|
280 |
|
281 |
|
282 } // void run() |
263 } // void run() |
283 |
264 |
284 |
265 |
285 |
266 |
286 /* |
267 /* |