1 // -*- C++ -*- |
1 // -*- C++ -*- |
2 |
2 |
3 /* |
3 /* |
4 Heuristics: |
4 Heuristics: |
5 2 phase |
5 2 phase |
6 gap |
6 gap |
7 list 'level_list' on the nodes on level i implemented by hand |
7 list 'level_list' on the nodes on level i implemented by hand |
8 stack 'active' on the active nodes on level i |
8 stack 'active' on the active nodes on level i |
9 runs heuristic 'highest label' for H1*n relabels |
9 runs heuristic 'highest label' for H1*n relabels |
10 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label' |
10 runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label' |
11 |
11 |
12 Parameters H0 and H1 are initialized to 20 and 1. |
12 Parameters H0 and H1 are initialized to 20 and 1. |
13 |
13 |
14 Constructors: |
14 Constructors: |
15 |
15 |
16 Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if |
16 Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if |
17 FlowMap is not constant zero, and should be true if it is |
17 FlowMap is not constant zero, and should be true if it is |
18 |
18 |
19 Members: |
19 Members: |
20 |
20 |
21 void run() |
21 void run() |
22 |
22 |
23 Num flowValue() : returns the value of a maximum flow |
23 Num flowValue() : returns the value of a maximum flow |
24 |
24 |
25 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 |
26 minimum min cut. M should be a map of bools initialized to false. ??Is it OK? |
26 minimum min cut. M should be a map of bools initialized to false. ??Is it OK? |
27 |
27 |
28 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 |
29 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. |
30 |
30 |
31 void minCut(CutMap& M) : sets M to the characteristic vector of |
31 void minCut(CutMap& M) : sets M to the characteristic vector of |
32 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. |
33 |
33 |
34 */ |
34 */ |
35 |
35 |
36 #ifndef HUGO_MAX_FLOW_H |
36 #ifndef HUGO_MAX_FLOW_H |
37 #define HUGO_MAX_FLOW_H |
37 #define HUGO_MAX_FLOW_H |
38 |
|
39 #define H0 20 |
|
40 #define H1 1 |
|
41 |
38 |
42 #include <vector> |
39 #include <vector> |
43 #include <queue> |
40 #include <queue> |
44 #include <stack> |
41 #include <stack> |
45 |
42 |
72 typedef typename std::vector<Node> VecNode; |
71 typedef typename std::vector<Node> VecNode; |
73 |
72 |
74 const Graph* g; |
73 const Graph* g; |
75 Node s; |
74 Node s; |
76 Node t; |
75 Node t; |
77 const CapMap* capacity; |
76 const CapMap* capacity; |
78 FlowMap* flow; |
77 FlowMap* flow; |
79 int n; //the number of nodes of G |
78 int n; //the number of nodes of G |
80 typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW; |
79 typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW; |
81 typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt; |
80 typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt; |
82 typedef typename ResGW::Edge ResGWEdge; |
81 typedef typename ResGW::Edge ResGWEdge; |
83 //typedef typename ResGW::template NodeMap<bool> ReachedMap; |
82 //typedef typename ResGW::template NodeMap<bool> ReachedMap; |
84 typedef typename Graph::template NodeMap<int> ReachedMap; |
83 typedef typename Graph::template NodeMap<int> ReachedMap; |
85 ReachedMap level; |
84 ReachedMap level; |
86 //level works as a bool map in augmenting path algorithms |
85 //level works as a bool map in augmenting path algorithms |
87 //and is used by bfs for storing reached information. |
86 //and is used by bfs for storing reached information. |
88 //In preflow, it shows levels of nodes. |
87 //In preflow, it shows levels of nodes. |
89 //typename Graph::template NodeMap<int> level; |
88 //typename Graph::template NodeMap<int> level; |
90 typename Graph::template NodeMap<Num> excess; |
89 typename Graph::template NodeMap<Num> excess; |
91 // protected: |
90 // protected: |
92 // MaxFlow() { } |
91 // MaxFlow() { } |
93 // void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, |
92 // void set(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, |
94 // FlowMap& _flow) |
93 // FlowMap& _flow) |
95 // { |
94 // { |
96 // g=&_G; |
95 // g=&_G; |
97 // s=_s; |
96 // s=_s; |
98 // t=_t; |
97 // t=_t; |
99 // capacity=&_capacity; |
98 // capacity=&_capacity; |
100 // flow=&_flow; |
99 // flow=&_flow; |
101 // n=_G.nodeNum; |
100 // n=_G.nodeNum; |
102 // level.set (_G); //kellene vmi ilyesmi fv |
101 // level.set (_G); //kellene vmi ilyesmi fv |
103 // excess(_G,0); //itt is |
102 // excess(_G,0); //itt is |
104 // } |
103 // } |
105 |
104 |
|
105 // constants used for heuristics |
|
106 static const int H0=20; |
|
107 static const int H1=1; |
|
108 |
106 public: |
109 public: |
107 |
110 |
108 ///\todo Document this. |
111 ///\todo Document this. |
109 ///\todo Maybe, it should be PRE_FLOW instead. |
112 ///\todo Maybe, it should be PRE_FLOW instead. |
|
113 ///- \c NO_FLOW means nothing, |
110 ///- \c ZERO_FLOW means something, |
114 ///- \c ZERO_FLOW means something, |
111 ///- \c GEN_FLOW means something else, |
115 ///- \c GEN_FLOW means something else, |
112 ///- \c PREFLOW is something different. |
116 ///- \c PRE_FLOW is something different. |
113 enum flowEnum{ |
117 enum flowEnum{ |
114 ZERO_FLOW=0, |
118 ZERO_FLOW, |
115 GEN_FLOW=1, |
119 GEN_FLOW, |
116 PREFLOW=2 |
120 PRE_FLOW, |
|
121 NO_FLOW |
117 }; |
122 }; |
118 |
123 |
119 MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, |
124 MaxFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, |
120 FlowMap& _flow) : |
125 FlowMap& _flow) : |
121 g(&_G), s(_s), t(_t), capacity(&_capacity), |
126 g(&_G), s(_s), t(_t), capacity(&_capacity), |
122 flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {} |
127 flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {} |
123 |
128 |
124 /// A max flow algorithm is run. |
129 /// A max flow algorithm is run. |
125 ///\pre the flow have to be 0 at the beginning. |
130 /// \pre The flow have to satisfy the requirements |
126 void run() { |
131 /// stated in fe. |
127 preflow(ZERO_FLOW); |
132 void run(flowEnum fe=ZERO_FLOW) { |
128 } |
133 preflow(fe); |
129 |
134 } |
130 /// A preflow algorithm is run. |
135 |
131 ///\pre The initial edge-map have to be a |
136 /// A preflow algorithm is run. |
|
137 /// \pre The initial edge-map have to be a |
132 /// zero flow if \c fe is \c ZERO_FLOW, |
138 /// zero flow if \c fe is \c ZERO_FLOW, |
133 /// a flow if \c fe is \c GEN_FLOW, |
139 /// a flow if \c fe is \c GEN_FLOW, |
134 /// and a pre-flow it is \c PREFLOW. |
140 /// a pre-flow if fe is \c PRE_FLOW and |
|
141 /// anything if fe is NO_FLOW. |
135 void preflow(flowEnum fe) { |
142 void preflow(flowEnum fe) { |
136 preflowPhase0(fe); |
143 preflowPhase0(fe); |
137 preflowPhase1(); |
144 preflowPhase1(); |
138 } |
145 } |
139 |
146 |
140 /// Run the first phase of preflow, starting from a 0 flow, from a flow, |
147 /// Run the first phase of preflow, starting from a 0 flow, from a flow, |
141 /// or from a preflow, according to \c fe. |
148 /// or from a preflow, of from undefined value according to \c fe. |
142 void preflowPhase0( flowEnum fe ); |
149 void preflowPhase0(flowEnum fe); |
143 |
150 |
144 /// Second phase of preflow. |
151 /// Second phase of preflow. |
145 void preflowPhase1(); |
152 void preflowPhase1(); |
146 |
153 |
147 /// Starting from a flow, this method searches for an augmenting path |
154 /// Starting from a flow, this method searches for an augmenting path |
148 /// according to the Edmonds-Karp algorithm |
155 /// according to the Edmonds-Karp algorithm |
149 /// and augments the flow on if any. |
156 /// and augments the flow on if any. |
150 /// The return value shows if the augmentation was succesful. |
157 /// The return value shows if the augmentation was succesful. |
151 bool augmentOnShortestPath(); |
158 bool augmentOnShortestPath(); |
152 |
159 |
153 /// Starting from a flow, this method searches for an augmenting blockin |
160 /// Starting from a flow, this method searches for an augmenting blocking |
154 /// flow according to Dinits' algorithm and augments the flow on if any. |
161 /// flow according to Dinits' algorithm and augments the flow on if any. |
155 /// The blocking flow is computed in a physically constructed |
162 /// The blocking flow is computed in a physically constructed |
156 /// residual graph of type \c Mutablegraph. |
163 /// residual graph of type \c Mutablegraph. |
157 /// The return value show sif the augmentation was succesful. |
164 /// The return value show sif the augmentation was succesful. |
158 template<typename MutableGraph> bool augmentOnBlockingFlow(); |
165 template<typename MutableGraph> bool augmentOnBlockingFlow(); |
159 |
166 |
160 /// The same as \c augmentOnBlockingFlow<MutableGraph> but the |
167 /// The same as \c augmentOnBlockingFlow<MutableGraph> but the |
161 /// residual graph is not constructed physically. |
168 /// residual graph is not constructed physically. |
162 /// The return value shows if the augmentation was succesful. |
169 /// The return value shows if the augmentation was succesful. |
163 bool augmentOnBlockingFlow2(); |
170 bool augmentOnBlockingFlow2(); |
164 |
171 |
165 /// Returns the actual flow value. |
172 /// Returns the actual flow value. |
166 /// More precisely, it returns the negative excess of s, thus |
173 /// More precisely, it returns the negative excess of s, thus |
167 /// this works also for preflows. |
174 /// this works also for preflows. |
168 Num flowValue() { |
175 Num flowValue() { |
169 Num a=0; |
176 Num a=0; |
170 FOR_EACH_INC_LOC(OutEdgeIt, e, *g, s) a+=(*flow)[e]; |
177 FOR_EACH_INC_LOC(InEdgeIt, e, *g, t) a+=(*flow)[e]; |
171 FOR_EACH_INC_LOC(InEdgeIt, e, *g, s) a-=(*flow)[e]; |
178 FOR_EACH_INC_LOC(OutEdgeIt, e, *g, t) a-=(*flow)[e]; |
172 return a; |
179 return a; |
173 } |
180 } |
174 |
181 |
175 /// Should be used between preflowPhase0 and preflowPhase1. |
182 /// Should be used between preflowPhase0 and preflowPhase1. |
176 ///\todo We have to make some status variable which shows the actual state |
183 /// \todo We have to make some status variable which shows the |
177 /// of the class. This enables us to determine which methods are valid |
184 /// actual state |
|
185 /// of the class. This enables us to determine which methods are valid |
178 /// for MinCut computation |
186 /// for MinCut computation |
179 template<typename _CutMap> |
187 template<typename _CutMap> |
180 void actMinCut(_CutMap& M) { |
188 void actMinCut(_CutMap& M) { |
181 NodeIt v; |
189 NodeIt v; |
182 for(g->first(v); g->valid(v); g->next(v)) { |
190 for(g->first(v); g->valid(v); g->next(v)) { |
272 |
279 |
273 /// |
280 /// |
274 void resetSource(Node _s) { s=_s; } |
281 void resetSource(Node _s) { s=_s; } |
275 /// |
282 /// |
276 void resetTarget(Node _t) { t=_t; } |
283 void resetTarget(Node _t) { t=_t; } |
277 |
284 |
278 /// capacity-map is changed. |
285 /// capacity-map is changed. |
279 void resetCap(const CapMap& _cap) { capacity=&_cap; } |
286 void resetCap(const CapMap& _cap) { capacity=&_cap; } |
280 |
287 |
281 /// flow-map is changed. |
288 /// flow-map is changed. |
282 void resetFlow(FlowMap& _flow) { flow=&_flow; } |
289 void resetFlow(FlowMap& _flow) { flow=&_flow; } |
283 |
290 |
284 |
291 |
285 private: |
292 private: |
286 |
293 |
287 int push(Node w, VecStack& active) { |
294 int push(Node w, VecStack& active) { |
288 |
295 |
289 int lev=level[w]; |
296 int lev=level[w]; |
290 Num exc=excess[w]; |
297 Num exc=excess[w]; |
291 int newlevel=n; //bound on the next level of w |
298 int newlevel=n; //bound on the next level of w |
292 |
299 |
293 OutEdgeIt e; |
300 OutEdgeIt e; |
294 for(g->first(e,w); g->valid(e); g->next(e)) { |
301 for(g->first(e,w); g->valid(e); g->next(e)) { |
295 |
302 |
296 if ( (*flow)[e] >= (*capacity)[e] ) continue; |
303 if ( (*flow)[e] >= (*capacity)[e] ) continue; |
297 Node v=g->head(e); |
304 Node v=g->head(e); |
298 |
305 |
299 if( lev > level[v] ) { //Push is allowed now |
306 if( lev > level[v] ) { //Push is allowed now |
300 |
307 |
301 if ( excess[v]<=0 && v!=t && v!=s ) { |
308 if ( excess[v]<=0 && v!=t && v!=s ) { |
302 int lev_v=level[v]; |
309 int lev_v=level[v]; |
303 active[lev_v].push(v); |
310 active[lev_v].push(v); |
304 } |
311 } |
305 |
312 |
306 Num cap=(*capacity)[e]; |
313 Num cap=(*capacity)[e]; |
307 Num flo=(*flow)[e]; |
314 Num flo=(*flow)[e]; |
308 Num remcap=cap-flo; |
315 Num remcap=cap-flo; |
309 |
316 |
310 if ( remcap >= exc ) { //A nonsaturating push. |
317 if ( remcap >= exc ) { //A nonsaturating push. |
311 |
318 |
312 flow->set(e, flo+exc); |
319 flow->set(e, flo+exc); |
313 excess.set(v, excess[v]+exc); |
320 excess.set(v, excess[v]+exc); |
314 exc=0; |
321 exc=0; |
315 break; |
322 break; |
316 |
323 |
317 } else { //A saturating push. |
324 } else { //A saturating push. |
318 flow->set(e, cap); |
325 flow->set(e, cap); |
319 excess.set(v, excess[v]+remcap); |
326 excess.set(v, excess[v]+remcap); |
320 exc-=remcap; |
327 exc-=remcap; |
321 } |
328 } |
322 } else if ( newlevel > level[v] ) newlevel = level[v]; |
329 } else if ( newlevel > level[v] ) newlevel = level[v]; |
323 } //for out edges wv |
330 } //for out edges wv |
324 |
331 |
325 if ( exc > 0 ) { |
332 if ( exc > 0 ) { |
326 InEdgeIt e; |
333 InEdgeIt e; |
327 for(g->first(e,w); g->valid(e); g->next(e)) { |
334 for(g->first(e,w); g->valid(e); g->next(e)) { |
328 |
335 |
329 if( (*flow)[e] <= 0 ) continue; |
336 if( (*flow)[e] <= 0 ) continue; |
330 Node v=g->tail(e); |
337 Node v=g->tail(e); |
331 |
338 |
332 if( lev > level[v] ) { //Push is allowed now |
339 if( lev > level[v] ) { //Push is allowed now |
333 |
340 |
334 if ( excess[v]<=0 && v!=t && v!=s ) { |
341 if ( excess[v]<=0 && v!=t && v!=s ) { |
335 int lev_v=level[v]; |
342 int lev_v=level[v]; |
336 active[lev_v].push(v); |
343 active[lev_v].push(v); |
337 } |
344 } |
338 |
345 |
339 Num flo=(*flow)[e]; |
346 Num flo=(*flow)[e]; |
340 |
347 |
341 if ( flo >= exc ) { //A nonsaturating push. |
348 if ( flo >= exc ) { //A nonsaturating push. |
342 |
349 |
343 flow->set(e, flo-exc); |
350 flow->set(e, flo-exc); |
344 excess.set(v, excess[v]+exc); |
351 excess.set(v, excess[v]+exc); |
345 exc=0; |
352 exc=0; |
346 break; |
353 break; |
347 } else { //A saturating push. |
354 } else { //A saturating push. |
348 |
355 |
349 excess.set(v, excess[v]+flo); |
356 excess.set(v, excess[v]+flo); |
350 exc-=flo; |
357 exc-=flo; |
351 flow->set(e,0); |
358 flow->set(e,0); |
352 } |
359 } |
353 } else if ( newlevel > level[v] ) newlevel = level[v]; |
360 } else if ( newlevel > level[v] ) newlevel = level[v]; |
354 } //for in edges vw |
361 } //for in edges vw |
355 |
362 |
356 } // if w still has excess after the out edge for cycle |
363 } // if w still has excess after the out edge for cycle |
357 |
364 |
358 excess.set(w, exc); |
365 excess.set(w, exc); |
359 |
366 |
360 return newlevel; |
367 return newlevel; |
361 } |
368 } |
362 |
369 |
363 |
370 |
364 void preflowPreproc ( flowEnum fe, VecStack& active, |
371 void preflowPreproc(flowEnum fe, VecStack& active, |
365 VecNode& level_list, NNMap& left, NNMap& right ) |
372 VecNode& level_list, NNMap& left, NNMap& right) |
366 { |
373 { |
367 |
|
368 std::queue<Node> bfs_queue; |
374 std::queue<Node> bfs_queue; |
369 |
375 |
370 switch ( fe ) { |
376 switch (fe) { |
371 case ZERO_FLOW: |
377 case ZERO_FLOW: |
372 { |
378 { |
373 //Reverse_bfs from t, to find the starting level. |
379 //Reverse_bfs from t, to find the starting level. |
374 level.set(t,0); |
380 level.set(t,0); |
375 bfs_queue.push(t); |
381 bfs_queue.push(t); |
376 |
382 |
377 while (!bfs_queue.empty()) { |
383 while (!bfs_queue.empty()) { |
378 |
384 |
379 Node v=bfs_queue.front(); |
385 Node v=bfs_queue.front(); |
380 bfs_queue.pop(); |
386 bfs_queue.pop(); |
381 int l=level[v]+1; |
387 int l=level[v]+1; |
382 |
388 |
383 InEdgeIt e; |
389 InEdgeIt e; |
384 for(g->first(e,v); g->valid(e); g->next(e)) { |
390 for(g->first(e,v); g->valid(e); g->next(e)) { |
385 Node w=g->tail(e); |
391 Node w=g->tail(e); |
386 if ( level[w] == n && w != s ) { |
392 if ( level[w] == n && w != s ) { |
387 bfs_queue.push(w); |
393 bfs_queue.push(w); |
450 level_list[l]=w; |
456 level_list[l]=w; |
451 level.set(w, l); |
457 level.set(w, l); |
452 } |
458 } |
453 } |
459 } |
454 } |
460 } |
455 |
461 |
456 |
462 |
457 //the starting flow |
463 //the starting flow |
458 OutEdgeIt e; |
464 OutEdgeIt e; |
459 for(g->first(e,s); g->valid(e); g->next(e)) |
465 for(g->first(e,s); g->valid(e); g->next(e)) |
460 { |
466 { |
461 Num rem=(*capacity)[e]-(*flow)[e]; |
467 Num rem=(*capacity)[e]-(*flow)[e]; |
462 if ( rem <= 0 ) continue; |
468 if ( rem <= 0 ) continue; |
463 Node w=g->head(e); |
469 Node w=g->head(e); |
464 if ( level[w] < n ) { |
470 if ( level[w] < n ) { |
465 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
471 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
466 flow->set(e, (*capacity)[e]); |
472 flow->set(e, (*capacity)[e]); |
467 excess.set(w, excess[w]+rem); |
473 excess.set(w, excess[w]+rem); |
468 } |
474 } |
469 } |
475 } |
470 |
476 |
471 InEdgeIt f; |
477 InEdgeIt f; |
472 for(g->first(f,s); g->valid(f); g->next(f)) |
478 for(g->first(f,s); g->valid(f); g->next(f)) |
473 { |
479 { |
474 if ( (*flow)[f] <= 0 ) continue; |
480 if ( (*flow)[f] <= 0 ) continue; |
475 Node w=g->tail(f); |
481 Node w=g->tail(f); |
476 if ( level[w] < n ) { |
482 if ( level[w] < n ) { |
477 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
483 if ( excess[w] <= 0 && w!=t ) active[level[w]].push(w); |
478 excess.set(w, excess[w]+(*flow)[f]); |
484 excess.set(w, excess[w]+(*flow)[f]); |
479 flow->set(f, 0); |
485 flow->set(f, 0); |
480 } |
486 } |
481 } |
487 } |
482 break; |
488 break; |
483 } //case PREFLOW |
489 } //case PRE_FLOW |
484 } |
490 } |
485 } //preflowPreproc |
491 } //preflowPreproc |
486 |
492 |
487 |
493 |
488 |
494 |
489 void relabel(Node w, int newlevel, VecStack& active, |
495 void relabel(Node w, int newlevel, VecStack& active, |
490 VecNode& level_list, NNMap& left, |
496 VecNode& level_list, NNMap& left, |
491 NNMap& right, int& b, int& k, bool what_heur ) |
497 NNMap& right, int& b, int& k, bool what_heur ) |
492 { |
498 { |
493 |
499 |
494 Num lev=level[w]; |
500 Num lev=level[w]; |
495 |
501 |
496 Node right_n=right[w]; |
502 Node right_n=right[w]; |
497 Node left_n=left[w]; |
503 Node left_n=left[w]; |
498 |
504 |
499 //unlacing starts |
505 //unlacing starts |
500 if ( g->valid(right_n) ) { |
506 if ( g->valid(right_n) ) { |
501 if ( g->valid(left_n) ) { |
507 if ( g->valid(left_n) ) { |
502 right.set(left_n, right_n); |
508 right.set(left_n, right_n); |
503 left.set(right_n, left_n); |
509 left.set(right_n, left_n); |
504 } else { |
510 } else { |
505 level_list[lev]=right_n; |
511 level_list[lev]=right_n; |
506 left.set(right_n, INVALID); |
512 left.set(right_n, INVALID); |
507 } |
513 } |
508 } else { |
514 } else { |
509 if ( g->valid(left_n) ) { |
515 if ( g->valid(left_n) ) { |
510 right.set(left_n, INVALID); |
516 right.set(left_n, INVALID); |
511 } else { |
517 } else { |
512 level_list[lev]=INVALID; |
518 level_list[lev]=INVALID; |
513 } |
519 } |
514 } |
520 } |
515 //unlacing ends |
521 //unlacing ends |
516 |
522 |
517 if ( !g->valid(level_list[lev]) ) { |
523 if ( !g->valid(level_list[lev]) ) { |
518 |
524 |
519 //gapping starts |
525 //gapping starts |
520 for (int i=lev; i!=k ; ) { |
526 for (int i=lev; i!=k ; ) { |
521 Node v=level_list[++i]; |
527 Node v=level_list[++i]; |
522 while ( g->valid(v) ) { |
528 while ( g->valid(v) ) { |
523 level.set(v,n); |
529 level.set(v,n); |
549 right.set(w,first); |
555 right.set(w,first); |
550 left.set(w,INVALID); |
556 left.set(w,INVALID); |
551 level_list[newlevel]=w; |
557 level_list[newlevel]=w; |
552 } |
558 } |
553 } |
559 } |
554 |
560 |
555 } //relabel |
561 } //relabel |
556 |
562 |
557 |
563 |
558 template<typename MapGraphWrapper> |
564 template<typename MapGraphWrapper> |
559 class DistanceMap { |
565 class DistanceMap { |
560 protected: |
566 protected: |
561 const MapGraphWrapper* g; |
567 const MapGraphWrapper* g; |
562 typename MapGraphWrapper::template NodeMap<int> dist; |
568 typename MapGraphWrapper::template NodeMap<int> dist; |
563 public: |
569 public: |
564 DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { } |
570 DistanceMap(MapGraphWrapper& _g) : g(&_g), dist(*g, g->nodeNum()) { } |
565 void set(const typename MapGraphWrapper::Node& n, int a) { |
571 void set(const typename MapGraphWrapper::Node& n, int a) { |
566 dist.set(n, a); |
572 dist.set(n, a); |
567 } |
573 } |
568 int operator[](const typename MapGraphWrapper::Node& n) |
574 int operator[](const typename MapGraphWrapper::Node& n) |
569 { return dist[n]; } |
575 { return dist[n]; } |
570 // int get(const typename MapGraphWrapper::Node& n) const { |
576 // int get(const typename MapGraphWrapper::Node& n) const { |
571 // return dist[n]; } |
577 // return dist[n]; } |
572 // bool get(const typename MapGraphWrapper::Edge& e) const { |
578 // bool get(const typename MapGraphWrapper::Edge& e) const { |
573 // return (dist.get(g->tail(e))<dist.get(g->head(e))); } |
579 // return (dist.get(g->tail(e))<dist.get(g->head(e))); } |
574 bool operator[](const typename MapGraphWrapper::Edge& e) const { |
580 bool operator[](const typename MapGraphWrapper::Edge& e) const { |
575 return (dist[g->tail(e)]<dist[g->head(e)]); |
581 return (dist[g->tail(e)]<dist[g->head(e)]); |
576 } |
582 } |
577 }; |
583 }; |
578 |
584 |
579 }; |
585 }; |
580 |
586 |
581 |
587 |
582 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
588 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
583 void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe ) |
589 void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase0( flowEnum fe ) |
584 { |
590 { |
585 |
591 |
586 int heur0=(int)(H0*n); //time while running 'bound decrease' |
592 int heur0=(int)(H0*n); //time while running 'bound decrease' |
587 int heur1=(int)(H1*n); //time while running 'highest label' |
593 int heur1=(int)(H1*n); //time while running 'highest label' |
588 int heur=heur1; //starting time interval (#of relabels) |
594 int heur=heur1; //starting time interval (#of relabels) |
589 int numrelabel=0; |
595 int numrelabel=0; |
590 |
596 |
591 bool what_heur=1; |
597 bool what_heur=1; |
592 //It is 0 in case 'bound decrease' and 1 in case 'highest label' |
598 //It is 0 in case 'bound decrease' and 1 in case 'highest label' |
593 |
599 |
594 bool end=false; |
600 bool end=false; |
595 //Needed for 'bound decrease', true means no active nodes are above bound b. |
601 //Needed for 'bound decrease', true means no active nodes are above bound |
|
602 //b. |
596 |
603 |
597 int k=n-2; //bound on the highest level under n containing a node |
604 int k=n-2; //bound on the highest level under n containing a node |
598 int b=k; //bound on the highest level under n of an active node |
605 int b=k; //bound on the highest level under n of an active node |
599 |
606 |
600 VecStack active(n); |
607 VecStack active(n); |
601 |
608 |
602 NNMap left(*g, INVALID); |
609 NNMap left(*g, INVALID); |
603 NNMap right(*g, INVALID); |
610 NNMap right(*g, INVALID); |
604 VecNode level_list(n,INVALID); |
611 VecNode level_list(n,INVALID); |
605 //List of the nodes in level i<n, set to n. |
612 //List of the nodes in level i<n, set to n. |
606 |
613 |
607 NodeIt v; |
614 NodeIt v; |
608 for(g->first(v); g->valid(v); g->next(v)) level.set(v,n); |
615 for(g->first(v); g->valid(v); g->next(v)) level.set(v,n); |
609 //setting each node to level n |
616 //setting each node to level n |
610 |
617 |
611 switch ( fe ) { |
618 switch (fe) { |
612 case PREFLOW: |
619 case PRE_FLOW: |
613 { |
620 { |
614 //counting the excess |
621 //computing the excess |
615 NodeIt v; |
622 NodeIt v; |
616 for(g->first(v); g->valid(v); g->next(v)) { |
623 for(g->first(v); g->valid(v); g->next(v)) { |
617 Num exc=0; |
624 Num exc=0; |
618 |
625 |
619 InEdgeIt e; |
626 InEdgeIt e; |
620 for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e]; |
627 for(g->first(e,v); g->valid(e); g->next(e)) exc+=(*flow)[e]; |
621 OutEdgeIt f; |
628 OutEdgeIt f; |
622 for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f]; |
629 for(g->first(f,v); g->valid(f); g->next(f)) exc-=(*flow)[f]; |
623 |
630 |
624 excess.set(v,exc); |
631 excess.set(v,exc); |
625 |
632 |
626 //putting the active nodes into the stack |
633 //putting the active nodes into the stack |
627 int lev=level[v]; |
634 int lev=level[v]; |
628 if ( exc > 0 && lev < n && v != t ) active[lev].push(v); |
635 if ( exc > 0 && lev < n && v != t ) active[lev].push(v); |
629 } |
636 } |
630 break; |
637 break; |
631 } |
638 } |
632 case GEN_FLOW: |
639 case GEN_FLOW: |
633 { |
640 { |
634 //Counting the excess of t |
641 //computing the excess of t |
635 Num exc=0; |
642 Num exc=0; |
636 |
643 |
637 InEdgeIt e; |
644 InEdgeIt e; |
638 for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e]; |
645 for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e]; |
639 OutEdgeIt f; |
646 OutEdgeIt f; |
640 for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f]; |
647 for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f]; |
641 |
648 |
642 excess.set(t,exc); |
649 excess.set(t,exc); |
643 |
650 |
644 break; |
651 break; |
645 } |
652 } |
646 default: |
653 default:; |
647 break; |
654 } |
648 } |
655 |
649 |
656 preflowPreproc(fe, active, level_list, left, right); |
650 preflowPreproc( fe, active, level_list, left, right ); |
657 //End of preprocessing |
651 //End of preprocessing |
658 |
652 |
659 |
653 |
|
654 //Push/relabel on the highest level active nodes. |
660 //Push/relabel on the highest level active nodes. |
655 while ( true ) { |
661 while ( true ) { |
656 if ( b == 0 ) { |
662 if ( b == 0 ) { |
657 if ( !what_heur && !end && k > 0 ) { |
663 if ( !what_heur && !end && k > 0 ) { |
658 b=k; |
664 b=k; |
659 end=true; |
665 end=true; |
660 } else break; |
666 } else break; |
661 } |
667 } |
662 |
668 |
663 if ( active[b].empty() ) --b; |
669 if ( active[b].empty() ) --b; |
664 else { |
670 else { |
665 end=false; |
671 end=false; |
666 Node w=active[b].top(); |
672 Node w=active[b].top(); |
667 active[b].pop(); |
673 active[b].pop(); |
668 int newlevel=push(w,active); |
674 int newlevel=push(w,active); |
669 if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, |
675 if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, |
670 left, right, b, k, what_heur); |
676 left, right, b, k, what_heur); |
671 |
677 |
672 ++numrelabel; |
678 ++numrelabel; |
673 if ( numrelabel >= heur ) { |
679 if ( numrelabel >= heur ) { |
674 numrelabel=0; |
680 numrelabel=0; |
675 if ( what_heur ) { |
681 if ( what_heur ) { |
676 what_heur=0; |
682 what_heur=0; |
677 heur=heur0; |
683 heur=heur0; |
678 end=false; |
684 end=false; |
679 } else { |
685 } else { |
680 what_heur=1; |
686 what_heur=1; |
681 heur=heur1; |
687 heur=heur1; |
682 b=k; |
688 b=k; |
683 } |
689 } |
684 } |
690 } |
685 } |
691 } |
686 } |
692 } |
687 } |
693 } |
688 |
694 |
689 |
695 |
690 |
696 |
691 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
697 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
692 void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1() |
698 void MaxFlow<Graph, Num, CapMap, FlowMap>::preflowPhase1() |
693 { |
699 { |
694 |
700 |
695 int k=n-2; //bound on the highest level under n containing a node |
701 int k=n-2; //bound on the highest level under n containing a node |
696 int b=k; //bound on the highest level under n of an active node |
702 int b=k; //bound on the highest level under n of an active node |
697 |
703 |
698 VecStack active(n); |
704 VecStack active(n); |
699 level.set(s,0); |
705 level.set(s,0); |
700 std::queue<Node> bfs_queue; |
706 std::queue<Node> bfs_queue; |
701 bfs_queue.push(s); |
707 bfs_queue.push(s); |
702 |
708 |
703 while (!bfs_queue.empty()) { |
709 while (!bfs_queue.empty()) { |
704 |
710 |
705 Node v=bfs_queue.front(); |
711 Node v=bfs_queue.front(); |
706 bfs_queue.pop(); |
712 bfs_queue.pop(); |
707 int l=level[v]+1; |
713 int l=level[v]+1; |
708 |
714 |
709 InEdgeIt e; |
715 InEdgeIt e; |
710 for(g->first(e,v); g->valid(e); g->next(e)) { |
716 for(g->first(e,v); g->valid(e); g->next(e)) { |
711 if ( (*capacity)[e] <= (*flow)[e] ) continue; |
717 if ( (*capacity)[e] <= (*flow)[e] ) continue; |
712 Node u=g->tail(e); |
718 Node u=g->tail(e); |
713 if ( level[u] >= n ) { |
719 if ( level[u] >= n ) { |
714 bfs_queue.push(u); |
720 bfs_queue.push(u); |
715 level.set(u, l); |
721 level.set(u, l); |
716 if ( excess[u] > 0 ) active[l].push(u); |
722 if ( excess[u] > 0 ) active[l].push(u); |
717 } |
723 } |
718 } |
724 } |
719 |
725 |
720 OutEdgeIt f; |
726 OutEdgeIt f; |
721 for(g->first(f,v); g->valid(f); g->next(f)) { |
727 for(g->first(f,v); g->valid(f); g->next(f)) { |
722 if ( 0 >= (*flow)[f] ) continue; |
728 if ( 0 >= (*flow)[f] ) continue; |
723 Node u=g->head(f); |
729 Node u=g->head(f); |
724 if ( level[u] >= n ) { |
730 if ( level[u] >= n ) { |
725 bfs_queue.push(u); |
731 bfs_queue.push(u); |
726 level.set(u, l); |
732 level.set(u, l); |
727 if ( excess[u] > 0 ) active[l].push(u); |
733 if ( excess[u] > 0 ) active[l].push(u); |
728 } |
734 } |
729 } |
735 } |
730 } |
736 } |
731 b=n-2; |
737 b=n-2; |
732 |
738 |
733 while ( true ) { |
739 while ( true ) { |
734 |
740 |
735 if ( b == 0 ) break; |
741 if ( b == 0 ) break; |
736 |
742 |
737 if ( active[b].empty() ) --b; |
743 if ( active[b].empty() ) --b; |
738 else { |
744 else { |
739 Node w=active[b].top(); |
745 Node w=active[b].top(); |
740 active[b].pop(); |
746 active[b].pop(); |
741 int newlevel=push(w,active); |
747 int newlevel=push(w,active); |
742 |
748 |
743 //relabel |
749 //relabel |
744 if ( excess[w] > 0 ) { |
750 if ( excess[w] > 0 ) { |
745 level.set(w,++newlevel); |
751 level.set(w,++newlevel); |
746 active[newlevel].push(w); |
752 active[newlevel].push(w); |
751 } |
757 } |
752 |
758 |
753 |
759 |
754 |
760 |
755 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
761 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
756 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() |
762 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath() |
757 { |
763 { |
758 ResGW res_graph(*g, *capacity, *flow); |
764 ResGW res_graph(*g, *capacity, *flow); |
759 bool _augment=false; |
765 bool _augment=false; |
760 |
766 |
761 //ReachedMap level(res_graph); |
767 //ReachedMap level(res_graph); |
762 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
768 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
763 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
769 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
764 bfs.pushAndSetReached(s); |
770 bfs.pushAndSetReached(s); |
765 |
771 |
766 typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); |
772 typename ResGW::template NodeMap<ResGWEdge> pred(res_graph); |
767 pred.set(s, INVALID); |
773 pred.set(s, INVALID); |
768 |
774 |
769 typename ResGW::template NodeMap<Num> free(res_graph); |
775 typename ResGW::template NodeMap<Num> free(res_graph); |
770 |
776 |
771 //searching for augmenting path |
777 //searching for augmenting path |
772 while ( !bfs.finished() ) { |
778 while ( !bfs.finished() ) { |
773 ResGWOutEdgeIt e=bfs; |
779 ResGWOutEdgeIt e=bfs; |
774 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { |
780 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { |
775 Node v=res_graph.tail(e); |
781 Node v=res_graph.tail(e); |
776 Node w=res_graph.head(e); |
782 Node w=res_graph.head(e); |
777 pred.set(w, e); |
783 pred.set(w, e); |
778 if (res_graph.valid(pred[v])) { |
784 if (res_graph.valid(pred[v])) { |
779 free.set(w, std::min(free[v], res_graph.resCap(e))); |
785 free.set(w, std::min(free[v], res_graph.resCap(e))); |
780 } else { |
786 } else { |
781 free.set(w, res_graph.resCap(e)); |
787 free.set(w, res_graph.resCap(e)); |
782 } |
788 } |
783 if (res_graph.head(e)==t) { _augment=true; break; } |
789 if (res_graph.head(e)==t) { _augment=true; break; } |
784 } |
790 } |
785 |
791 |
786 ++bfs; |
792 ++bfs; |
787 } //end of searching augmenting path |
793 } //end of searching augmenting path |
788 |
794 |
789 if (_augment) { |
795 if (_augment) { |
790 Node n=t; |
796 Node n=t; |
791 Num augment_value=free[t]; |
797 Num augment_value=free[t]; |
792 while (res_graph.valid(pred[n])) { |
798 while (res_graph.valid(pred[n])) { |
793 ResGWEdge e=pred[n]; |
799 ResGWEdge e=pred[n]; |
794 res_graph.augment(e, augment_value); |
800 res_graph.augment(e, augment_value); |
795 n=res_graph.tail(e); |
801 n=res_graph.tail(e); |
796 } |
802 } |
797 } |
803 } |
798 |
804 |
799 return _augment; |
805 return _augment; |
803 |
809 |
804 |
810 |
805 |
811 |
806 |
812 |
807 |
813 |
808 |
|
809 |
|
810 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
814 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
811 template<typename MutableGraph> |
815 template<typename MutableGraph> |
812 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() |
816 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow() |
813 { |
817 { |
814 typedef MutableGraph MG; |
818 typedef MutableGraph MG; |
815 bool _augment=false; |
819 bool _augment=false; |
816 |
820 |
817 ResGW res_graph(*g, *capacity, *flow); |
821 ResGW res_graph(*g, *capacity, *flow); |
818 |
822 |
819 //bfs for distances on the residual graph |
823 //bfs for distances on the residual graph |
820 //ReachedMap level(res_graph); |
824 //ReachedMap level(res_graph); |
821 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
825 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
822 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
826 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
823 bfs.pushAndSetReached(s); |
827 bfs.pushAndSetReached(s); |
824 typename ResGW::template NodeMap<int> |
828 typename ResGW::template NodeMap<int> |
825 dist(res_graph); //filled up with 0's |
829 dist(res_graph); //filled up with 0's |
826 |
830 |
827 //F will contain the physical copy of the residual graph |
831 //F will contain the physical copy of the residual graph |
828 //with the set of edges which are on shortest paths |
832 //with the set of edges which are on shortest paths |
829 MG F; |
833 MG F; |
830 typename ResGW::template NodeMap<typename MG::Node> |
834 typename ResGW::template NodeMap<typename MG::Node> |
831 res_graph_to_F(res_graph); |
835 res_graph_to_F(res_graph); |
832 { |
836 { |
833 typename ResGW::NodeIt n; |
837 typename ResGW::NodeIt n; |
834 for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) { |
838 for(res_graph.first(n); res_graph.valid(n); res_graph.next(n)) { |
835 res_graph_to_F.set(n, F.addNode()); |
839 res_graph_to_F.set(n, F.addNode()); |
839 typename MG::Node sF=res_graph_to_F[s]; |
843 typename MG::Node sF=res_graph_to_F[s]; |
840 typename MG::Node tF=res_graph_to_F[t]; |
844 typename MG::Node tF=res_graph_to_F[t]; |
841 typename MG::template EdgeMap<ResGWEdge> original_edge(F); |
845 typename MG::template EdgeMap<ResGWEdge> original_edge(F); |
842 typename MG::template EdgeMap<Num> residual_capacity(F); |
846 typename MG::template EdgeMap<Num> residual_capacity(F); |
843 |
847 |
844 while ( !bfs.finished() ) { |
848 while ( !bfs.finished() ) { |
845 ResGWOutEdgeIt e=bfs; |
849 ResGWOutEdgeIt e=bfs; |
846 if (res_graph.valid(e)) { |
850 if (res_graph.valid(e)) { |
847 if (bfs.isBNodeNewlyReached()) { |
851 if (bfs.isBNodeNewlyReached()) { |
848 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); |
852 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); |
849 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]); |
853 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], |
|
854 res_graph_to_F[res_graph.head(e)]); |
850 original_edge.update(); |
855 original_edge.update(); |
851 original_edge.set(f, e); |
856 original_edge.set(f, e); |
852 residual_capacity.update(); |
857 residual_capacity.update(); |
853 residual_capacity.set(f, res_graph.resCap(e)); |
858 residual_capacity.set(f, res_graph.resCap(e)); |
854 } else { |
859 } else { |
855 if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) { |
860 if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) { |
856 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], res_graph_to_F[res_graph.head(e)]); |
861 typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)], |
|
862 res_graph_to_F[res_graph.head(e)]); |
857 original_edge.update(); |
863 original_edge.update(); |
858 original_edge.set(f, e); |
864 original_edge.set(f, e); |
859 residual_capacity.update(); |
865 residual_capacity.update(); |
860 residual_capacity.set(f, res_graph.resCap(e)); |
866 residual_capacity.set(f, res_graph.resCap(e)); |
861 } |
867 } |
874 pred.set(sF, INVALID); |
880 pred.set(sF, INVALID); |
875 //invalid iterators for sources |
881 //invalid iterators for sources |
876 |
882 |
877 typename MG::template NodeMap<Num> free(F); |
883 typename MG::template NodeMap<Num> free(F); |
878 |
884 |
879 dfs.pushAndSetReached(sF); |
885 dfs.pushAndSetReached(sF); |
880 while (!dfs.finished()) { |
886 while (!dfs.finished()) { |
881 ++dfs; |
887 ++dfs; |
882 if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) { |
888 if (F.valid(/*typename MG::OutEdgeIt*/(dfs))) { |
883 if (dfs.isBNodeNewlyReached()) { |
889 if (dfs.isBNodeNewlyReached()) { |
884 typename MG::Node v=F.aNode(dfs); |
890 typename MG::Node v=F.aNode(dfs); |
885 typename MG::Node w=F.bNode(dfs); |
891 typename MG::Node w=F.bNode(dfs); |
886 pred.set(w, dfs); |
892 pred.set(w, dfs); |
887 if (F.valid(pred[v])) { |
893 if (F.valid(pred[v])) { |
888 free.set(w, std::min(free[v], residual_capacity[dfs])); |
894 free.set(w, std::min(free[v], residual_capacity[dfs])); |
889 } else { |
895 } else { |
890 free.set(w, residual_capacity[dfs]); |
896 free.set(w, residual_capacity[dfs]); |
891 } |
897 } |
892 if (w==tF) { |
898 if (w==tF) { |
893 __augment=true; |
899 __augment=true; |
894 _augment=true; |
900 _augment=true; |
895 break; |
901 break; |
896 } |
902 } |
897 |
903 |
898 } else { |
904 } else { |
899 F.erase(/*typename MG::OutEdgeIt*/(dfs)); |
905 F.erase(/*typename MG::OutEdgeIt*/(dfs)); |
900 } |
906 } |
901 } |
907 } |
902 } |
908 } |
903 |
909 |
904 if (__augment) { |
910 if (__augment) { |
905 typename MG::Node n=tF; |
911 typename MG::Node n=tF; |
906 Num augment_value=free[tF]; |
912 Num augment_value=free[tF]; |
907 while (F.valid(pred[n])) { |
913 while (F.valid(pred[n])) { |
908 typename MG::Edge e=pred[n]; |
914 typename MG::Edge e=pred[n]; |
909 res_graph.augment(original_edge[e], augment_value); |
915 res_graph.augment(original_edge[e], augment_value); |
910 n=F.tail(e); |
916 n=F.tail(e); |
911 if (residual_capacity[e]==augment_value) |
917 if (residual_capacity[e]==augment_value) |
912 F.erase(e); |
918 F.erase(e); |
913 else |
919 else |
914 residual_capacity.set(e, residual_capacity[e]-augment_value); |
920 residual_capacity.set(e, residual_capacity[e]-augment_value); |
915 } |
921 } |
916 } |
922 } |
917 |
923 |
918 } |
924 } |
919 |
925 |
920 return _augment; |
926 return _augment; |
921 } |
927 } |
922 |
928 |
923 |
929 |
924 |
930 |
925 |
931 |
926 |
|
927 |
|
928 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
932 template <typename Graph, typename Num, typename CapMap, typename FlowMap> |
929 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() |
933 bool MaxFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2() |
930 { |
934 { |
931 bool _augment=false; |
935 bool _augment=false; |
932 |
936 |
933 ResGW res_graph(*g, *capacity, *flow); |
937 ResGW res_graph(*g, *capacity, *flow); |
934 |
938 |
935 //ReachedMap level(res_graph); |
939 //ReachedMap level(res_graph); |
936 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
940 FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0); |
937 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
941 BfsIterator<ResGW, ReachedMap> bfs(res_graph, level); |
938 |
942 |
939 bfs.pushAndSetReached(s); |
943 bfs.pushAndSetReached(s); |
940 DistanceMap<ResGW> dist(res_graph); |
944 DistanceMap<ResGW> dist(res_graph); |
941 while ( !bfs.finished() ) { |
945 while ( !bfs.finished() ) { |
942 ResGWOutEdgeIt e=bfs; |
946 ResGWOutEdgeIt e=bfs; |
943 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { |
947 if (res_graph.valid(e) && bfs.isBNodeNewlyReached()) { |
944 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); |
948 dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1); |
945 } |
949 } |
946 ++bfs; |
950 ++bfs; |
947 } //computing distances from s in the residual graph |
951 } //computing distances from s in the residual graph |
948 |
952 |
949 //Subgraph containing the edges on some shortest paths |
953 //Subgraph containing the edges on some shortest paths |
950 ConstMap<typename ResGW::Node, bool> true_map(true); |
954 ConstMap<typename ResGW::Node, bool> true_map(true); |
951 typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, |
955 typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>, |
952 DistanceMap<ResGW> > FilterResGW; |
956 DistanceMap<ResGW> > FilterResGW; |
953 FilterResGW filter_res_graph(res_graph, true_map, dist); |
957 FilterResGW filter_res_graph(res_graph, true_map, dist); |
954 |
958 |
955 //Subgraph, which is able to delete edges which are already |
959 //Subgraph, which is able to delete edges which are already |
956 //met by the dfs |
960 //met by the dfs |
957 typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> |
961 typename FilterResGW::template NodeMap<typename FilterResGW::OutEdgeIt> |
958 first_out_edges(filter_res_graph); |
962 first_out_edges(filter_res_graph); |
959 typename FilterResGW::NodeIt v; |
963 typename FilterResGW::NodeIt v; |
960 for(filter_res_graph.first(v); filter_res_graph.valid(v); |
964 for(filter_res_graph.first(v); filter_res_graph.valid(v); |
961 filter_res_graph.next(v)) |
965 filter_res_graph.next(v)) |
962 { |
966 { |
963 typename FilterResGW::OutEdgeIt e; |
967 typename FilterResGW::OutEdgeIt e; |
964 filter_res_graph.first(e, v); |
968 filter_res_graph.first(e, v); |
965 first_out_edges.set(v, e); |
969 first_out_edges.set(v, e); |
966 } |
970 } |
972 |
976 |
973 while (__augment) { |
977 while (__augment) { |
974 |
978 |
975 __augment=false; |
979 __augment=false; |
976 //computing blocking flow with dfs |
980 //computing blocking flow with dfs |
977 DfsIterator< ErasingResGW, |
981 DfsIterator< ErasingResGW, |
978 typename ErasingResGW::template NodeMap<bool> > |
982 typename ErasingResGW::template NodeMap<bool> > |
979 dfs(erasing_res_graph); |
983 dfs(erasing_res_graph); |
980 typename ErasingResGW:: |
984 typename ErasingResGW:: |
981 template NodeMap<typename ErasingResGW::OutEdgeIt> |
985 template NodeMap<typename ErasingResGW::OutEdgeIt> |
982 pred(erasing_res_graph); |
986 pred(erasing_res_graph); |
983 pred.set(s, INVALID); |
987 pred.set(s, INVALID); |
984 //invalid iterators for sources |
988 //invalid iterators for sources |
985 |
989 |
986 typename ErasingResGW::template NodeMap<Num> |
990 typename ErasingResGW::template NodeMap<Num> |
987 free1(erasing_res_graph); |
991 free1(erasing_res_graph); |
988 |
992 |
989 dfs.pushAndSetReached( |
993 dfs.pushAndSetReached |
990 typename ErasingResGW::Node( |
994 ///\bug hugo 0.2 |
991 typename FilterResGW::Node( |
995 (typename ErasingResGW::Node |
992 typename ResGW::Node(s) |
996 (typename FilterResGW::Node |
993 ) |
997 (typename ResGW::Node(s) |
994 ) |
998 ) |
995 ); |
999 ) |
|
1000 ); |
996 while (!dfs.finished()) { |
1001 while (!dfs.finished()) { |
997 ++dfs; |
1002 ++dfs; |
998 if (erasing_res_graph.valid( |
1003 if (erasing_res_graph.valid(typename ErasingResGW::OutEdgeIt(dfs))) |
999 typename ErasingResGW::OutEdgeIt(dfs))) |
1004 { |
1000 { |
|
1001 if (dfs.isBNodeNewlyReached()) { |
1005 if (dfs.isBNodeNewlyReached()) { |
1002 |
1006 |
1003 typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs); |
1007 typename ErasingResGW::Node v=erasing_res_graph.aNode(dfs); |
1004 typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs); |
1008 typename ErasingResGW::Node w=erasing_res_graph.bNode(dfs); |
1005 |
1009 |
1006 pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs)); |
1010 pred.set(w, /*typename ErasingResGW::OutEdgeIt*/(dfs)); |
1007 if (erasing_res_graph.valid(pred[v])) { |
1011 if (erasing_res_graph.valid(pred[v])) { |
1008 free1.set(w, std::min(free1[v], res_graph.resCap( |
1012 free1.set |
1009 typename ErasingResGW::OutEdgeIt(dfs)))); |
1013 (w, std::min(free1[v], res_graph.resCap |
|
1014 (typename ErasingResGW::OutEdgeIt(dfs)))); |
1010 } else { |
1015 } else { |
1011 free1.set(w, res_graph.resCap( |
1016 free1.set |
1012 typename ErasingResGW::OutEdgeIt(dfs))); |
1017 (w, res_graph.resCap |
|
1018 (typename ErasingResGW::OutEdgeIt(dfs))); |
1013 } |
1019 } |
1014 |
1020 |
1015 if (w==t) { |
1021 if (w==t) { |
1016 __augment=true; |
1022 __augment=true; |
1017 _augment=true; |
1023 _augment=true; |
1018 break; |
1024 break; |
1019 } |
1025 } |
1020 } else { |
1026 } else { |
1021 erasing_res_graph.erase(dfs); |
1027 erasing_res_graph.erase(dfs); |
1022 } |
1028 } |
1023 } |
1029 } |
1024 } |
1030 } |
1025 |
1031 |
1026 if (__augment) { |
1032 if (__augment) { |
1027 typename ErasingResGW::Node n=typename FilterResGW::Node(typename ResGW::Node(t)); |
1033 typename ErasingResGW::Node |
|
1034 n=typename FilterResGW::Node(typename ResGW::Node(t)); |
1028 // typename ResGW::NodeMap<Num> a(res_graph); |
1035 // typename ResGW::NodeMap<Num> a(res_graph); |
1029 // typename ResGW::Node b; |
1036 // typename ResGW::Node b; |
1030 // Num j=a[b]; |
1037 // Num j=a[b]; |
1031 // typename FilterResGW::NodeMap<Num> a1(filter_res_graph); |
1038 // typename FilterResGW::NodeMap<Num> a1(filter_res_graph); |
1032 // typename FilterResGW::Node b1; |
1039 // typename FilterResGW::Node b1; |
1033 // Num j1=a1[b1]; |
1040 // Num j1=a1[b1]; |
1034 // typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph); |
1041 // typename ErasingResGW::NodeMap<Num> a2(erasing_res_graph); |
1035 // typename ErasingResGW::Node b2; |
1042 // typename ErasingResGW::Node b2; |
1036 // Num j2=a2[b2]; |
1043 // Num j2=a2[b2]; |
1037 Num augment_value=free1[n]; |
1044 Num augment_value=free1[n]; |
1038 while (erasing_res_graph.valid(pred[n])) { |
1045 while (erasing_res_graph.valid(pred[n])) { |
1039 typename ErasingResGW::OutEdgeIt e=pred[n]; |
1046 typename ErasingResGW::OutEdgeIt e=pred[n]; |
1040 res_graph.augment(e, augment_value); |
1047 res_graph.augment(e, augment_value); |
1041 n=erasing_res_graph.tail(e); |
1048 n=erasing_res_graph.tail(e); |
1042 if (res_graph.resCap(e)==0) |
1049 if (res_graph.resCap(e)==0) |
1043 erasing_res_graph.erase(e); |
1050 erasing_res_graph.erase(e); |
1044 } |
1051 } |
1045 } |
1052 } |
1046 |
1053 |
1047 } //while (__augment) |
1054 } //while (__augment) |
1048 |
1055 |
1049 return _augment; |
1056 return _augment; |
1050 } |
1057 } |
1051 |
1058 |
1052 |
1059 |
1053 |
|
1054 |
|
1055 } //namespace hugo |
1060 } //namespace hugo |
1056 |
1061 |
1057 #endif //HUGO_MAX_FLOW_H |
1062 #endif //HUGO_MAX_FLOW_H |
1058 |
1063 |
1059 |
1064 |