src/work/jacint/preflow.h
changeset 452 6636be9bc35e
parent 389 770cc1f4861f
child 465 d72e56f1730d
equal deleted inserted replaced
7:d70a9f19ba43 8:7efbde2c8673
     1 // -*- C++ -*-
     1 // -*- C++ -*-
     2 
       
     3 //run gyorsan tudna adni a minmincutot a 2 fazis elejen , ne vegyuk be konstruktorba egy cutmapet?
       
     4 //constzero jo igy?
       
     5 
       
     6 //majd marci megmondja betegyem-e bfs-t meg resgraphot
       
     7 
     2 
     8 /*
     3 /*
     9 Heuristics: 
     4 Heuristics: 
    10  2 phase
     5  2 phase
    11  gap
     6  gap
    12  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
    13  stack 'active' on the active nodes on level i implemented by hand
     8  stack 'active' on the active nodes on level i
    14  runs heuristic 'highest label' for H1*n relabels
     9  runs heuristic 'highest label' for H1*n relabels
    15  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'
    16  
    11  
    17 Parameters H0 and H1 are initialized to 20 and 10.
    12 Parameters H0 and H1 are initialized to 20 and 1.
    18 
    13 
    19 Constructors:
    14 Constructors:
    20 
    15 
    21 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 
    22      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
    26 void run()
    21 void run()
    27 
    22 
    28 T flowValue() : returns the value of a maximum flow
    23 T flowValue() : returns the value of a maximum flow
    29 
    24 
    30 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 
    31      minimum min cut. M should be a map of bools initialized to false.
    26      minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
    32 
    27 
    33 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 
    34      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.
    35 
    30 
    36 void minCut(CutMap& M) : sets M to the characteristic vector of 
    31 void minCut(CutMap& M) : sets M to the characteristic vector of 
    37      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.
    38 
    33 
    39 FIXME reset
       
    40 
       
    41 */
    34 */
    42 
    35 
    43 #ifndef HUGO_PREFLOW_H
    36 #ifndef HUGO_PREFLOW_H
    44 #define HUGO_PREFLOW_H
    37 #define HUGO_PREFLOW_H
    45 
    38 
    46 #define H0 20
    39 #define H0 20
    47 #define H1 1
    40 #define H1 1
    48 
    41 
    49 #include <vector>
    42 #include <vector>
    50 #include <queue>
    43 #include <queue>
       
    44 #include <stack>
    51 
    45 
    52 namespace hugo {
    46 namespace hugo {
    53 
    47 
    54   template <typename Graph, typename T, 
    48   template <typename Graph, typename T, 
    55 	    typename CapMap=typename Graph::template EdgeMap<T>, 
    49 	    typename CapMap=typename Graph::template EdgeMap<T>, 
    56             typename FlowMap=typename Graph::template EdgeMap<T> >
    50             typename FlowMap=typename Graph::template EdgeMap<T> >
    57   class Preflow {
    51   class Preflow {
    58     
    52     
    59     typedef typename Graph::Node Node;
    53     typedef typename Graph::Node Node;
    60     typedef typename Graph::Edge Edge;
       
    61     typedef typename Graph::NodeIt NodeIt;
    54     typedef typename Graph::NodeIt NodeIt;
    62     typedef typename Graph::OutEdgeIt OutEdgeIt;
    55     typedef typename Graph::OutEdgeIt OutEdgeIt;
    63     typedef typename Graph::InEdgeIt InEdgeIt;
    56     typedef typename Graph::InEdgeIt InEdgeIt;
    64     
    57 
       
    58     typedef typename std::vector<std::stack<Node> > VecStack;
       
    59     typedef typename Graph::template NodeMap<Node> NNMap;
       
    60     typedef typename std::vector<Node> VecNode;
       
    61 
    65     const Graph& G;
    62     const Graph& G;
    66     Node s;
    63     Node s;
    67     Node t;
    64     Node t;
    68     const CapMap& capacity;  
    65     CapMap* capacity;  
    69     FlowMap& flow;
    66     FlowMap* flow;
    70     T value;
    67     int n;      //the number of nodes of G
    71     bool constzero;
    68     typename Graph::template NodeMap<int> level;      
       
    69     typename Graph::template NodeMap<T> excess; 
       
    70 
    72 
    71 
    73   public:
    72   public:
       
    73  
       
    74     enum flowEnum{
       
    75       ZERO_FLOW=0,
       
    76       GEN_FLOW=1,
       
    77       PREFLOW=2
       
    78     };
       
    79 
    74     Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity, 
    80     Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity, 
    75 	    FlowMap& _flow, bool _constzero ) :
    81 	    FlowMap& _flow) :
    76       G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero) {}
    82       G(_G), s(_s), t(_t), capacity(&_capacity), 
       
    83       flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
       
    84 
       
    85     void run() {
       
    86       preflow( ZERO_FLOW );
       
    87     }
    77     
    88     
    78     
    89     void preflow( flowEnum fe ) {
    79     void run() {
    90       preflowPhase0(fe);
    80       
    91       preflowPhase1();
    81       value=0;                //for the subsequent runs
    92     }
    82 
    93 
    83       bool phase=0;        //phase 0 is the 1st phase, phase 1 is the 2nd
    94     void preflowPhase0( flowEnum fe ) {
    84       int n=G.nodeNum(); 
    95       
    85       int heur0=(int)(H0*n);  //time while running 'bound decrease' 
    96       int heur0=(int)(H0*n);  //time while running 'bound decrease' 
    86       int heur1=(int)(H1*n);  //time while running 'highest label'
    97       int heur1=(int)(H1*n);  //time while running 'highest label'
    87       int heur=heur1;         //starting time interval (#of relabels)
    98       int heur=heur1;         //starting time interval (#of relabels)
       
    99       int numrelabel=0;
       
   100      
    88       bool what_heur=1;       
   101       bool what_heur=1;       
    89       /*
   102       //It is 0 in case 'bound decrease' and 1 in case 'highest label'
    90 	what_heur is 0 in case 'bound decrease' 
   103 
    91 	and 1 in case 'highest label'
       
    92       */
       
    93       bool end=false;     
   104       bool end=false;     
    94       /*
   105       //Needed for 'bound decrease', true means no active nodes are above bound b.
    95 	Needed for 'bound decrease', 'true'
   106 
    96 	means no active nodes are above bound b.
       
    97       */
       
    98       int relabel=0;
       
    99       int k=n-2;  //bound on the highest level under n containing a node
   107       int k=n-2;  //bound on the highest level under n containing a node
   100       int b=k;    //bound on the highest level under n of an active node
   108       int b=k;    //bound on the highest level under n of an active node
   101       
   109       
   102       typename Graph::template NodeMap<int> level(G,n);      
   110       VecStack active(n);
   103       typename Graph::template NodeMap<T> excess(G); 
   111       
   104 
   112       NNMap left(G,INVALID);
   105       std::vector<Node> active(n-1,INVALID);
   113       NNMap right(G,INVALID);
   106       typename Graph::template NodeMap<Node> next(G,INVALID);
   114       VecNode level_list(n,INVALID);
   107       //Stack of the active nodes in level i < n.
   115       //List of the nodes in level i<n, set to n.
   108       //We use it in both phases.
   116 
   109 
   117       NodeIt v;
   110       typename Graph::template NodeMap<Node> left(G,INVALID);
   118       for(G.first(v); G.valid(v); G.next(v)) level.set(v,n);
   111       typename Graph::template NodeMap<Node> right(G,INVALID);
   119       //setting each node to level n
   112       std::vector<Node> level_list(n,INVALID);
   120       
   113       /*
   121       switch ( fe ) {
   114 	List of the nodes in level i<n.
   122       case PREFLOW:
   115       */
   123 	{
   116 
   124 	  //counting the excess
   117 
   125 	  NodeIt v;
   118       if ( constzero ) {
   126 	  for(G.first(v); G.valid(v); G.next(v)) {
   119      
   127 	    T exc=0;
   120 	/*Reverse_bfs from t, to find the starting level.*/
   128 	  
   121 	level.set(t,0);
   129 	    InEdgeIt e;
   122 	std::queue<Node> bfs_queue;
   130 	    for(G.first(e,v); G.valid(e); G.next(e)) exc+=(*flow)[e];
   123 	bfs_queue.push(t);
   131 	    OutEdgeIt f;
   124 	
   132 	    for(G.first(f,v); G.valid(f); G.next(f)) exc-=(*flow)[f];
   125 	while (!bfs_queue.empty()) {
   133 	    
   126 	  
   134 	    excess.set(v,exc);	  
   127 	  Node v=bfs_queue.front();	
   135 	    
   128 	  bfs_queue.pop();
   136 	    //putting the active nodes into the stack
   129 	  int l=level[v]+1;
   137 	    int lev=level[v];
       
   138 	    if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
       
   139 	  }
       
   140 	  break;
       
   141 	}
       
   142       case GEN_FLOW:
       
   143 	{
       
   144 	  //Counting the excess of t
       
   145 	  T exc=0;
   130 	  
   146 	  
   131 	  InEdgeIt e;
   147 	  InEdgeIt e;
   132 	  for(G.first(e,v); G.valid(e); G.next(e)) {
   148 	  for(G.first(e,t); G.valid(e); G.next(e)) exc+=(*flow)[e];
   133 	    Node w=G.tail(e);
   149 	  OutEdgeIt f;
   134 	    if ( level[w] == n && w != s ) {
   150 	  for(G.first(f,t); G.valid(f); G.next(f)) exc-=(*flow)[f];
   135 	      bfs_queue.push(w);
   151 	  
   136 	      Node first=level_list[l];
   152 	  excess.set(t,exc);	
   137 	      if ( G.valid(first) ) left.set(first,w);
   153 	  
   138 	      right.set(w,first);
   154 	  break;
   139 	      level_list[l]=w;
       
   140 	      level.set(w, l);
       
   141 	    }
       
   142 	  }
       
   143 	}
       
   144 
       
   145 	//the starting flow
       
   146 	OutEdgeIt e;
       
   147 	for(G.first(e,s); G.valid(e); G.next(e)) 
       
   148 	{
       
   149 	  T c=capacity[e];
       
   150 	  if ( c == 0 ) continue;
       
   151 	  Node w=G.head(e);
       
   152 	  if ( level[w] < n ) {	  
       
   153 	    if ( excess[w] == 0 && w!=t ) {
       
   154 	      next.set(w,active[level[w]]);
       
   155 	      active[level[w]]=w;
       
   156 	    }
       
   157 	    flow.set(e, c); 
       
   158 	    excess.set(w, excess[w]+c);
       
   159 	  }
       
   160 	}
   155 	}
   161       }
   156       }
   162       else 
   157       
   163       {
   158       preflowPreproc( fe, active, level_list, left, right );
   164 	
   159       //End of preprocessing 
   165 	/*
   160       
   166 	  Reverse_bfs from t in the residual graph, 
   161       
   167 	  to find the starting level.
   162       //Push/relabel on the highest level active nodes.
   168 	*/
       
   169 	level.set(t,0);
       
   170 	std::queue<Node> bfs_queue;
       
   171 	bfs_queue.push(t);
       
   172 	
       
   173 	while (!bfs_queue.empty()) {
       
   174 	  
       
   175 	  Node v=bfs_queue.front();	
       
   176 	  bfs_queue.pop();
       
   177 	  int l=level[v]+1;
       
   178 	  
       
   179 	  InEdgeIt e;
       
   180 	  for(G.first(e,v); G.valid(e); G.next(e)) {
       
   181 	    if ( capacity[e] == flow[e] ) continue;
       
   182 	    Node w=G.tail(e);
       
   183 	    if ( level[w] == n && w != s ) {
       
   184 	      bfs_queue.push(w);
       
   185 	      Node first=level_list[l];
       
   186 	      if ( G.valid(first) ) left.set(first,w);
       
   187 	      right.set(w,first);
       
   188 	      level_list[l]=w;
       
   189 	      level.set(w, l);
       
   190 	    }
       
   191 	  }
       
   192 	    
       
   193 	  OutEdgeIt f;
       
   194 	  for(G.first(f,v); G.valid(f); G.next(f)) {
       
   195 	    if ( 0 == flow[f] ) continue;
       
   196 	    Node w=G.head(f);
       
   197 	    if ( level[w] == n && w != s ) {
       
   198 	      bfs_queue.push(w);
       
   199 	      Node first=level_list[l];
       
   200 	      if ( G.valid(first) ) left.set(first,w);
       
   201 	      right.set(w,first);
       
   202 	      level_list[l]=w;
       
   203 	      level.set(w, l);
       
   204 	    }
       
   205 	  }
       
   206 	}
       
   207       
       
   208 	
       
   209 	/*
       
   210 	  Counting the excess
       
   211 	*/
       
   212 	NodeIt v;
       
   213 	for(G.first(v); G.valid(v); G.next(v)) {
       
   214 	  T exc=0;
       
   215 
       
   216 	  InEdgeIt e;
       
   217 	  for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e];
       
   218 	  OutEdgeIt f;
       
   219 	  for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[e];
       
   220 
       
   221 	  excess.set(v,exc);	  
       
   222 
       
   223 	  //putting the active nodes into the stack
       
   224 	  int lev=level[v];
       
   225 	  if ( exc > 0 && lev < n ) {
       
   226 	    next.set(v,active[lev]);
       
   227 	    active[lev]=v;
       
   228 	  }
       
   229 	}
       
   230 
       
   231 
       
   232 	//the starting flow
       
   233 	OutEdgeIt e;
       
   234 	for(G.first(e,s); G.valid(e); G.next(e)) 
       
   235 	{
       
   236 	  T rem=capacity[e]-flow[e];
       
   237 	  if ( rem == 0 ) continue;
       
   238 	  Node w=G.head(e);
       
   239 	  if ( level[w] < n ) {	  
       
   240 	    if ( excess[w] == 0 && w!=t ) {
       
   241 	      next.set(w,active[level[w]]);
       
   242 	      active[level[w]]=w;
       
   243 	    }
       
   244 	    flow.set(e, capacity[e]); 
       
   245 	    excess.set(w, excess[w]+rem);
       
   246 	  }
       
   247 	}
       
   248 	
       
   249 	InEdgeIt f;
       
   250 	for(G.first(f,s); G.valid(f); G.next(f)) 
       
   251 	{
       
   252 	  if ( flow[f] == 0 ) continue;
       
   253 	  Node w=G.head(f);
       
   254 	  if ( level[w] < n ) {	  
       
   255 	    if ( excess[w] == 0 && w!=t ) {
       
   256 	      next.set(w,active[level[w]]);
       
   257 	      active[level[w]]=w;
       
   258 	    }
       
   259 	    excess.set(w, excess[w]+flow[f]);
       
   260 	    flow.set(f, 0); 
       
   261 	  }
       
   262 	}
       
   263       }
       
   264 
       
   265 
       
   266 
       
   267 
       
   268       /* 
       
   269 	 End of preprocessing 
       
   270       */
       
   271 
       
   272 
       
   273 
       
   274       /*
       
   275 	Push/relabel on the highest level active nodes.
       
   276       */	
       
   277       while ( true ) {
   163       while ( true ) {
   278 	
       
   279 	if ( b == 0 ) {
   164 	if ( b == 0 ) {
   280 	  if ( phase ) break;
       
   281 	  
       
   282 	  if ( !what_heur && !end && k > 0 ) {
   165 	  if ( !what_heur && !end && k > 0 ) {
   283 	    b=k;
   166 	    b=k;
   284 	    end=true;
   167 	    end=true;
   285 	  } else {
   168 	  } else break;
   286 	    phase=1;
   169 	}
   287 	    level.set(s,0);
   170 	
   288 	    std::queue<Node> bfs_queue;
   171 	if ( active[b].empty() ) --b; 
   289 	    bfs_queue.push(s);
       
   290 	    
       
   291 	    while (!bfs_queue.empty()) {
       
   292 	      
       
   293 	      Node v=bfs_queue.front();	
       
   294 	      bfs_queue.pop();
       
   295 	      int l=level[v]+1;
       
   296 	      
       
   297 	      InEdgeIt e;
       
   298 	      for(G.first(e,v); G.valid(e); G.next(e)) {
       
   299 		if ( capacity[e] == flow[e] ) continue;
       
   300 		Node u=G.tail(e);
       
   301 		if ( level[u] >= n ) { 
       
   302 		  bfs_queue.push(u);
       
   303 		  level.set(u, l);
       
   304 		  if ( excess[u] > 0 ) {
       
   305 		    next.set(u,active[l]);
       
   306 		    active[l]=u;
       
   307 		  }
       
   308 		}
       
   309 	      }
       
   310 	    
       
   311 	      OutEdgeIt f;
       
   312 	      for(G.first(f,v); G.valid(f); G.next(f)) {
       
   313 		if ( 0 == flow[f] ) continue;
       
   314 		Node u=G.head(f);
       
   315 		if ( level[u] >= n ) { 
       
   316 		  bfs_queue.push(u);
       
   317 		  level.set(u, l);
       
   318 		  if ( excess[u] > 0 ) {
       
   319 		    next.set(u,active[l]);
       
   320 		    active[l]=u;
       
   321 		  }
       
   322 		}
       
   323 	      }
       
   324 	    }
       
   325 	    b=n-2;
       
   326 	    }
       
   327 	    
       
   328 	}
       
   329 	  
       
   330 	  
       
   331 	if ( !G.valid(active[b]) ) --b; 
       
   332 	else {
   172 	else {
   333 	  end=false;  
   173 	  end=false;  
   334 
   174 	  Node w=active[b].top();
   335 	  Node w=active[b];
   175 	  active[b].pop();
   336 	  active[b]=next[w];
   176 	  int newlevel=push(w,active);
   337 	  int lev=level[w];
   177 	  if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, 
   338 	  T exc=excess[w];
   178 				       left, right, b, k, what_heur);
   339 	  int newlevel=n;       //bound on the next level of w
   179 	  
   340 	  
   180 	  ++numrelabel; 
   341 	  OutEdgeIt e;
   181 	  if ( numrelabel >= heur ) {
   342 	  for(G.first(e,w); G.valid(e); G.next(e)) {
   182 	    numrelabel=0;
   343 	    
   183 	    if ( what_heur ) {
   344 	    if ( flow[e] == capacity[e] ) continue; 
   184 	      what_heur=0;
   345 	    Node v=G.head(e);            
   185 	      heur=heur0;
   346 	    //e=wv	    
   186 	      end=false;
   347 	    
   187 	    } else {
   348 	    if( lev > level[v] ) {      
   188 	      what_heur=1;
   349 	      /*Push is allowed now*/
   189 	      heur=heur1;
       
   190 	      b=k; 
       
   191 	    }
       
   192 	  }
       
   193 	} 
       
   194       } 
       
   195     }
       
   196 
       
   197 
       
   198     void preflowPhase1() {
       
   199       
       
   200       int k=n-2;  //bound on the highest level under n containing a node
       
   201       int b=k;    //bound on the highest level under n of an active node
       
   202       
       
   203       VecStack active(n);
       
   204       level.set(s,0);
       
   205       std::queue<Node> bfs_queue;
       
   206       bfs_queue.push(s);
       
   207 	    
       
   208       while (!bfs_queue.empty()) {
       
   209 	
       
   210 	Node v=bfs_queue.front();	
       
   211 	bfs_queue.pop();
       
   212 	int l=level[v]+1;
   350 	      
   213 	      
   351 	      if ( excess[v]==0 && v!=t && v!=s ) {
   214 	InEdgeIt e;
   352 		int lev_v=level[v];
   215 	for(G.first(e,v); G.valid(e); G.next(e)) {
   353 		next.set(v,active[lev_v]);
   216 	  if ( (*capacity)[e] == (*flow)[e] ) continue;
   354 		active[lev_v]=v;
   217 	  Node u=G.tail(e);
   355 	      }
   218 	  if ( level[u] >= n ) { 
   356 	      
   219 	    bfs_queue.push(u);
   357 	      T cap=capacity[e];
   220 	    level.set(u, l);
   358 	      T flo=flow[e];
   221 	    if ( excess[u] > 0 ) active[l].push(u);
   359 	      T remcap=cap-flo;
   222 	  }
   360 	      
   223 	}
   361 	      if ( remcap >= exc ) {       
   224 	
   362 		/*A nonsaturating push.*/
   225 	OutEdgeIt f;
   363 		
   226 	for(G.first(f,v); G.valid(f); G.next(f)) {
   364 		flow.set(e, flo+exc);
   227 	  if ( 0 == (*flow)[f] ) continue;
   365 		excess.set(v, excess[v]+exc);
   228 	  Node u=G.head(f);
   366 		exc=0;
   229 	  if ( level[u] >= n ) { 
   367 		break; 
   230 	    bfs_queue.push(u);
   368 		
   231 	    level.set(u, l);
   369 	      } else { 
   232 	    if ( excess[u] > 0 ) active[l].push(u);
   370 		/*A saturating push.*/
   233 	  }
   371 		
   234 	}
   372 		flow.set(e, cap);
   235       }
   373 		excess.set(v, excess[v]+remcap);
   236       b=n-2;
   374 		exc-=remcap;
   237 
   375 	      }
   238       while ( true ) {
   376 	    } else if ( newlevel > level[v] ){
   239 	
   377 	      newlevel = level[v];
   240 	if ( b == 0 ) break;
   378 	    }	    
   241 
   379 	    
   242 	if ( active[b].empty() ) --b; 
   380 	  } //for out edges wv 
   243 	else {
   381 	
   244 	  Node w=active[b].top();
   382 	
   245 	  active[b].pop();
   383 	if ( exc > 0 ) {	
   246 	  int newlevel=push(w,active);	  
   384 	  InEdgeIt e;
   247 
   385 	  for(G.first(e,w); G.valid(e); G.next(e)) {
   248 	  //relabel
   386 	    
   249 	  if ( excess[w] > 0 ) {
   387 	    if( flow[e] == 0 ) continue; 
       
   388 	    Node v=G.tail(e);  
       
   389 	    //e=vw
       
   390 	    
       
   391 	    if( lev > level[v] ) {  
       
   392 	      /*Push is allowed now*/
       
   393 	      
       
   394 	      if ( excess[v]==0 && v!=t && v!=s ) {
       
   395 		int lev_v=level[v];
       
   396 		next.set(v,active[lev_v]);
       
   397 		active[lev_v]=v;
       
   398 	      }
       
   399 	      
       
   400 	      T flo=flow[e];
       
   401 	      
       
   402 	      if ( flo >= exc ) { 
       
   403 		/*A nonsaturating push.*/
       
   404 		
       
   405 		flow.set(e, flo-exc);
       
   406 		excess.set(v, excess[v]+exc);
       
   407 		exc=0;
       
   408 		break; 
       
   409 	      } else {                                               
       
   410 		/*A saturating push.*/
       
   411 		
       
   412 		excess.set(v, excess[v]+flo);
       
   413 		exc-=flo;
       
   414 		flow.set(e,0);
       
   415 	      }  
       
   416 	    } else if ( newlevel > level[v] ) {
       
   417 	      newlevel = level[v];
       
   418 	    }	    
       
   419 	  } //for in edges vw
       
   420 	  
       
   421 	} // if w still has excess after the out edge for cycle
       
   422 	
       
   423 	excess.set(w, exc);
       
   424 	 
       
   425 	/*
       
   426 	  Relabel
       
   427 	*/
       
   428 	
       
   429 
       
   430 	if ( exc > 0 ) {
       
   431 	  //now 'lev' is the old level of w
       
   432 	
       
   433 	  if ( phase ) {
       
   434 	    level.set(w,++newlevel);
   250 	    level.set(w,++newlevel);
   435 	    next.set(w,active[newlevel]);
   251 	    active[newlevel].push(w);
   436 	    active[newlevel]=w;
       
   437 	    b=newlevel;
   252 	    b=newlevel;
   438 	  } else {
   253 	  }
   439 	    //unlacing starts
       
   440 	    Node right_n=right[w];
       
   441 	    Node left_n=left[w];
       
   442 
       
   443 	    if ( G.valid(right_n) ) {
       
   444 	      if ( G.valid(left_n) ) {
       
   445 		right.set(left_n, right_n);
       
   446 		left.set(right_n, left_n);
       
   447 	      } else {
       
   448 		level_list[lev]=right_n;   
       
   449 		left.set(right_n, INVALID);
       
   450 	      } 
       
   451 	    } else {
       
   452 	      if ( G.valid(left_n) ) {
       
   453 		right.set(left_n, INVALID);
       
   454 	      } else { 
       
   455 		level_list[lev]=INVALID;   
       
   456 	      } 
       
   457 	    } 
       
   458 	    //unlacing ends
       
   459 		
       
   460 	    if ( !G.valid(level_list[lev]) ) {
       
   461 	      
       
   462 	       //gapping starts
       
   463 	      for (int i=lev; i!=k ; ) {
       
   464 		Node v=level_list[++i];
       
   465 		while ( G.valid(v) ) {
       
   466 		  level.set(v,n);
       
   467 		  v=right[v];
       
   468 		}
       
   469 		level_list[i]=INVALID;
       
   470 		if ( !what_heur ) active[i]=INVALID;
       
   471 	      }	     
       
   472 
       
   473 	      level.set(w,n);
       
   474 	      b=lev-1;
       
   475 	      k=b;
       
   476 	      //gapping ends
       
   477 	    
       
   478 	    } else {
       
   479 	      
       
   480 	      if ( newlevel == n ) level.set(w,n); 
       
   481 	      else {
       
   482 		level.set(w,++newlevel);
       
   483 		next.set(w,active[newlevel]);
       
   484 		active[newlevel]=w;
       
   485 		if ( what_heur ) b=newlevel;
       
   486 		if ( k < newlevel ) ++k;      //now k=newlevel
       
   487 		Node first=level_list[newlevel];
       
   488 		if ( G.valid(first) ) left.set(first,w);
       
   489 		right.set(w,first);
       
   490 		left.set(w,INVALID);
       
   491 		level_list[newlevel]=w;
       
   492 	      }
       
   493 	    }
       
   494 
       
   495 
       
   496 	    ++relabel; 
       
   497 	    if ( relabel >= heur ) {
       
   498 	      relabel=0;
       
   499 	      if ( what_heur ) {
       
   500 		what_heur=0;
       
   501 		heur=heur0;
       
   502 		end=false;
       
   503 	      } else {
       
   504 		what_heur=1;
       
   505 		heur=heur1;
       
   506 		b=k; 
       
   507 	      }
       
   508 	    }
       
   509 	  } //phase 0
       
   510 	  
       
   511 	  
       
   512 	} // if ( exc > 0 )
       
   513 	  
       
   514 	
       
   515 	}  // if stack[b] is nonempty
   254 	}  // if stack[b] is nonempty
   516 	
       
   517       } // while(true)
   255       } // while(true)
   518 
   256     }
   519 
   257 
   520       value = excess[t];
   258 
   521       /*Max flow value.*/
   259     //Returns the maximum value of a flow.
   522      
       
   523     } //void run()
       
   524 
       
   525 
       
   526 
       
   527 
       
   528 
       
   529     /*
       
   530       Returns the maximum value of a flow.
       
   531      */
       
   532 
       
   533     T flowValue() {
   260     T flowValue() {
   534       return value;
   261       return excess[t];
   535     }
   262     }
   536 
   263 
   537 
   264     //should be used only between preflowPhase0 and preflowPhase1
   538     FlowMap Flow() {
   265     template<typename _CutMap>
   539       return flow;
   266     void actMinCut(_CutMap& M) {
   540       }
       
   541 
       
   542 
       
   543     
       
   544     void Flow(FlowMap& _flow ) {
       
   545       NodeIt v;
   267       NodeIt v;
   546       for(G.first(v) ; G.valid(v); G.next(v))
   268       for(G.first(v); G.valid(v); G.next(v)) 
   547 	_flow.set(v,flow[v]);
   269 	if ( level[v] < n ) M.set(v,false);
       
   270 	else M.set(v,true);
   548     }
   271     }
   549 
   272 
   550 
   273 
   551 
   274 
   552     /*
   275     /*
   553       Returns the minimum min cut, by a bfs from s in the residual graph.
   276       Returns the minimum min cut, by a bfs from s in the residual graph.
   554     */
   277     */
   555    
       
   556     template<typename _CutMap>
   278     template<typename _CutMap>
   557     void minMinCut(_CutMap& M) {
   279     void minMinCut(_CutMap& M) {
   558     
   280     
   559       std::queue<Node> queue;
   281       std::queue<Node> queue;
   560       
   282       
   566 	queue.pop();
   288 	queue.pop();
   567 
   289 
   568 	OutEdgeIt e;
   290 	OutEdgeIt e;
   569 	for(G.first(e,w) ; G.valid(e); G.next(e)) {
   291 	for(G.first(e,w) ; G.valid(e); G.next(e)) {
   570 	  Node v=G.head(e);
   292 	  Node v=G.head(e);
   571 	  if (!M[v] && flow[e] < capacity[e] ) {
   293 	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
   572 	    queue.push(v);
   294 	    queue.push(v);
   573 	    M.set(v, true);
   295 	    M.set(v, true);
   574 	  }
   296 	  }
   575 	} 
   297 	} 
   576 
   298 
   577 	InEdgeIt f;
   299 	InEdgeIt f;
   578 	for(G.first(f,w) ; G.valid(f); G.next(f)) {
   300 	for(G.first(f,w) ; G.valid(f); G.next(f)) {
   579 	  Node v=G.tail(f);
   301 	  Node v=G.tail(f);
   580 	  if (!M[v] && flow[f] > 0 ) {
   302 	  if (!M[v] && (*flow)[f] > 0 ) {
   581 	    queue.push(v);
   303 	    queue.push(v);
   582 	    M.set(v, true);
   304 	    M.set(v, true);
   583 	  }
   305 	  }
   584 	} 
   306 	} 
   585       }
   307       }
   592       from t in the residual graph.
   314       from t in the residual graph.
   593     */
   315     */
   594     
   316     
   595     template<typename _CutMap>
   317     template<typename _CutMap>
   596     void maxMinCut(_CutMap& M) {
   318     void maxMinCut(_CutMap& M) {
   597     
   319 
       
   320       NodeIt v;
       
   321       for(G.first(v) ; G.valid(v); G.next(v)) {
       
   322 	M.set(v, true);
       
   323       }
       
   324 
   598       std::queue<Node> queue;
   325       std::queue<Node> queue;
   599       
   326       
   600       M.set(t,true);        
   327       M.set(t,false);        
   601       queue.push(t);
   328       queue.push(t);
   602 
   329 
   603       while (!queue.empty()) {
   330       while (!queue.empty()) {
   604         Node w=queue.front();
   331         Node w=queue.front();
   605 	queue.pop();
   332 	queue.pop();
   606 
   333 
   607 
   334 
   608 	InEdgeIt e;
   335 	InEdgeIt e;
   609 	for(G.first(e,w) ; G.valid(e); G.next(e)) {
   336 	for(G.first(e,w) ; G.valid(e); G.next(e)) {
   610 	  Node v=G.tail(e);
   337 	  Node v=G.tail(e);
   611 	  if (!M[v] && flow[e] < capacity[e] ) {
   338 	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
   612 	    queue.push(v);
   339 	    queue.push(v);
   613 	    M.set(v, true);
   340 	    M.set(v, false);
   614 	  }
   341 	  }
   615 	}
   342 	}
   616 	
   343 	
   617 	OutEdgeIt f;
   344 	OutEdgeIt f;
   618 	for(G.first(f,w) ; G.valid(f); G.next(f)) {
   345 	for(G.first(f,w) ; G.valid(f); G.next(f)) {
   619 	  Node v=G.head(f);
   346 	  Node v=G.head(f);
   620 	  if (!M[v] && flow[f] > 0 ) {
   347 	  if (M[v] && (*flow)[f] > 0 ) {
   621 	    queue.push(v);
   348 	    queue.push(v);
   622 	    M.set(v, true);
   349 	    M.set(v, false);
   623 	  }
   350 	  }
   624 	}
   351 	}
   625       }
   352       }
   626 
   353     }
   627       NodeIt v;
       
   628       for(G.first(v) ; G.valid(v); G.next(v)) {
       
   629 	M.set(v, !M[v]);
       
   630       }
       
   631 
       
   632     }
       
   633 
       
   634 
   354 
   635 
   355 
   636     template<typename CutMap>
   356     template<typename CutMap>
   637     void minCut(CutMap& M) {
   357     void minCut(CutMap& M) {
   638       minMinCut(M);
   358       minMinCut(M);
   639     }
   359     }
   640 
   360 
   641     
   361     
   642     void reset_target (Node _t) {t=_t;}
   362     void resetTarget (const Node _t) {t=_t;}
   643     void reset_source (Node _s) {s=_s;}
   363 
       
   364     void resetSource (const Node _s) {s=_s;}
   644    
   365    
   645     template<typename _CapMap>   
   366     void resetCap (const CapMap& _cap) {
   646     void reset_cap (_CapMap _cap) {capacity=_cap;}
   367       capacity=&_cap;
   647 
   368     }
   648     template<typename _FlowMap>   
   369     
   649     void reset_cap (_FlowMap _flow, bool _constzero) {
   370     void resetFlow (FlowMap& _flow) {
   650       flow=_flow;
   371       flow=&_flow;
   651       constzero=_constzero;
   372     }
   652     }
   373 
   653 
   374 
   654 
   375   private:
       
   376 
       
   377     int push(const Node w, VecStack& active) {
       
   378       
       
   379       int lev=level[w];
       
   380       T exc=excess[w];
       
   381       int newlevel=n;       //bound on the next level of w
       
   382 	  
       
   383       OutEdgeIt e;
       
   384       for(G.first(e,w); G.valid(e); G.next(e)) {
       
   385 	    
       
   386 	if ( (*flow)[e] == (*capacity)[e] ) continue; 
       
   387 	Node v=G.head(e);            
       
   388 	    
       
   389 	if( lev > level[v] ) { //Push is allowed now
       
   390 	  
       
   391 	  if ( excess[v]==0 && v!=t && v!=s ) {
       
   392 	    int lev_v=level[v];
       
   393 	    active[lev_v].push(v);
       
   394 	  }
       
   395 	  
       
   396 	  T cap=(*capacity)[e];
       
   397 	  T flo=(*flow)[e];
       
   398 	  T remcap=cap-flo;
       
   399 	  
       
   400 	  if ( remcap >= exc ) { //A nonsaturating push.
       
   401 	    
       
   402 	    flow->set(e, flo+exc);
       
   403 	    excess.set(v, excess[v]+exc);
       
   404 	    exc=0;
       
   405 	    break; 
       
   406 	    
       
   407 	  } else { //A saturating push.
       
   408 	    flow->set(e, cap);
       
   409 	    excess.set(v, excess[v]+remcap);
       
   410 	    exc-=remcap;
       
   411 	  }
       
   412 	} else if ( newlevel > level[v] ) newlevel = level[v];
       
   413       } //for out edges wv 
       
   414       
       
   415       if ( exc > 0 ) {	
       
   416 	InEdgeIt e;
       
   417 	for(G.first(e,w); G.valid(e); G.next(e)) {
       
   418 	  
       
   419 	  if( (*flow)[e] == 0 ) continue; 
       
   420 	  Node v=G.tail(e); 
       
   421 	  
       
   422 	  if( lev > level[v] ) { //Push is allowed now
       
   423 	    
       
   424 	    if ( excess[v]==0 && v!=t && v!=s ) {
       
   425 	      int lev_v=level[v];
       
   426 	      active[lev_v].push(v);
       
   427 	    }
       
   428 	    
       
   429 	    T flo=(*flow)[e];
       
   430 	    
       
   431 	    if ( flo >= exc ) { //A nonsaturating push.
       
   432 	      
       
   433 	      flow->set(e, flo-exc);
       
   434 	      excess.set(v, excess[v]+exc);
       
   435 	      exc=0;
       
   436 	      break; 
       
   437 	    } else {  //A saturating push.
       
   438 	      
       
   439 	      excess.set(v, excess[v]+flo);
       
   440 	      exc-=flo;
       
   441 	      flow->set(e,0);
       
   442 	    }  
       
   443 	  } else if ( newlevel > level[v] ) newlevel = level[v];
       
   444 	} //for in edges vw
       
   445 	
       
   446       } // if w still has excess after the out edge for cycle
       
   447       
       
   448       excess.set(w, exc);
       
   449       
       
   450       return newlevel;
       
   451      }
       
   452 
       
   453 
       
   454     void preflowPreproc ( flowEnum fe, VecStack& active, 
       
   455 			  VecNode& level_list, NNMap& left, NNMap& right ) {
       
   456 
       
   457       std::queue<Node> bfs_queue;
       
   458       
       
   459       switch ( fe ) {
       
   460       case ZERO_FLOW: 
       
   461 	{
       
   462 	  //Reverse_bfs from t, to find the starting level.
       
   463 	  level.set(t,0);
       
   464 	  bfs_queue.push(t);
       
   465 	
       
   466 	  while (!bfs_queue.empty()) {
       
   467 	    
       
   468 	    Node v=bfs_queue.front();	
       
   469 	    bfs_queue.pop();
       
   470 	    int l=level[v]+1;
       
   471 	    
       
   472 	    InEdgeIt e;
       
   473 	    for(G.first(e,v); G.valid(e); G.next(e)) {
       
   474 	      Node w=G.tail(e);
       
   475 	      if ( level[w] == n && w != s ) {
       
   476 		bfs_queue.push(w);
       
   477 		Node first=level_list[l];
       
   478 		if ( G.valid(first) ) left.set(first,w);
       
   479 		right.set(w,first);
       
   480 		level_list[l]=w;
       
   481 		level.set(w, l);
       
   482 	      }
       
   483 	    }
       
   484 	  }
       
   485 	  
       
   486 	  //the starting flow
       
   487 	  OutEdgeIt e;
       
   488 	  for(G.first(e,s); G.valid(e); G.next(e)) 
       
   489 	    {
       
   490 	      T c=(*capacity)[e];
       
   491 	      if ( c == 0 ) continue;
       
   492 	      Node w=G.head(e);
       
   493 	      if ( level[w] < n ) {	  
       
   494 		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
       
   495 		flow->set(e, c); 
       
   496 		excess.set(w, excess[w]+c);
       
   497 	      }
       
   498 	    }
       
   499 	  break;
       
   500 	}
       
   501 	
       
   502       case GEN_FLOW:
       
   503       case PREFLOW:
       
   504 	{
       
   505 	  //Reverse_bfs from t in the residual graph, 
       
   506 	  //to find the starting level.
       
   507 	  level.set(t,0);
       
   508 	  bfs_queue.push(t);
       
   509 	  
       
   510 	  while (!bfs_queue.empty()) {
       
   511 	    
       
   512 	    Node v=bfs_queue.front();	
       
   513 	    bfs_queue.pop();
       
   514 	    int l=level[v]+1;
       
   515 	    
       
   516 	    InEdgeIt e;
       
   517 	    for(G.first(e,v); G.valid(e); G.next(e)) {
       
   518 	      if ( (*capacity)[e] == (*flow)[e] ) continue;
       
   519 	      Node w=G.tail(e);
       
   520 	      if ( level[w] == n && w != s ) {
       
   521 		bfs_queue.push(w);
       
   522 		Node first=level_list[l];
       
   523 		if ( G.valid(first) ) left.set(first,w);
       
   524 		right.set(w,first);
       
   525 		level_list[l]=w;
       
   526 		level.set(w, l);
       
   527 	      }
       
   528 	    }
       
   529 	    
       
   530 	    OutEdgeIt f;
       
   531 	    for(G.first(f,v); G.valid(f); G.next(f)) {
       
   532 	      if ( 0 == (*flow)[f] ) continue;
       
   533 	      Node w=G.head(f);
       
   534 	      if ( level[w] == n && w != s ) {
       
   535 		bfs_queue.push(w);
       
   536 		Node first=level_list[l];
       
   537 		if ( G.valid(first) ) left.set(first,w);
       
   538 		right.set(w,first);
       
   539 		level_list[l]=w;
       
   540 		level.set(w, l);
       
   541 	      }
       
   542 	    }
       
   543 	  }
       
   544 	  
       
   545 	  
       
   546 	  //the starting flow
       
   547 	  OutEdgeIt e;
       
   548 	  for(G.first(e,s); G.valid(e); G.next(e)) 
       
   549 	    {
       
   550 	      T rem=(*capacity)[e]-(*flow)[e];
       
   551 	      if ( rem == 0 ) continue;
       
   552 	      Node w=G.head(e);
       
   553 	      if ( level[w] < n ) {	  
       
   554 		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
       
   555 		flow->set(e, (*capacity)[e]); 
       
   556 		excess.set(w, excess[w]+rem);
       
   557 	      }
       
   558 	    }
       
   559 	  
       
   560 	  InEdgeIt f;
       
   561 	  for(G.first(f,s); G.valid(f); G.next(f)) 
       
   562 	    {
       
   563 	      if ( (*flow)[f] == 0 ) continue;
       
   564 	      Node w=G.tail(f);
       
   565 	      if ( level[w] < n ) {	  
       
   566 		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
       
   567 		excess.set(w, excess[w]+(*flow)[f]);
       
   568 		flow->set(f, 0); 
       
   569 	      }
       
   570 	    }  
       
   571 	  break;
       
   572 	} //case PREFLOW
       
   573       }
       
   574     } //preflowPreproc
       
   575 
       
   576 
       
   577 
       
   578     void relabel( const Node w, int newlevel, VecStack& active,  
       
   579 		  VecNode& level_list, NNMap& left, 
       
   580 		  NNMap& right, int& b, int& k, const bool what_heur ) {
       
   581 
       
   582       T lev=level[w];	
       
   583       
       
   584       Node right_n=right[w];
       
   585       Node left_n=left[w];
       
   586       
       
   587       //unlacing starts
       
   588       if ( G.valid(right_n) ) {
       
   589 	if ( G.valid(left_n) ) {
       
   590 	  right.set(left_n, right_n);
       
   591 	  left.set(right_n, left_n);
       
   592 	} else {
       
   593 	  level_list[lev]=right_n;   
       
   594 	  left.set(right_n, INVALID);
       
   595 	} 
       
   596       } else {
       
   597 	if ( G.valid(left_n) ) {
       
   598 	  right.set(left_n, INVALID);
       
   599 	} else { 
       
   600 	  level_list[lev]=INVALID;   
       
   601 	} 
       
   602       } 
       
   603       //unlacing ends
       
   604 		
       
   605       if ( !G.valid(level_list[lev]) ) {
       
   606 	      
       
   607 	//gapping starts
       
   608 	for (int i=lev; i!=k ; ) {
       
   609 	  Node v=level_list[++i];
       
   610 	  while ( G.valid(v) ) {
       
   611 	    level.set(v,n);
       
   612 	    v=right[v];
       
   613 	  }
       
   614 	  level_list[i]=INVALID;
       
   615 	  if ( !what_heur ) {
       
   616 	    while ( !active[i].empty() ) {
       
   617 	      active[i].pop();    //FIXME: ezt szebben kene
       
   618 	    }
       
   619 	  }	     
       
   620 	}
       
   621 	
       
   622 	level.set(w,n);
       
   623 	b=lev-1;
       
   624 	k=b;
       
   625 	//gapping ends
       
   626 	
       
   627       } else {
       
   628 	
       
   629 	if ( newlevel == n ) level.set(w,n); 
       
   630 	else {
       
   631 	  level.set(w,++newlevel);
       
   632 	  active[newlevel].push(w);
       
   633 	  if ( what_heur ) b=newlevel;
       
   634 	  if ( k < newlevel ) ++k;      //now k=newlevel
       
   635 	  Node first=level_list[newlevel];
       
   636 	  if ( G.valid(first) ) left.set(first,w);
       
   637 	  right.set(w,first);
       
   638 	  left.set(w,INVALID);
       
   639 	  level_list[newlevel]=w;
       
   640 	}
       
   641       }
       
   642       
       
   643     } //relabel
       
   644     
   655 
   645 
   656   };
   646   };
   657 
   647 
   658 } //namespace hugo
   648 } //namespace hugo
   659 
   649 
   660 #endif //PREFLOW_H
   650 #endif //HUGO_PREFLOW_H
   661 
   651 
   662 
   652 
   663 
   653 
   664 
   654