src/hugo/max_flow.h
changeset 835 eb9587f09b42
parent 773 ce9438c5a82d
equal deleted inserted replaced
7:9f1a4a2a8c63 8:9a220d72378d
     3 #define HUGO_MAX_FLOW_H
     3 #define HUGO_MAX_FLOW_H
     4 
     4 
     5 #include <vector>
     5 #include <vector>
     6 #include <queue>
     6 #include <queue>
     7 
     7 
     8 #include <hugo/graph_wrapper.h>
     8 //#include <hugo/graph_wrapper.h>
     9 #include <hugo/invalid.h>
     9 #include <hugo/invalid.h>
    10 #include <hugo/maps.h>
    10 #include <hugo/maps.h>
    11 
    11 
    12 /// \file
    12 /// \file
    13 /// \ingroup flowalgs
    13 /// \ingroup flowalgs
    60     Node s;
    60     Node s;
    61     Node t;
    61     Node t;
    62     const CapMap* capacity;
    62     const CapMap* capacity;
    63     FlowMap* flow;
    63     FlowMap* flow;
    64     int n;      //the number of nodes of G
    64     int n;      //the number of nodes of G
    65     typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
    65     //    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
    66     //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    66     //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
    67     typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    67     //    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
    68     typedef typename ResGW::Edge ResGWEdge;
    68     //    typedef typename ResGW::Edge ResGWEdge;
    69     typedef typename Graph::template NodeMap<int> ReachedMap;
    69     typedef typename Graph::template NodeMap<int> ReachedMap;
    70 
    70 
    71 
    71 
    72     //level works as a bool map in augmenting path algorithms and is
    72     //level works as a bool map in augmenting path algorithms and is
    73     //used by bfs for storing reached information.  In preflow, it
    73     //used by bfs for storing reached information.  In preflow, it
   110     };
   110     };
   111 
   111 
   112     /// Do not needle this flag only if necessary.
   112     /// Do not needle this flag only if necessary.
   113     StatusEnum status;
   113     StatusEnum status;
   114 
   114 
   115 //     int number_of_augmentations;
   115     //     int number_of_augmentations;
   116 
   116 
   117 
   117 
   118 //     template<typename IntMap>
   118     //     template<typename IntMap>
   119 //     class TrickyReachedMap {
   119     //     class TrickyReachedMap {
   120 //     protected:
   120     //     protected:
   121 //       IntMap* map;
   121     //       IntMap* map;
   122 //       int* number_of_augmentations;
   122     //       int* number_of_augmentations;
   123 //     public:
   123     //     public:
   124 //       TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) : 
   124     //       TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) : 
   125 // 	map(&_map), number_of_augmentations(&_number_of_augmentations) { }
   125     // 	map(&_map), number_of_augmentations(&_number_of_augmentations) { }
   126 //       void set(const Node& n, bool b) {
   126     //       void set(const Node& n, bool b) {
   127 // 	if (b)
   127     // 	if (b)
   128 // 	  map->set(n, *number_of_augmentations);
   128     // 	  map->set(n, *number_of_augmentations);
   129 // 	else 
   129     // 	else 
   130 // 	  map->set(n, *number_of_augmentations-1);
   130     // 	  map->set(n, *number_of_augmentations-1);
   131 //       }
   131     //       }
   132 //       bool operator[](const Node& n) const { 
   132     //       bool operator[](const Node& n) const { 
   133 // 	return (*map)[n]==*number_of_augmentations; 
   133     // 	return (*map)[n]==*number_of_augmentations; 
   134 //       }
   134     //       }
   135 //     };
   135     //     };
   136     
   136     
   137     ///Constructor
   137     ///Constructor
   138 
   138 
   139     ///\todo Document, please.
   139     ///\todo Document, please.
   140     ///
   140     ///
   232 	    b=k;
   232 	    b=k;
   233 	    end=true;
   233 	    end=true;
   234 	  } else break;
   234 	  } else break;
   235 	}
   235 	}
   236 
   236 
   237 	if ( !g->valid(first[b]) ) --b;
   237 	if ( first[b]==INVALID ) --b;
   238 	else {
   238 	else {
   239 	  end=false;
   239 	  end=false;
   240 	  Node w=first[b];
   240 	  Node w=first[b];
   241 	  first[b]=next[w];
   241 	  first[b]=next[w];
   242 	  int newlevel=push(w, next, first);
   242 	  int newlevel=push(w, next, first);
   287 
   287 
   288 	Node v=bfs_queue.front();
   288 	Node v=bfs_queue.front();
   289 	bfs_queue.pop();
   289 	bfs_queue.pop();
   290 	int l=level[v]+1;
   290 	int l=level[v]+1;
   291 
   291 
   292 	InEdgeIt e;
   292 	for(InEdgeIt e(*g,v); e!=INVALID; ++e) {
   293 	for(g->first(e,v); g->valid(e); g->next(e)) {
       
   294 	  if ( (*capacity)[e] <= (*flow)[e] ) continue;
   293 	  if ( (*capacity)[e] <= (*flow)[e] ) continue;
   295 	  Node u=g->tail(e);
   294 	  Node u=g->tail(e);
   296 	  if ( level[u] >= n ) {
   295 	  if ( level[u] >= n ) {
   297 	    bfs_queue.push(u);
   296 	    bfs_queue.push(u);
   298 	    level.set(u, l);
   297 	    level.set(u, l);
   301 	      first[l]=u;
   300 	      first[l]=u;
   302 	    }
   301 	    }
   303 	  }
   302 	  }
   304 	}
   303 	}
   305 
   304 
   306 	OutEdgeIt f;
   305 	for(OutEdgeIt e(*g,v); e!=INVALID; ++e) {
   307 	for(g->first(f,v); g->valid(f); g->next(f)) {
   306 	  if ( 0 >= (*flow)[e] ) continue;
   308 	  if ( 0 >= (*flow)[f] ) continue;
   307 	  Node u=g->head(e);
   309 	  Node u=g->head(f);
       
   310 	  if ( level[u] >= n ) {
   308 	  if ( level[u] >= n ) {
   311 	    bfs_queue.push(u);
   309 	    bfs_queue.push(u);
   312 	    level.set(u, l);
   310 	    level.set(u, l);
   313 	    if ( excess[u] > 0 ) {
   311 	    if ( excess[u] > 0 ) {
   314 	      next.set(u,first[l]);
   312 	      next.set(u,first[l]);
   321 
   319 
   322       while ( true ) {
   320       while ( true ) {
   323 
   321 
   324 	if ( b == 0 ) break;
   322 	if ( b == 0 ) break;
   325 
   323 
   326 	if ( !g->valid(first[b]) ) --b;
   324 	if ( first[b]==INVALID ) --b;
   327 	else {
   325 	else {
   328 
   326 
   329 	  Node w=first[b];
   327 	  Node w=first[b];
   330 	  first[b]=next[w];
   328 	  first[b]=next[w];
   331 	  int newlevel=push(w,next, first/*active*/);
   329 	  int newlevel=push(w,next, first/*active*/);
   349     /// Returns the excess of the target node \ref t. 
   347     /// Returns the excess of the target node \ref t. 
   350     /// After running \ref preflowPhase1, this is the value of 
   348     /// After running \ref preflowPhase1, this is the value of 
   351     /// the maximum flow.
   349     /// the maximum flow.
   352     /// It can be called already after running \ref preflowPhase1.
   350     /// It can be called already after running \ref preflowPhase1.
   353     Num flowValue() const {
   351     Num flowValue() const {
   354 //       Num a=0;
   352       //       Num a=0;
   355 //       for(InEdgeIt e(*g,t);g->valid(e);g->next(e)) a+=(*flow)[e];
   353       //       for(InEdgeIt e(*g,t);g->valid(e);g->next(e)) a+=(*flow)[e];
   356 //       for(OutEdgeIt e(*g,t);g->valid(e);g->next(e)) a-=(*flow)[e];
   354       //       for(OutEdgeIt e(*g,t);g->valid(e);g->next(e)) a-=(*flow)[e];
   357 //       return a;
   355       //       return a;
   358       return excess[t];
   356       return excess[t];
   359       //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
   357       //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
   360     }
   358     }
   361 
   359 
   362 
   360 
   372     /// actual state
   370     /// actual state
   373     /// of the class. This enables us to determine which methods are valid
   371     /// of the class. This enables us to determine which methods are valid
   374     /// for MinCut computation
   372     /// for MinCut computation
   375     template<typename _CutMap>
   373     template<typename _CutMap>
   376     void actMinCut(_CutMap& M) const {
   374     void actMinCut(_CutMap& M) const {
   377       NodeIt v;
       
   378       switch (status) {
   375       switch (status) {
   379       case AFTER_PRE_FLOW_PHASE_1:
   376 	case AFTER_PRE_FLOW_PHASE_1:
   380 	for(g->first(v); g->valid(v); g->next(v)) {
   377 	for(NodeIt v(*g); v!=INVALID; ++v) {
   381 	  if (level[v] < n) {
   378 	  if (level[v] < n) {
   382 	    M.set(v, false);
   379 	    M.set(v, false);
   383 	  } else {
   380 	  } else {
   384 	    M.set(v, true);
   381 	    M.set(v, true);
   385 	  }
   382 	  }
   386 	}
   383 	}
   387 	break;
   384 	break;
   388       case AFTER_PRE_FLOW_PHASE_2:
   385 	case AFTER_PRE_FLOW_PHASE_2:
   389       case AFTER_NOTHING:
   386 	case AFTER_NOTHING:
   390       case AFTER_AUGMENTING:
   387 	case AFTER_AUGMENTING:
   391       case AFTER_FAST_AUGMENTING:
   388 	case AFTER_FAST_AUGMENTING:
   392 	minMinCut(M);
   389 	minMinCut(M);
   393 	break;
   390 	break;
   394       }
   391       }
   395     }
   392     }
   396 
   393 
   410 
   407 
   411       while (!queue.empty()) {
   408       while (!queue.empty()) {
   412         Node w=queue.front();
   409         Node w=queue.front();
   413 	queue.pop();
   410 	queue.pop();
   414 
   411 
   415 	OutEdgeIt e;
   412 	for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   416 	for(g->first(e,w) ; g->valid(e); g->next(e)) {
       
   417 	  Node v=g->head(e);
   413 	  Node v=g->head(e);
   418 	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   414 	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   419 	    queue.push(v);
   415 	    queue.push(v);
   420 	    M.set(v, true);
   416 	    M.set(v, true);
   421 	  }
   417 	  }
   422 	}
   418 	}
   423 
   419 
   424 	InEdgeIt f;
   420 	for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   425 	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   421 	  Node v=g->tail(e);
   426 	  Node v=g->tail(f);
   422 	  if (!M[v] && (*flow)[e] > 0 ) {
   427 	  if (!M[v] && (*flow)[f] > 0 ) {
       
   428 	    queue.push(v);
   423 	    queue.push(v);
   429 	    M.set(v, true);
   424 	    M.set(v, true);
   430 	  }
   425 	  }
   431 	}
   426 	}
   432       }
   427       }
   440     ///\pre M should be a node map of bools initialized to false.
   435     ///\pre M should be a node map of bools initialized to false.
   441     ///\pre \c flow must be a maximum flow. 
   436     ///\pre \c flow must be a maximum flow. 
   442     template<typename _CutMap>
   437     template<typename _CutMap>
   443     void maxMinCut(_CutMap& M) const {
   438     void maxMinCut(_CutMap& M) const {
   444 
   439 
   445       NodeIt v;
   440       for(NodeIt v(*g) ; v!=INVALID; ++v) M.set(v, true);
   446       for(g->first(v) ; g->valid(v); g->next(v)) {
       
   447 	M.set(v, true);
       
   448       }
       
   449 
   441 
   450       std::queue<Node> queue;
   442       std::queue<Node> queue;
   451 
   443 
   452       M.set(t,false);
   444       M.set(t,false);
   453       queue.push(t);
   445       queue.push(t);
   454 
   446 
   455       while (!queue.empty()) {
   447       while (!queue.empty()) {
   456         Node w=queue.front();
   448         Node w=queue.front();
   457 	queue.pop();
   449 	queue.pop();
   458 
   450 
   459 	InEdgeIt e;
   451 	for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   460 	for(g->first(e,w) ; g->valid(e); g->next(e)) {
       
   461 	  Node v=g->tail(e);
   452 	  Node v=g->tail(e);
   462 	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   453 	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   463 	    queue.push(v);
   454 	    queue.push(v);
   464 	    M.set(v, false);
   455 	    M.set(v, false);
   465 	  }
   456 	  }
   466 	}
   457 	}
   467 
   458 
   468 	OutEdgeIt f;
   459 	for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   469 	for(g->first(f,w) ; g->valid(f); g->next(f)) {
   460 	  Node v=g->head(e);
   470 	  Node v=g->head(f);
   461 	  if (M[v] && (*flow)[e] > 0 ) {
   471 	  if (M[v] && (*flow)[f] > 0 ) {
       
   472 	    queue.push(v);
   462 	    queue.push(v);
   473 	    M.set(v, false);
   463 	    M.set(v, false);
   474 	  }
   464 	  }
   475 	}
   465 	}
   476       }
   466       }
   516 
   506 
   517       int lev=level[w];
   507       int lev=level[w];
   518       Num exc=excess[w];
   508       Num exc=excess[w];
   519       int newlevel=n;       //bound on the next level of w
   509       int newlevel=n;       //bound on the next level of w
   520 
   510 
   521       OutEdgeIt e;
   511       for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   522       for(g->first(e,w); g->valid(e); g->next(e)) {
       
   523 
       
   524 	if ( (*flow)[e] >= (*capacity)[e] ) continue;
   512 	if ( (*flow)[e] >= (*capacity)[e] ) continue;
   525 	Node v=g->head(e);
   513 	Node v=g->head(e);
   526 
   514 
   527 	if( lev > level[v] ) { //Push is allowed now
   515 	if( lev > level[v] ) { //Push is allowed now
   528 
   516 	  
   529 	  if ( excess[v]<=0 && v!=t && v!=s ) {
   517 	  if ( excess[v]<=0 && v!=t && v!=s ) {
   530 	    next.set(v,first[level[v]]);
   518 	    next.set(v,first[level[v]]);
   531 	    first[level[v]]=v;
   519 	    first[level[v]]=v;
   532 	  }
   520 	  }
   533 
   521 
   534 	  Num cap=(*capacity)[e];
   522 	  Num cap=(*capacity)[e];
   535 	  Num flo=(*flow)[e];
   523 	  Num flo=(*flow)[e];
   536 	  Num remcap=cap-flo;
   524 	  Num remcap=cap-flo;
   537 
   525 	  
   538 	  if ( remcap >= exc ) { //A nonsaturating push.
   526 	  if ( remcap >= exc ) { //A nonsaturating push.
   539 
   527 	    
   540 	    flow->set(e, flo+exc);
   528 	    flow->set(e, flo+exc);
   541 	    excess.set(v, excess[v]+exc);
   529 	    excess.set(v, excess[v]+exc);
   542 	    exc=0;
   530 	    exc=0;
   543 	    break;
   531 	    break;
   544 
   532 
   549 	  }
   537 	  }
   550 	} else if ( newlevel > level[v] ) newlevel = level[v];
   538 	} else if ( newlevel > level[v] ) newlevel = level[v];
   551       } //for out edges wv
   539       } //for out edges wv
   552 
   540 
   553       if ( exc > 0 ) {
   541       if ( exc > 0 ) {
   554 	InEdgeIt e;
   542 	for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) {
   555 	for(g->first(e,w); g->valid(e); g->next(e)) {
   543 	  
   556 
       
   557 	  if( (*flow)[e] <= 0 ) continue;
   544 	  if( (*flow)[e] <= 0 ) continue;
   558 	  Node v=g->tail(e);
   545 	  Node v=g->tail(e);
   559 
   546 
   560 	  if( lev > level[v] ) { //Push is allowed now
   547 	  if( lev > level[v] ) { //Push is allowed now
   561 
   548 
   582 	} //for in edges vw
   569 	} //for in edges vw
   583 
   570 
   584       } // if w still has excess after the out edge for cycle
   571       } // if w still has excess after the out edge for cycle
   585 
   572 
   586       excess.set(w, exc);
   573       excess.set(w, exc);
   587 
   574       
   588       return newlevel;
   575       return newlevel;
   589     }
   576     }
   590 
   577     
   591 
   578     
   592 
   579     
   593     void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
   580     void preflowPreproc(FlowEnum fe, NNMap& next, VecFirst& first,
   594 			VecNode& level_list, NNMap& left, NNMap& right)
   581 			VecNode& level_list, NNMap& left, NNMap& right)
   595     {
   582     {
   596       switch (fe) { //setting excess
   583       switch (fe) {  //setting excess
   597 	case NO_FLOW: 
   584 	case NO_FLOW: 
       
   585 	for(EdgeIt e(*g); e!=INVALID; ++e) flow->set(e,0);
       
   586 	for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
       
   587 	break;
       
   588 	case ZERO_FLOW: 
       
   589 	for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
       
   590 	break;
       
   591 	case GEN_FLOW:
       
   592 	for(NodeIt v(*g); v!=INVALID; ++v) excess.set(v,0);
   598 	{
   593 	{
   599 	  EdgeIt e;
       
   600 	  for(g->first(e); g->valid(e); g->next(e)) flow->set(e,0);
       
   601 	  
       
   602 	  NodeIt v;
       
   603 	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
       
   604 	  break;
       
   605 	}
       
   606 	case ZERO_FLOW: 
       
   607 	{
       
   608 	  NodeIt v;
       
   609 	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
       
   610 	  break;
       
   611 	}
       
   612 	case GEN_FLOW:
       
   613 	{
       
   614 	  NodeIt v;
       
   615 	  for(g->first(v); g->valid(v); g->next(v)) excess.set(v,0);
       
   616 
       
   617 	  Num exc=0;
   594 	  Num exc=0;
   618 	  InEdgeIt e;
   595 	  for(InEdgeIt e(*g,t) ; e!=INVALID; ++e) exc+=(*flow)[e];
   619 	  for(g->first(e,t); g->valid(e); g->next(e)) exc+=(*flow)[e];
   596 	  for(OutEdgeIt e(*g,t) ; e!=INVALID; ++e) exc-=(*flow)[e];
   620 	  OutEdgeIt f;
       
   621 	  for(g->first(f,t); g->valid(f); g->next(f)) exc-=(*flow)[f];
       
   622 	  excess.set(t,exc);
   597 	  excess.set(t,exc);
   623 	  break;
   598 	}
   624 	}
   599 	break;
   625 	default: break;
   600 	default:
   626       }
   601 	break;
   627 
   602       }
   628       NodeIt v;
   603       
   629       for(g->first(v); g->valid(v); g->next(v)) level.set(v,n);
   604       for(NodeIt v(*g); v!=INVALID; ++v) level.set(v,n);
   630       //setting each node to level n
   605       //setting each node to level n
   631       
   606       
   632       std::queue<Node> bfs_queue;
   607       std::queue<Node> bfs_queue;
   633 
   608 
   634 
   609 
   635       switch (fe) {
   610       switch (fe) {
   636       case NO_FLOW:   //flow is already set to const zero
   611       case NO_FLOW:   //flow is already set to const zero
   637       case ZERO_FLOW:
   612       case ZERO_FLOW:
   638 	{
   613 	//Reverse_bfs from t, to find the starting level.
   639 	  //Reverse_bfs from t, to find the starting level.
   614 	level.set(t,0);
   640 	  level.set(t,0);
   615 	bfs_queue.push(t);
   641 	  bfs_queue.push(t);
   616 	
   642 
   617 	while (!bfs_queue.empty()) {
   643 	  while (!bfs_queue.empty()) {
   618 	  
   644 
   619 	  Node v=bfs_queue.front();
   645 	    Node v=bfs_queue.front();
   620 	  bfs_queue.pop();
   646 	    bfs_queue.pop();
   621 	  int l=level[v]+1;
   647 	    int l=level[v]+1;
   622 	  
   648 
   623 	  for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
   649 	    InEdgeIt e;
   624 	    Node w=g->tail(e);
   650 	    for(g->first(e,v); g->valid(e); g->next(e)) {
   625 	    if ( level[w] == n && w != s ) {
   651 	      Node w=g->tail(e);
   626 	      bfs_queue.push(w);
   652 	      if ( level[w] == n && w != s ) {
   627 	      Node z=level_list[l];
   653 		bfs_queue.push(w);
   628 	      if ( z!=INVALID ) left.set(z,w);
   654 		Node z=level_list[l];
   629 	      right.set(w,z);
   655 		if ( g->valid(z) ) left.set(z,w);
   630 	      level_list[l]=w;
   656 		right.set(w,z);
   631 	      level.set(w, l);
   657 		level_list[l]=w;
   632 	    }
   658 		level.set(w, l);
   633 	  }
   659 	      }
   634 	}
   660 	    }
   635 	
   661 	  }
   636 	//the starting flow
   662 
   637 	for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e)
   663 	  //the starting flow
   638 	  {
   664 	  OutEdgeIt e;
   639 	    Num c=(*capacity)[e];
   665 	  for(g->first(e,s); g->valid(e); g->next(e))
   640 	    if ( c <= 0 ) continue;
       
   641 	    Node w=g->head(e);
       
   642 	    if ( level[w] < n ) {
       
   643 	      if ( excess[w] <= 0 && w!=t ) //putting into the stack
       
   644 		{ 
       
   645 		  next.set(w,first[level[w]]);
       
   646 		  first[level[w]]=w;
       
   647 		}
       
   648 	      flow->set(e, c);
       
   649 	      excess.set(w, excess[w]+c);
       
   650 	    }
       
   651 	  }
       
   652 	break;
       
   653       case GEN_FLOW:
       
   654 	//Reverse_bfs from t in the residual graph,
       
   655 	//to find the starting level.
       
   656 	level.set(t,0);
       
   657 	bfs_queue.push(t);
       
   658 	
       
   659 	while (!bfs_queue.empty()) {
       
   660 	  
       
   661 	  Node v=bfs_queue.front();
       
   662 	  bfs_queue.pop();
       
   663 	  int l=level[v]+1;
       
   664 	  
       
   665 	  for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
       
   666 	    if ( (*capacity)[e] <= (*flow)[e] ) continue;
       
   667 	    Node w=g->tail(e);
       
   668 	    if ( level[w] == n && w != s ) {
       
   669 	      bfs_queue.push(w);
       
   670 	      Node z=level_list[l];
       
   671 	      if ( z!=INVALID ) left.set(z,w);
       
   672 	      right.set(w,z);
       
   673 	      level_list[l]=w;
       
   674 	      level.set(w, l);
       
   675 	    }
       
   676 	  }
       
   677 	  
       
   678 	  for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
       
   679 	    if ( 0 >= (*flow)[e] ) continue;
       
   680 	    Node w=g->head(e);
       
   681 	    if ( level[w] == n && w != s ) {
       
   682 	      bfs_queue.push(w);
       
   683 	      Node z=level_list[l];
       
   684 	      if ( z!=INVALID ) left.set(z,w);
       
   685 	      right.set(w,z);
       
   686 	      level_list[l]=w;
       
   687 	      level.set(w, l);
       
   688 	    }
       
   689 	  }
       
   690 	}
       
   691 	
       
   692 	//the starting flow
       
   693 	for(OutEdgeIt e(*g,s); e!=INVALID; ++e)
       
   694 	  {
       
   695 	    Num rem=(*capacity)[e]-(*flow)[e];
       
   696 	    if ( rem <= 0 ) continue;
       
   697 	    Node w=g->head(e);
       
   698 	    if ( level[w] < n ) {
       
   699 	      if ( excess[w] <= 0 && w!=t ) //putting into the stack
       
   700 		{
       
   701 		  next.set(w,first[level[w]]);
       
   702 		  first[level[w]]=w;
       
   703 		}   
       
   704 	      flow->set(e, (*capacity)[e]);
       
   705 	      excess.set(w, excess[w]+rem);
       
   706 	    }
       
   707 	  }
       
   708 	
       
   709 	for(InEdgeIt e(*g,s); e!=INVALID; ++e)
       
   710 	  {
       
   711 	    if ( (*flow)[e] <= 0 ) continue;
       
   712 	    Node w=g->tail(e);
       
   713 	    if ( level[w] < n ) {
       
   714 	      if ( excess[w] <= 0 && w!=t )
       
   715 		{
       
   716 		  next.set(w,first[level[w]]);
       
   717 		  first[level[w]]=w;
       
   718 		}  
       
   719 	      excess.set(w, excess[w]+(*flow)[e]);
       
   720 	      flow->set(e, 0);
       
   721 	    }
       
   722 	  }
       
   723 	break;
       
   724       case PRE_FLOW:
       
   725 	//Reverse_bfs from t in the residual graph,
       
   726 	//to find the starting level.
       
   727 	level.set(t,0);
       
   728 	bfs_queue.push(t);
       
   729 	
       
   730 	while (!bfs_queue.empty()) {
       
   731 	  
       
   732 	  Node v=bfs_queue.front();
       
   733 	  bfs_queue.pop();
       
   734 	  int l=level[v]+1;
       
   735 	  
       
   736 	  for(InEdgeIt e(*g,v) ; e!=INVALID; ++e) {
       
   737 	    if ( (*capacity)[e] <= (*flow)[e] ) continue;
       
   738 	    Node w=g->tail(e);
       
   739 	    if ( level[w] == n && w != s ) {
       
   740 	      bfs_queue.push(w);
       
   741 	      Node z=level_list[l];
       
   742 	      if ( z!=INVALID ) left.set(z,w);
       
   743 	      right.set(w,z);
       
   744 	      level_list[l]=w;
       
   745 	      level.set(w, l);
       
   746 	    }
       
   747 	  }
       
   748 	  
       
   749 	  for(OutEdgeIt e(*g,v) ; e!=INVALID; ++e) {
       
   750 	    if ( 0 >= (*flow)[e] ) continue;
       
   751 	    Node w=g->head(e);
       
   752 	    if ( level[w] == n && w != s ) {
       
   753 	      bfs_queue.push(w);
       
   754 	      Node z=level_list[l];
       
   755 	      if ( z!=INVALID ) left.set(z,w);
       
   756 	      right.set(w,z);
       
   757 	      level_list[l]=w;
       
   758 	      level.set(w, l);
       
   759 	    }
       
   760 	  }
       
   761 	}
       
   762 	
       
   763 	
       
   764 	//the starting flow
       
   765 	for(OutEdgeIt e(*g,s) ; e!=INVALID; ++e) {
       
   766 	  Num rem=(*capacity)[e]-(*flow)[e];
       
   767 	  if ( rem <= 0 ) continue;
       
   768 	  Node w=g->head(e);
       
   769 	  if ( level[w] < n ) {
       
   770 	    flow->set(e, (*capacity)[e]);
       
   771 	    excess.set(w, excess[w]+rem);
       
   772 	  }
       
   773 	}
       
   774 	
       
   775 	for(InEdgeIt e(*g,s) ; e!=INVALID; ++e) {
       
   776 	  if ( (*flow)[e] <= 0 ) continue;
       
   777 	  Node w=g->tail(e);
       
   778 	  if ( level[w] < n ) {
       
   779 	    excess.set(w, excess[w]+(*flow)[e]);
       
   780 	    flow->set(e, 0);
       
   781 	  }
       
   782 	}
       
   783 	
       
   784 	//computing the excess
       
   785 	for(NodeIt w(*g); w!=INVALID; ++w) {
       
   786 	  Num exc=0;
       
   787 	  
       
   788 	  for(InEdgeIt e(*g,w) ; e!=INVALID; ++e) exc+=(*flow)[e];
       
   789 	  for(OutEdgeIt e(*g,w) ; e!=INVALID; ++e) exc-=(*flow)[e];
       
   790 	  
       
   791 	  excess.set(w,exc);
       
   792 	  
       
   793 	  //putting the active nodes into the stack
       
   794 	  int lev=level[w];
       
   795 	    if ( exc > 0 && lev < n && Node(w) != t ) 
       
   796 	      ///\bug	    if ( exc > 0 && lev < n && w != t ) temporarily for working with wrappers. 
   666 	    {
   797 	    {
   667 	      Num c=(*capacity)[e];
   798 	      next.set(w,first[lev]);
   668 	      if ( c <= 0 ) continue;
   799 	      first[lev]=w;
   669 	      Node w=g->head(e);
   800 	    }
   670 	      if ( level[w] < n ) {
   801 	}
   671 		if ( excess[w] <= 0 && w!=t ) //putting into the stack
   802 	break;
   672 		  { 
   803       } //switch
   673 		    next.set(w,first[level[w]]);
       
   674 		    first[level[w]]=w;
       
   675 		  }
       
   676 		flow->set(e, c);
       
   677 		excess.set(w, excess[w]+c);
       
   678 	      }
       
   679 	    }
       
   680 	  break;
       
   681 	}
       
   682 
       
   683       case GEN_FLOW:
       
   684 	{
       
   685 	  //Reverse_bfs from t in the residual graph,
       
   686 	  //to find the starting level.
       
   687 	  level.set(t,0);
       
   688 	  bfs_queue.push(t);
       
   689 
       
   690 	  while (!bfs_queue.empty()) {
       
   691 
       
   692 	    Node v=bfs_queue.front();
       
   693 	    bfs_queue.pop();
       
   694 	    int l=level[v]+1;
       
   695 
       
   696 	    InEdgeIt e;
       
   697 	    for(g->first(e,v); g->valid(e); g->next(e)) {
       
   698 	      if ( (*capacity)[e] <= (*flow)[e] ) continue;
       
   699 	      Node w=g->tail(e);
       
   700 	      if ( level[w] == n && w != s ) {
       
   701 		bfs_queue.push(w);
       
   702 		Node z=level_list[l];
       
   703 		if ( g->valid(z) ) left.set(z,w);
       
   704 		right.set(w,z);
       
   705 		level_list[l]=w;
       
   706 		level.set(w, l);
       
   707 	      }
       
   708 	    }
       
   709 
       
   710 	    OutEdgeIt f;
       
   711 	    for(g->first(f,v); g->valid(f); g->next(f)) {
       
   712 	      if ( 0 >= (*flow)[f] ) continue;
       
   713 	      Node w=g->head(f);
       
   714 	      if ( level[w] == n && w != s ) {
       
   715 		bfs_queue.push(w);
       
   716 		Node z=level_list[l];
       
   717 		if ( g->valid(z) ) left.set(z,w);
       
   718 		right.set(w,z);
       
   719 		level_list[l]=w;
       
   720 		level.set(w, l);
       
   721 	      }
       
   722 	    }
       
   723 	  }
       
   724 
       
   725 	  //the starting flow
       
   726 	  OutEdgeIt e;
       
   727 	  for(g->first(e,s); g->valid(e); g->next(e))
       
   728 	    {
       
   729 	      Num rem=(*capacity)[e]-(*flow)[e];
       
   730 	      if ( rem <= 0 ) continue;
       
   731 	      Node w=g->head(e);
       
   732 	      if ( level[w] < n ) {
       
   733 		if ( excess[w] <= 0 && w!=t ) //putting into the stack
       
   734 		  {
       
   735 		    next.set(w,first[level[w]]);
       
   736 		    first[level[w]]=w;
       
   737 		  }   
       
   738 		flow->set(e, (*capacity)[e]);
       
   739 		excess.set(w, excess[w]+rem);
       
   740 	      }
       
   741 	    }
       
   742 
       
   743 	  InEdgeIt f;
       
   744 	  for(g->first(f,s); g->valid(f); g->next(f))
       
   745 	    {
       
   746 	      if ( (*flow)[f] <= 0 ) continue;
       
   747 	      Node w=g->tail(f);
       
   748 	      if ( level[w] < n ) {
       
   749 		if ( excess[w] <= 0 && w!=t )
       
   750 		  {
       
   751 		    next.set(w,first[level[w]]);
       
   752 		    first[level[w]]=w;
       
   753 		  }  
       
   754 		excess.set(w, excess[w]+(*flow)[f]);
       
   755 		flow->set(f, 0);
       
   756 	      }
       
   757 	    }
       
   758 	  break;
       
   759 	} //case GEN_FLOW
       
   760     
       
   761       case PRE_FLOW:
       
   762 	{
       
   763 	  //Reverse_bfs from t in the residual graph,
       
   764 	  //to find the starting level.
       
   765 	  level.set(t,0);
       
   766 	  bfs_queue.push(t);
       
   767 
       
   768 	  while (!bfs_queue.empty()) {
       
   769 
       
   770 	    Node v=bfs_queue.front();
       
   771 	    bfs_queue.pop();
       
   772 	    int l=level[v]+1;
       
   773 
       
   774 	    InEdgeIt e;
       
   775 	    for(g->first(e,v); g->valid(e); g->next(e)) {
       
   776 	      if ( (*capacity)[e] <= (*flow)[e] ) continue;
       
   777 	      Node w=g->tail(e);
       
   778 	      if ( level[w] == n && w != s ) {
       
   779 		bfs_queue.push(w);
       
   780 		Node z=level_list[l];
       
   781 		if ( g->valid(z) ) left.set(z,w);
       
   782 		right.set(w,z);
       
   783 		level_list[l]=w;
       
   784 		level.set(w, l);
       
   785 	      }
       
   786 	    }
       
   787 
       
   788 	    OutEdgeIt f;
       
   789 	    for(g->first(f,v); g->valid(f); g->next(f)) {
       
   790 	      if ( 0 >= (*flow)[f] ) continue;
       
   791 	      Node w=g->head(f);
       
   792 	      if ( level[w] == n && w != s ) {
       
   793 		bfs_queue.push(w);
       
   794 		Node z=level_list[l];
       
   795 		if ( g->valid(z) ) left.set(z,w);
       
   796 		right.set(w,z);
       
   797 		level_list[l]=w;
       
   798 		level.set(w, l);
       
   799 	      }
       
   800 	    }
       
   801 	  }
       
   802 
       
   803 
       
   804 	  //the starting flow
       
   805 	  OutEdgeIt e;
       
   806 	  for(g->first(e,s); g->valid(e); g->next(e))
       
   807 	    {
       
   808 	      Num rem=(*capacity)[e]-(*flow)[e];
       
   809 	      if ( rem <= 0 ) continue;
       
   810 	      Node w=g->head(e);
       
   811 	      if ( level[w] < n ) {
       
   812 		flow->set(e, (*capacity)[e]);
       
   813 		excess.set(w, excess[w]+rem);
       
   814 	      }
       
   815 	    }
       
   816 
       
   817 	  InEdgeIt f;
       
   818 	  for(g->first(f,s); g->valid(f); g->next(f))
       
   819 	    {
       
   820 	      if ( (*flow)[f] <= 0 ) continue;
       
   821 	      Node w=g->tail(f);
       
   822 	      if ( level[w] < n ) {
       
   823 		excess.set(w, excess[w]+(*flow)[f]);
       
   824 		flow->set(f, 0);
       
   825 	      }
       
   826 	    }
       
   827 	  
       
   828 	  NodeIt w; //computing the excess
       
   829 	  for(g->first(w); g->valid(w); g->next(w)) {
       
   830 	    Num exc=0;
       
   831 
       
   832 	    InEdgeIt e;
       
   833 	    for(g->first(e,w); g->valid(e); g->next(e)) exc+=(*flow)[e];
       
   834 	    OutEdgeIt f;
       
   835 	    for(g->first(f,w); g->valid(f); g->next(f)) exc-=(*flow)[f];
       
   836 
       
   837 	    excess.set(w,exc);
       
   838 
       
   839 	    //putting the active nodes into the stack
       
   840 	    int lev=level[w];
       
   841 	    if ( exc > 0 && lev < n && Node(w) != t ) 
       
   842 ///\bug	    if ( exc > 0 && lev < n && w != t ) tempararily for working with wrappers. in hugo 0.2 it will work. Amugy mukodik sage_graph-fal, de smart_graph-fal nem, azt hozzatennem.
       
   843 	      {
       
   844 		next.set(w,first[lev]);
       
   845 		first[lev]=w;
       
   846 	      }
       
   847 	  }
       
   848 	  break;
       
   849 	} //case PRE_FLOW
       
   850       }
       
   851     } //preflowPreproc
   804     } //preflowPreproc
   852 
   805 
   853 
   806 
   854     void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
   807     void relabel(Node w, int newlevel, NNMap& next, VecFirst& first,
   855 		 VecNode& level_list, NNMap& left,
   808 		 VecNode& level_list, NNMap& left,
   860 
   813 
   861       Node right_n=right[w];
   814       Node right_n=right[w];
   862       Node left_n=left[w];
   815       Node left_n=left[w];
   863 
   816 
   864       //unlacing starts
   817       //unlacing starts
   865       if ( g->valid(right_n) ) {
   818       if ( right_n!=INVALID ) {
   866 	if ( g->valid(left_n) ) {
   819 	if ( left_n!=INVALID ) {
   867 	  right.set(left_n, right_n);
   820 	  right.set(left_n, right_n);
   868 	  left.set(right_n, left_n);
   821 	  left.set(right_n, left_n);
   869 	} else {
   822 	} else {
   870 	  level_list[lev]=right_n;
   823 	  level_list[lev]=right_n;
   871 	  left.set(right_n, INVALID);
   824 	  left.set(right_n, INVALID);
   872 	}
   825 	}
   873       } else {
   826       } else {
   874 	if ( g->valid(left_n) ) {
   827 	if ( left_n!=INVALID ) {
   875 	  right.set(left_n, INVALID);
   828 	  right.set(left_n, INVALID);
   876 	} else {
   829 	} else {
   877 	  level_list[lev]=INVALID;
   830 	  level_list[lev]=INVALID;
   878 	}
   831 	}
   879       }
   832       }
   880       //unlacing ends
   833       //unlacing ends
   881 
   834 
   882       if ( !g->valid(level_list[lev]) ) {
   835       if ( level_list[lev]==INVALID ) {
   883 
   836 
   884 	//gapping starts
   837 	//gapping starts
   885 	for (int i=lev; i!=k ; ) {
   838 	for (int i=lev; i!=k ; ) {
   886 	  Node v=level_list[++i];
   839 	  Node v=level_list[++i];
   887 	  while ( g->valid(v) ) {
   840 	  while ( v!=INVALID ) {
   888 	    level.set(v,n);
   841 	    level.set(v,n);
   889 	    v=right[v];
   842 	    v=right[v];
   890 	  }
   843 	  }
   891 	  level_list[i]=INVALID;
   844 	  level_list[i]=INVALID;
   892 	  if ( !what_heur ) first[i]=INVALID;
   845 	  if ( !what_heur ) first[i]=INVALID;
   905 	  next.set(w,first[newlevel]);
   858 	  next.set(w,first[newlevel]);
   906 	  first[newlevel]=w;
   859 	  first[newlevel]=w;
   907 	  if ( what_heur ) b=newlevel;
   860 	  if ( what_heur ) b=newlevel;
   908 	  if ( k < newlevel ) ++k;      //now k=newlevel
   861 	  if ( k < newlevel ) ++k;      //now k=newlevel
   909 	  Node z=level_list[newlevel];
   862 	  Node z=level_list[newlevel];
   910 	  if ( g->valid(z) ) left.set(z,w);
   863 	  if ( z!=INVALID ) left.set(z,w);
   911 	  right.set(w,z);
   864 	  right.set(w,z);
   912 	  left.set(w,INVALID);
   865 	  left.set(w,INVALID);
   913 	  level_list[newlevel]=w;
   866 	  level_list[newlevel]=w;
   914 	}
   867 	}
   915       }
   868       }
   916     } //relabel
   869     } //relabel
   917 
   870 
   918     void printexcess() {////
   871     void printexcess() {////
   919       std::cout << "Excesses:" <<std::endl;
   872       std::cout << "Excesses:" <<std::endl;
   920 
   873 
   921       NodeIt v;
   874       for(NodeIt v(*g); v!=INVALID ; ++v) {
   922       for(g->first(v); g->valid(v); g->next(v)) {
       
   923 	std::cout << 1+(g->id(v)) << ":" << excess[v]<<std::endl; 
   875 	std::cout << 1+(g->id(v)) << ":" << excess[v]<<std::endl; 
   924       }
   876       }
   925     }
   877     }
   926 
   878 
   927  void printlevel() {////
   879     void printlevel() {////
   928       std::cout << "Levels:" <<std::endl;
   880       std::cout << "Levels:" <<std::endl;
   929 
   881 
   930       NodeIt v;
   882       for(NodeIt v(*g); v!=INVALID ; ++v) {
   931       for(g->first(v); g->valid(v); g->next(v)) {
       
   932 	std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl; 
   883 	std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl; 
   933       }
   884       }
   934     }
   885     }
   935 
   886 
   936 void printactive() {////
   887     void printactive() {////
   937       std::cout << "Levels:" <<std::endl;
   888       std::cout << "Levels:" <<std::endl;
   938 
   889 
   939       NodeIt v;
   890       for(NodeIt v(*g); v!=INVALID ; ++v) {
   940       for(g->first(v); g->valid(v); g->next(v)) {
       
   941 	std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl; 
   891 	std::cout << 1+(g->id(v)) << ":" << level[v]<<std::endl; 
   942       }
   892       }
   943     }
   893     }
   944 
   894 
   945 
   895