35 |
35 |
36 ///%Preflow algorithms class. |
36 ///%Preflow algorithms class. |
37 |
37 |
38 ///This class provides an implementation of the \e preflow \e |
38 ///This class provides an implementation of the \e preflow \e |
39 ///algorithm producing a flow of maximum value in a directed |
39 ///algorithm producing a flow of maximum value in a directed |
40 ///graph. The preflow algorithms are the fastest max flow algorithms |
40 ///graph. The preflow algorithms are the fastest known max flow algorithms |
41 ///up to now. The \e source node, the \e target node, the \e |
41 ///up to now. The \e source node, the \e target node, the \e |
42 ///capacity of the edges and the \e starting \e flow value of the |
42 ///capacity of the edges and the \e starting \e flow value of the |
43 ///edges should be passed to the algorithm through the |
43 ///edges should be passed to the algorithm through the |
44 ///constructor. It is possible to change these quantities using the |
44 ///constructor. It is possible to change these quantities using the |
45 ///functions \ref setSource, \ref setTarget, \ref setCap and \ref |
45 ///functions \ref source, \ref target, \ref setCap and \ref |
46 ///setFlow. |
46 ///setFlow. |
47 /// |
47 /// |
48 ///After running \ref lemon::Preflow::phase1() "phase1()" |
48 ///After running \ref lemon::Preflow::phase1() "phase1()" |
49 ///or \ref lemon::Preflow::run() "run()", the maximal flow |
49 ///or \ref lemon::Preflow::run() "run()", the maximal flow |
50 ///value can be obtained by calling \ref flowValue(). The minimum |
50 ///value can be obtained by calling \ref flowValue(). The minimum |
53 ///the inclusionwise minimum and maximum of the minimum value cuts, |
53 ///the inclusionwise minimum and maximum of the minimum value cuts, |
54 ///resp.) |
54 ///resp.) |
55 /// |
55 /// |
56 ///\param Graph The directed graph type the algorithm runs on. |
56 ///\param Graph The directed graph type the algorithm runs on. |
57 ///\param Num The number type of the capacities and the flow values. |
57 ///\param Num The number type of the capacities and the flow values. |
58 ///\param CapMap The capacity map type. |
58 ///\param CapacityMap The capacity map type. |
59 ///\param FlowMap The flow map type. |
59 ///\param FlowMap The flow map type. |
60 /// |
60 /// |
61 ///\author Jacint Szabo |
61 ///\author Jacint Szabo |
62 template <typename Graph, typename Num, |
62 template <typename Graph, typename Num, |
63 typename CapMap=typename Graph::template EdgeMap<Num>, |
63 typename CapacityMap=typename Graph::template EdgeMap<Num>, |
64 typename FlowMap=typename Graph::template EdgeMap<Num> > |
64 typename FlowMap=typename Graph::template EdgeMap<Num> > |
65 class Preflow { |
65 class Preflow { |
66 protected: |
66 protected: |
67 typedef typename Graph::Node Node; |
67 typedef typename Graph::Node Node; |
68 typedef typename Graph::NodeIt NodeIt; |
68 typedef typename Graph::NodeIt NodeIt; |
131 |
134 |
132 ///The constructor of the class. |
135 ///The constructor of the class. |
133 ///\param _G The directed graph the algorithm runs on. |
136 ///\param _G The directed graph the algorithm runs on. |
134 ///\param _s The source node. |
137 ///\param _s The source node. |
135 ///\param _t The target node. |
138 ///\param _t The target node. |
136 ///\param _capacity The capacity of the edges. |
139 ///\param _cap The capacity of the edges. |
137 ///\param _flow The flow of the edges. |
140 ///\param _f The flow of the edges. |
138 ///Except the graph, all of these parameters can be reset by |
141 ///Except the graph, all of these parameters can be reset by |
139 ///calling \ref setSource, \ref setTarget, \ref setCap and \ref |
142 ///calling \ref source, \ref target, \ref setCap and \ref |
140 ///setFlow, resp. |
143 ///setFlow, resp. |
141 Preflow(const Graph& _G, Node _s, Node _t, |
144 Preflow(const Graph& _gr, Node _s, Node _t, |
142 const CapMap& _capacity, FlowMap& _flow) : |
145 const CapacityMap& _cap, FlowMap& _f) : |
143 g(&_G), s(_s), t(_t), capacity(&_capacity), |
146 _g(&_gr), _source(_s), _target(_t), _capacity(&_cap), |
144 flow(&_flow), n(countNodes(_G)), level(_G), excess(_G,0), |
147 _flow(&_f), _node_num(countNodes(_gr)), level(_gr), excess(_gr,0), |
145 flow_prop(NO_FLOW), status(AFTER_NOTHING) { } |
148 flow_prop(NO_FLOW), status(AFTER_NOTHING) { } |
146 |
149 |
147 |
150 |
148 |
151 |
149 ///Runs the preflow algorithm. |
152 ///Runs the preflow algorithm. |
202 ///minCut returns a minimum cut. |
205 ///minCut returns a minimum cut. |
203 ///\warning \ref minCut(), \ref minMinCut() and \ref maxMinCut() do not |
206 ///\warning \ref minCut(), \ref minMinCut() and \ref maxMinCut() do not |
204 ///give minimum value cuts unless calling \ref phase2(). |
207 ///give minimum value cuts unless calling \ref phase2(). |
205 void phase1() |
208 void phase1() |
206 { |
209 { |
207 int heur0=(int)(H0*n); //time while running 'bound decrease' |
210 int heur0=(int)(H0*_node_num); //time while running 'bound decrease' |
208 int heur1=(int)(H1*n); //time while running 'highest label' |
211 int heur1=(int)(H1*_node_num); //time while running 'highest label' |
209 int heur=heur1; //starting time interval (#of relabels) |
212 int heur=heur1; //starting time interval (#of relabels) |
210 int numrelabel=0; |
213 int numrelabel=0; |
211 |
214 |
212 bool what_heur=1; |
215 bool what_heur=1; |
213 //It is 0 in case 'bound decrease' and 1 in case 'highest label' |
216 //It is 0 in case 'bound decrease' and 1 in case 'highest label' |
214 |
217 |
215 bool end=false; |
218 bool end=false; |
216 //Needed for 'bound decrease', true means no active |
219 //Needed for 'bound decrease', true means no active |
217 //nodes are above bound b. |
220 //nodes are above bound b. |
218 |
221 |
219 int k=n-2; //bound on the highest level under n containing a node |
222 int k=_node_num-2; //bound on the highest level under n containing a node |
220 int b=k; //bound on the highest level under n of an active node |
223 int b=k; //bound on the highest level under n of an active node |
221 |
224 |
222 VecNode first(n, INVALID); |
225 VecNode first(_node_num, INVALID); |
223 NNMap next(*g, INVALID); |
226 NNMap next(*_g, INVALID); |
224 |
227 |
225 NNMap left(*g, INVALID); |
228 NNMap left(*_g, INVALID); |
226 NNMap right(*g, INVALID); |
229 NNMap right(*_g, INVALID); |
227 VecNode level_list(n,INVALID); |
230 VecNode level_list(_node_num,INVALID); |
228 //List of the nodes in level i<n, set to n. |
231 //List of the nodes in level i<n, set to n. |
229 |
232 |
230 preflowPreproc(first, next, level_list, left, right); |
233 preflowPreproc(first, next, level_list, left, right); |
231 |
234 |
232 //Push/relabel on the highest level active nodes. |
235 //Push/relabel on the highest level active nodes. |
285 ///maxMinCut return the inclusionwise minimum and maximum cuts of |
289 ///maxMinCut return the inclusionwise minimum and maximum cuts of |
286 ///minimum value, resp. \pre \ref phase1 must be called before. |
290 ///minimum value, resp. \pre \ref phase1 must be called before. |
287 void phase2() |
291 void phase2() |
288 { |
292 { |
289 |
293 |
290 int k=n-2; //bound on the highest level under n containing a node |
294 int k=_node_num-2; //bound on the highest level under n containing a node |
291 int b=k; //bound on the highest level under n of an active node |
295 int b=k; //bound on the highest level under n of an active node |
292 |
296 |
293 |
297 |
294 VecNode first(n, INVALID); |
298 VecNode first(_node_num, INVALID); |
295 NNMap next(*g, INVALID); |
299 NNMap next(*_g, INVALID); |
296 level.set(s,0); |
300 level.set(_source,0); |
297 std::queue<Node> bfs_queue; |
301 std::queue<Node> bfs_queue; |
298 bfs_queue.push(s); |
302 bfs_queue.push(_source); |
299 |
303 |
300 while ( !bfs_queue.empty() ) { |
304 while ( !bfs_queue.empty() ) { |
301 |
305 |
302 Node v=bfs_queue.front(); |
306 Node v=bfs_queue.front(); |
303 bfs_queue.pop(); |
307 bfs_queue.pop(); |
304 int l=level[v]+1; |
308 int l=level[v]+1; |
305 |
309 |
306 for(InEdgeIt e(*g,v); e!=INVALID; ++e) { |
310 for(InEdgeIt e(*_g,v); e!=INVALID; ++e) { |
307 if ( (*capacity)[e] <= (*flow)[e] ) continue; |
311 if ( (*_capacity)[e] <= (*_flow)[e] ) continue; |
308 Node u=g->source(e); |
312 Node u=_g->source(e); |
309 if ( level[u] >= n ) { |
313 if ( level[u] >= _node_num ) { |
310 bfs_queue.push(u); |
314 bfs_queue.push(u); |
311 level.set(u, l); |
315 level.set(u, l); |
312 if ( excess[u] > 0 ) { |
316 if ( excess[u] > 0 ) { |
313 next.set(u,first[l]); |
317 next.set(u,first[l]); |
314 first[l]=u; |
318 first[l]=u; |
315 } |
319 } |
316 } |
320 } |
317 } |
321 } |
318 |
322 |
319 for(OutEdgeIt e(*g,v); e!=INVALID; ++e) { |
323 for(OutEdgeIt e(*_g,v); e!=INVALID; ++e) { |
320 if ( 0 >= (*flow)[e] ) continue; |
324 if ( 0 >= (*_flow)[e] ) continue; |
321 Node u=g->target(e); |
325 Node u=_g->target(e); |
322 if ( level[u] >= n ) { |
326 if ( level[u] >= _node_num ) { |
323 bfs_queue.push(u); |
327 bfs_queue.push(u); |
324 level.set(u, l); |
328 level.set(u, l); |
325 if ( excess[u] > 0 ) { |
329 if ( excess[u] > 0 ) { |
326 next.set(u,first[l]); |
330 next.set(u,first[l]); |
327 first[l]=u; |
331 first[l]=u; |
328 } |
332 } |
329 } |
333 } |
330 } |
334 } |
331 } |
335 } |
332 b=n-2; |
336 b=_node_num-2; |
333 |
337 |
334 while ( true ) { |
338 while ( true ) { |
335 |
339 |
336 if ( b == 0 ) break; |
340 if ( b == 0 ) break; |
337 if ( first[b]==INVALID ) --b; |
341 if ( first[b]==INVALID ) --b; |
400 ///phase2 should already be run. |
404 ///phase2 should already be run. |
401 template<typename _CutMap> |
405 template<typename _CutMap> |
402 void minMinCut(_CutMap& M) const { |
406 void minMinCut(_CutMap& M) const { |
403 |
407 |
404 std::queue<Node> queue; |
408 std::queue<Node> queue; |
405 M.set(s,true); |
409 M.set(_source,true); |
406 queue.push(s); |
410 queue.push(s); |
407 |
411 |
408 while (!queue.empty()) { |
412 while (!queue.empty()) { |
409 Node w=queue.front(); |
413 Node w=queue.front(); |
410 queue.pop(); |
414 queue.pop(); |
411 |
415 |
412 for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) { |
416 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
413 Node v=g->target(e); |
417 Node v=_g->target(e); |
414 if (!M[v] && (*flow)[e] < (*capacity)[e] ) { |
418 if (!M[v] && (*_flow)[e] < (*_capacity)[e] ) { |
415 queue.push(v); |
419 queue.push(v); |
416 M.set(v, true); |
420 M.set(v, true); |
417 } |
421 } |
418 } |
422 } |
419 |
423 |
420 for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) { |
424 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
421 Node v=g->source(e); |
425 Node v=_g->source(e); |
422 if (!M[v] && (*flow)[e] > 0 ) { |
426 if (!M[v] && (*_flow)[e] > 0 ) { |
423 queue.push(v); |
427 queue.push(v); |
424 M.set(v, true); |
428 M.set(v, true); |
425 } |
429 } |
426 } |
430 } |
427 } |
431 } |
434 ///backward bfs from the target node \c t in the residual graph. |
438 ///backward bfs from the target node \c t in the residual graph. |
435 ///\pre \ref phase2() or run() should already be run. |
439 ///\pre \ref phase2() or run() should already be run. |
436 template<typename _CutMap> |
440 template<typename _CutMap> |
437 void maxMinCut(_CutMap& M) const { |
441 void maxMinCut(_CutMap& M) const { |
438 |
442 |
439 for(NodeIt v(*g) ; v!=INVALID; ++v) M.set(v, true); |
443 for(NodeIt v(*_g) ; v!=INVALID; ++v) M.set(v, true); |
440 |
444 |
441 std::queue<Node> queue; |
445 std::queue<Node> queue; |
442 |
446 |
443 M.set(t,false); |
447 M.set(_target,false); |
444 queue.push(t); |
448 queue.push(_target); |
445 |
449 |
446 while (!queue.empty()) { |
450 while (!queue.empty()) { |
447 Node w=queue.front(); |
451 Node w=queue.front(); |
448 queue.pop(); |
452 queue.pop(); |
449 |
453 |
450 for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) { |
454 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
451 Node v=g->source(e); |
455 Node v=_g->source(e); |
452 if (M[v] && (*flow)[e] < (*capacity)[e] ) { |
456 if (M[v] && (*_flow)[e] < (*_capacity)[e] ) { |
453 queue.push(v); |
457 queue.push(v); |
454 M.set(v, false); |
458 M.set(v, false); |
455 } |
459 } |
456 } |
460 } |
457 |
461 |
458 for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) { |
462 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
459 Node v=g->target(e); |
463 Node v=_g->target(e); |
460 if (M[v] && (*flow)[e] > 0 ) { |
464 if (M[v] && (*_flow)[e] > 0 ) { |
461 queue.push(v); |
465 queue.push(v); |
462 M.set(v, false); |
466 M.set(v, false); |
463 } |
467 } |
464 } |
468 } |
465 } |
469 } |
467 |
471 |
468 ///Sets the source node to \c _s. |
472 ///Sets the source node to \c _s. |
469 |
473 |
470 ///Sets the source node to \c _s. |
474 ///Sets the source node to \c _s. |
471 /// |
475 /// |
472 void setSource(Node _s) { |
476 void source(Node _s) { |
473 s=_s; |
477 _source=_s; |
474 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW; |
478 if ( flow_prop != ZERO_FLOW ) flow_prop=NO_FLOW; |
475 status=AFTER_NOTHING; |
479 status=AFTER_NOTHING; |
476 } |
480 } |
477 |
481 |
|
482 ///Returns the source node. |
|
483 |
|
484 ///Returns the source node. |
|
485 /// |
|
486 Node source() const { |
|
487 return _source; |
|
488 } |
|
489 |
478 ///Sets the target node to \c _t. |
490 ///Sets the target node to \c _t. |
479 |
491 |
480 ///Sets the target node to \c _t. |
492 ///Sets the target node to \c _t. |
481 /// |
493 /// |
482 void setTarget(Node _t) { |
494 void target(Node _t) { |
483 t=_t; |
495 _target=_t; |
484 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW; |
496 if ( flow_prop == GEN_FLOW ) flow_prop=PRE_FLOW; |
485 status=AFTER_NOTHING; |
497 status=AFTER_NOTHING; |
486 } |
498 } |
487 |
499 |
|
500 ///Returns the target node. |
|
501 |
|
502 ///Returns the target node. |
|
503 /// |
|
504 Node target() const { |
|
505 return _target; |
|
506 } |
|
507 |
488 /// Sets the edge map of the capacities to _cap. |
508 /// Sets the edge map of the capacities to _cap. |
489 |
509 |
490 /// Sets the edge map of the capacities to _cap. |
510 /// Sets the edge map of the capacities to _cap. |
491 /// |
511 /// |
492 void setCap(const CapMap& _cap) { |
512 void capacityMap(const CapacityMap& _cap) { |
493 capacity=&_cap; |
513 _capacity=&_cap; |
494 status=AFTER_NOTHING; |
514 status=AFTER_NOTHING; |
|
515 } |
|
516 /// Returns a reference to to capacity map. |
|
517 |
|
518 /// Returns a reference to to capacity map. |
|
519 /// |
|
520 const CapacityMap &capacityMap() const { |
|
521 return *_capacity; |
495 } |
522 } |
496 |
523 |
497 /// Sets the edge map of the flows to _flow. |
524 /// Sets the edge map of the flows to _flow. |
498 |
525 |
499 /// Sets the edge map of the flows to _flow. |
526 /// Sets the edge map of the flows to _flow. |
500 /// |
527 /// |
501 void setFlow(FlowMap& _flow) { |
528 void flowMap(FlowMap& _f) { |
502 flow=&_flow; |
529 _flow=&_f; |
503 flow_prop=NO_FLOW; |
530 flow_prop=NO_FLOW; |
504 status=AFTER_NOTHING; |
531 status=AFTER_NOTHING; |
505 } |
532 } |
506 |
533 |
|
534 /// Returns a reference to to flow map. |
|
535 |
|
536 /// Returns a reference to to flow map. |
|
537 /// |
|
538 const FlowMap &flowMap() const { |
|
539 return *_flow; |
|
540 } |
507 |
541 |
508 private: |
542 private: |
509 |
543 |
510 int push(Node w, NNMap& next, VecNode& first) { |
544 int push(Node w, NNMap& next, VecNode& first) { |
511 |
545 |
512 int lev=level[w]; |
546 int lev=level[w]; |
513 Num exc=excess[w]; |
547 Num exc=excess[w]; |
514 int newlevel=n; //bound on the next level of w |
548 int newlevel=_node_num; //bound on the next level of w |
515 |
549 |
516 for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) { |
550 for(OutEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
517 if ( (*flow)[e] >= (*capacity)[e] ) continue; |
551 if ( (*_flow)[e] >= (*_capacity)[e] ) continue; |
518 Node v=g->target(e); |
552 Node v=_g->target(e); |
519 |
553 |
520 if( lev > level[v] ) { //Push is allowed now |
554 if( lev > level[v] ) { //Push is allowed now |
521 |
555 |
522 if ( excess[v]<=0 && v!=t && v!=s ) { |
556 if ( excess[v]<=0 && v!=_target && v!=_source ) { |
523 next.set(v,first[level[v]]); |
557 next.set(v,first[level[v]]); |
524 first[level[v]]=v; |
558 first[level[v]]=v; |
525 } |
559 } |
526 |
560 |
527 Num cap=(*capacity)[e]; |
561 Num cap=(*_capacity)[e]; |
528 Num flo=(*flow)[e]; |
562 Num flo=(*_flow)[e]; |
529 Num remcap=cap-flo; |
563 Num remcap=cap-flo; |
530 |
564 |
531 if ( remcap >= exc ) { //A nonsaturating push. |
565 if ( remcap >= exc ) { //A nonsaturating push. |
532 |
566 |
533 flow->set(e, flo+exc); |
567 _flow->set(e, flo+exc); |
534 excess.set(v, excess[v]+exc); |
568 excess.set(v, excess[v]+exc); |
535 exc=0; |
569 exc=0; |
536 break; |
570 break; |
537 |
571 |
538 } else { //A saturating push. |
572 } else { //A saturating push. |
539 flow->set(e, cap); |
573 _flow->set(e, cap); |
540 excess.set(v, excess[v]+remcap); |
574 excess.set(v, excess[v]+remcap); |
541 exc-=remcap; |
575 exc-=remcap; |
542 } |
576 } |
543 } else if ( newlevel > level[v] ) newlevel = level[v]; |
577 } else if ( newlevel > level[v] ) newlevel = level[v]; |
544 } //for out edges wv |
578 } //for out edges wv |
545 |
579 |
546 if ( exc > 0 ) { |
580 if ( exc > 0 ) { |
547 for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) { |
581 for(InEdgeIt e(*_g,w) ; e!=INVALID; ++e) { |
548 |
582 |
549 if( (*flow)[e] <= 0 ) continue; |
583 if( (*_flow)[e] <= 0 ) continue; |
550 Node v=g->source(e); |
584 Node v=_g->source(e); |
551 |
585 |
552 if( lev > level[v] ) { //Push is allowed now |
586 if( lev > level[v] ) { //Push is allowed now |
553 |
587 |
554 if ( excess[v]<=0 && v!=t && v!=s ) { |
588 if ( excess[v]<=0 && v!=_target && v!=_source ) { |
555 next.set(v,first[level[v]]); |
589 next.set(v,first[level[v]]); |
556 first[level[v]]=v; |
590 first[level[v]]=v; |
557 } |
591 } |
558 |
592 |
559 Num flo=(*flow)[e]; |
593 Num flo=(*_flow)[e]; |
560 |
594 |
561 if ( flo >= exc ) { //A nonsaturating push. |
595 if ( flo >= exc ) { //A nonsaturating push. |
562 |
596 |
563 flow->set(e, flo-exc); |
597 _flow->set(e, flo-exc); |
564 excess.set(v, excess[v]+exc); |
598 excess.set(v, excess[v]+exc); |
565 exc=0; |
599 exc=0; |
566 break; |
600 break; |
567 } else { //A saturating push. |
601 } else { //A saturating push. |
568 |
602 |
569 excess.set(v, excess[v]+flo); |
603 excess.set(v, excess[v]+flo); |
570 exc-=flo; |
604 exc-=flo; |
571 flow->set(e,0); |
605 _flow->set(e,0); |
572 } |
606 } |
573 } else if ( newlevel > level[v] ) newlevel = level[v]; |
607 } else if ( newlevel > level[v] ) newlevel = level[v]; |
574 } //for in edges vw |
608 } //for in edges vw |
575 |
609 |
576 } // if w still has excess after the out edge for cycle |
610 } // if w still has excess after the out edge for cycle |
583 |
617 |
584 |
618 |
585 void preflowPreproc(VecNode& first, NNMap& next, |
619 void preflowPreproc(VecNode& first, NNMap& next, |
586 VecNode& level_list, NNMap& left, NNMap& right) |
620 VecNode& level_list, NNMap& left, NNMap& right) |
587 { |
621 { |
588 for(NodeIt v(*g); v!=INVALID; ++v) level.set(v,n); |
622 for(NodeIt v(*_g); v!=INVALID; ++v) level.set(v,_node_num); |
589 std::queue<Node> bfs_queue; |
623 std::queue<Node> bfs_queue; |
590 |
624 |
591 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) { |
625 if ( flow_prop == GEN_FLOW || flow_prop == PRE_FLOW ) { |
592 //Reverse_bfs from t in the residual graph, |
626 //Reverse_bfs from t in the residual graph, |
593 //to find the starting level. |
627 //to find the starting level. |
594 level.set(t,0); |
628 level.set(_target,0); |
595 bfs_queue.push(t); |
629 bfs_queue.push(_target); |
596 |
630 |
597 while ( !bfs_queue.empty() ) { |
631 while ( !bfs_queue.empty() ) { |
598 |
632 |
599 Node v=bfs_queue.front(); |
633 Node v=bfs_queue.front(); |
600 bfs_queue.pop(); |
634 bfs_queue.pop(); |
601 int l=level[v]+1; |
635 int l=level[v]+1; |
602 |
636 |
603 for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) { |
637 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
604 if ( (*capacity)[e] <= (*flow)[e] ) continue; |
638 if ( (*_capacity)[e] <= (*_flow)[e] ) continue; |
605 Node w=g->source(e); |
639 Node w=_g->source(e); |
606 if ( level[w] == n && w != s ) { |
640 if ( level[w] == _node_num && w != _source ) { |
607 bfs_queue.push(w); |
641 bfs_queue.push(w); |
608 Node z=level_list[l]; |
642 Node z=level_list[l]; |
609 if ( z!=INVALID ) left.set(z,w); |
643 if ( z!=INVALID ) left.set(z,w); |
610 right.set(w,z); |
644 right.set(w,z); |
611 level_list[l]=w; |
645 level_list[l]=w; |
612 level.set(w, l); |
646 level.set(w, l); |
613 } |
647 } |
614 } |
648 } |
615 |
649 |
616 for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) { |
650 for(OutEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
617 if ( 0 >= (*flow)[e] ) continue; |
651 if ( 0 >= (*_flow)[e] ) continue; |
618 Node w=g->target(e); |
652 Node w=_g->target(e); |
619 if ( level[w] == n && w != s ) { |
653 if ( level[w] == _node_num && w != _source ) { |
620 bfs_queue.push(w); |
654 bfs_queue.push(w); |
621 Node z=level_list[l]; |
655 Node z=level_list[l]; |
622 if ( z!=INVALID ) left.set(z,w); |
656 if ( z!=INVALID ) left.set(z,w); |
623 right.set(w,z); |
657 right.set(w,z); |
624 level_list[l]=w; |
658 level_list[l]=w; |
629 } //if |
663 } //if |
630 |
664 |
631 |
665 |
632 switch (flow_prop) { |
666 switch (flow_prop) { |
633 case NO_FLOW: |
667 case NO_FLOW: |
634 for(EdgeIt e(*g); e!=INVALID; ++e) flow->set(e,0); |
668 for(EdgeIt e(*_g); e!=INVALID; ++e) _flow->set(e,0); |
635 case ZERO_FLOW: |
669 case ZERO_FLOW: |
636 for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0); |
670 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0); |
637 |
671 |
638 //Reverse_bfs from t, to find the starting level. |
672 //Reverse_bfs from t, to find the starting level. |
639 level.set(t,0); |
673 level.set(_target,0); |
640 bfs_queue.push(t); |
674 bfs_queue.push(_target); |
641 |
675 |
642 while ( !bfs_queue.empty() ) { |
676 while ( !bfs_queue.empty() ) { |
643 |
677 |
644 Node v=bfs_queue.front(); |
678 Node v=bfs_queue.front(); |
645 bfs_queue.pop(); |
679 bfs_queue.pop(); |
646 int l=level[v]+1; |
680 int l=level[v]+1; |
647 |
681 |
648 for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) { |
682 for(InEdgeIt e(*_g,v) ; e!=INVALID; ++e) { |
649 Node w=g->source(e); |
683 Node w=_g->source(e); |
650 if ( level[w] == n && w != s ) { |
684 if ( level[w] == _node_num && w != _source ) { |
651 bfs_queue.push(w); |
685 bfs_queue.push(w); |
652 Node z=level_list[l]; |
686 Node z=level_list[l]; |
653 if ( z!=INVALID ) left.set(z,w); |
687 if ( z!=INVALID ) left.set(z,w); |
654 right.set(w,z); |
688 right.set(w,z); |
655 level_list[l]=w; |
689 level_list[l]=w; |
657 } |
691 } |
658 } |
692 } |
659 } |
693 } |
660 |
694 |
661 //the starting flow |
695 //the starting flow |
662 for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) { |
696 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
663 Num c=(*capacity)[e]; |
697 Num c=(*_capacity)[e]; |
664 if ( c <= 0 ) continue; |
698 if ( c <= 0 ) continue; |
665 Node w=g->target(e); |
699 Node w=_g->target(e); |
666 if ( level[w] < n ) { |
700 if ( level[w] < _node_num ) { |
667 if ( excess[w] <= 0 && w!=t ) { //putting into the stack |
701 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack |
668 next.set(w,first[level[w]]); |
702 next.set(w,first[level[w]]); |
669 first[level[w]]=w; |
703 first[level[w]]=w; |
670 } |
704 } |
671 flow->set(e, c); |
705 _flow->set(e, c); |
672 excess.set(w, excess[w]+c); |
706 excess.set(w, excess[w]+c); |
673 } |
707 } |
674 } |
708 } |
675 break; |
709 break; |
676 |
710 |
677 case GEN_FLOW: |
711 case GEN_FLOW: |
678 for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0); |
712 for(NodeIt v(*_g); v!=INVALID; ++v) excess.set(v,0); |
679 { |
713 { |
680 Num exc=0; |
714 Num exc=0; |
681 for(InEdgeIt e(*g,t) ; e!=INVALID; ++e) exc+=(*flow)[e]; |
715 for(InEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc+=(*_flow)[e]; |
682 for(OutEdgeIt e(*g,t) ; e!=INVALID; ++e) exc-=(*flow)[e]; |
716 for(OutEdgeIt e(*_g,_target) ; e!=INVALID; ++e) exc-=(*_flow)[e]; |
683 excess.set(t,exc); |
717 excess.set(_target,exc); |
684 } |
718 } |
685 |
719 |
686 //the starting flow |
720 //the starting flow |
687 for(OutEdgeIt e(*g,s); e!=INVALID; ++e) { |
721 for(OutEdgeIt e(*_g,_source); e!=INVALID; ++e) { |
688 Num rem=(*capacity)[e]-(*flow)[e]; |
722 Num rem=(*_capacity)[e]-(*_flow)[e]; |
689 if ( rem <= 0 ) continue; |
723 if ( rem <= 0 ) continue; |
690 Node w=g->target(e); |
724 Node w=_g->target(e); |
691 if ( level[w] < n ) { |
725 if ( level[w] < _node_num ) { |
692 if ( excess[w] <= 0 && w!=t ) { //putting into the stack |
726 if ( excess[w] <= 0 && w!=_target ) { //putting into the stack |
693 next.set(w,first[level[w]]); |
727 next.set(w,first[level[w]]); |
694 first[level[w]]=w; |
728 first[level[w]]=w; |
695 } |
729 } |
696 flow->set(e, (*capacity)[e]); |
730 _flow->set(e, (*_capacity)[e]); |
697 excess.set(w, excess[w]+rem); |
731 excess.set(w, excess[w]+rem); |
698 } |
732 } |
699 } |
733 } |
700 |
734 |
701 for(InEdgeIt e(*g,s); e!=INVALID; ++e) { |
735 for(InEdgeIt e(*_g,_source); e!=INVALID; ++e) { |
702 if ( (*flow)[e] <= 0 ) continue; |
736 if ( (*_flow)[e] <= 0 ) continue; |
703 Node w=g->source(e); |
737 Node w=_g->source(e); |
704 if ( level[w] < n ) { |
738 if ( level[w] < _node_num ) { |
705 if ( excess[w] <= 0 && w!=t ) { |
739 if ( excess[w] <= 0 && w!=_target ) { |
706 next.set(w,first[level[w]]); |
740 next.set(w,first[level[w]]); |
707 first[level[w]]=w; |
741 first[level[w]]=w; |
708 } |
742 } |
709 excess.set(w, excess[w]+(*flow)[e]); |
743 excess.set(w, excess[w]+(*_flow)[e]); |
710 flow->set(e, 0); |
744 _flow->set(e, 0); |
711 } |
745 } |
712 } |
746 } |
713 break; |
747 break; |
714 |
748 |
715 case PRE_FLOW: |
749 case PRE_FLOW: |
716 //the starting flow |
750 //the starting flow |
717 for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) { |
751 for(OutEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
718 Num rem=(*capacity)[e]-(*flow)[e]; |
752 Num rem=(*_capacity)[e]-(*_flow)[e]; |
719 if ( rem <= 0 ) continue; |
753 if ( rem <= 0 ) continue; |
720 Node w=g->target(e); |
754 Node w=_g->target(e); |
721 if ( level[w] < n ) flow->set(e, (*capacity)[e]); |
755 if ( level[w] < _node_num ) _flow->set(e, (*_capacity)[e]); |
722 } |
756 } |
723 |
757 |
724 for(InEdgeIt e(*g,s) ; e!=INVALID; ++e) { |
758 for(InEdgeIt e(*_g,_source) ; e!=INVALID; ++e) { |
725 if ( (*flow)[e] <= 0 ) continue; |
759 if ( (*_flow)[e] <= 0 ) continue; |
726 Node w=g->source(e); |
760 Node w=_g->source(e); |
727 if ( level[w] < n ) flow->set(e, 0); |
761 if ( level[w] < _node_num ) _flow->set(e, 0); |
728 } |
762 } |
729 |
763 |
730 //computing the excess |
764 //computing the excess |
731 for(NodeIt w(*g); w!=INVALID; ++w) { |
765 for(NodeIt w(*_g); w!=INVALID; ++w) { |
732 Num exc=0; |
766 Num exc=0; |
733 for(InEdgeIt e(*g,w); e!=INVALID; ++e) exc+=(*flow)[e]; |
767 for(InEdgeIt e(*_g,w); e!=INVALID; ++e) exc+=(*_flow)[e]; |
734 for(OutEdgeIt e(*g,w); e!=INVALID; ++e) exc-=(*flow)[e]; |
768 for(OutEdgeIt e(*_g,w); e!=INVALID; ++e) exc-=(*_flow)[e]; |
735 excess.set(w,exc); |
769 excess.set(w,exc); |
736 |
770 |
737 //putting the active nodes into the stack |
771 //putting the active nodes into the stack |
738 int lev=level[w]; |
772 int lev=level[w]; |
739 if ( exc > 0 && lev < n && Node(w) != t ) { |
773 if ( exc > 0 && lev < _node_num && Node(w) != _target ) { |
740 next.set(w,first[lev]); |
774 next.set(w,first[lev]); |
741 first[lev]=w; |
775 first[lev]=w; |
742 } |
776 } |
743 } |
777 } |
744 break; |
778 break; |