another heuristic
authorjacint
Wed, 18 Feb 2004 21:50:45 +0000
changeset 101d2ac583ed195
parent 100 f1de2ab64e1c
child 102 294cb99af985
another heuristic
src/work/jacint/preflow_hl2.h
src/work/jacint/preflow_hl3.h
src/work/jacint/preflow_push_hl.h
     1.1 --- a/src/work/jacint/preflow_hl2.h	Wed Feb 18 17:27:13 2004 +0000
     1.2 +++ b/src/work/jacint/preflow_hl2.h	Wed Feb 18 21:50:45 2004 +0000
     1.3 @@ -80,7 +80,7 @@
     1.4        IntMap level(G,n);      
     1.5        TMap excess(G); 
     1.6        
     1.7 -      std::vector<int> numb(n+1);    
     1.8 +      std::vector<int> numb(n);    
     1.9        /*
    1.10  	The number of nodes on level i < n. It is
    1.11  	initialized to n+1, because of the reverse_bfs-part.
    1.12 @@ -118,13 +118,14 @@
    1.13        /* Starting flow. It is everywhere 0 at the moment. */     
    1.14        for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) 
    1.15  	{
    1.16 -	  if ( capacity.get(e) == 0 ) continue; 
    1.17 +	  T c=capacity.get(e);
    1.18 +	  if ( c == 0 ) continue;
    1.19  	  NodeIt w=G.head(e);
    1.20  	  if ( w!=s ) {	  
    1.21  	    if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); 
    1.22 -	    flow.set(e, capacity.get(e)); 
    1.23 -	    excess.set(w, excess.get(w)+capacity.get(e));
    1.24 -	  } 
    1.25 +	    flow.set(e, c); 
    1.26 +	    excess.set(w, excess.get(w)+c);
    1.27 +	  }
    1.28  	}
    1.29  
    1.30        /* 
    1.31 @@ -155,7 +156,7 @@
    1.32  	  stack[b].pop();           
    1.33  	  int lev=level.get(w);
    1.34  	  int exc=excess.get(w);
    1.35 -	  int newlevel=2*n-2;      //In newlevel we bound the next level of w.
    1.36 +	  int newlevel=2*n;      //In newlevel we bound the next level of w.
    1.37  	  
    1.38  	  //  if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
    1.39  	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/src/work/jacint/preflow_hl3.h	Wed Feb 18 21:50:45 2004 +0000
     2.3 @@ -0,0 +1,472 @@
     2.4 +// -*- C++ -*-
     2.5 +/*
     2.6 +preflow_hl3.h
     2.7 +by jacint. 
     2.8 +Runs the highest label variant of the preflow push algorithm with 
     2.9 +running time O(n^2\sqrt(m)), with the felszippantos 'empty level' 
    2.10 +and with the two-phase heuristic: if there is no active node of
    2.11 +level at most n, then we go into phase 1, do a bfs
    2.12 +from s, and flow the excess back to s.
    2.13 +
    2.14 +In phase 1 we shift everything downwards by n.
    2.15 +
    2.16 +'A' is a parameter for the empty_level heuristic
    2.17 +
    2.18 +Member functions:
    2.19 +
    2.20 +void run() : runs the algorithm
    2.21 +
    2.22 + The following functions should be used after run() was already run.
    2.23 +
    2.24 +T maxflow() : returns the value of a maximum flow
    2.25 +
    2.26 +T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e) 
    2.27 +
    2.28 +FlowMap allflow() : returns the fixed maximum flow x
    2.29 +
    2.30 +void mincut(CutMap& M) : sets M to the characteristic vector of a 
    2.31 +     minimum cut. M should be a map of bools initialized to false.
    2.32 +
    2.33 +void min_mincut(CutMap& M) : sets M to the characteristic vector of the 
    2.34 +     minimum min cut. M should be a map of bools initialized to false.
    2.35 +
    2.36 +void max_mincut(CutMap& M) : sets M to the characteristic vector of the 
    2.37 +     maximum min cut. M should be a map of bools initialized to false.
    2.38 +
    2.39 +*/
    2.40 +
    2.41 +#ifndef PREFLOW_HL3_H
    2.42 +#define PREFLOW_HL3_H
    2.43 +
    2.44 +#define A 1
    2.45 +
    2.46 +#include <vector>
    2.47 +#include <stack>
    2.48 +#include <queue>
    2.49 +
    2.50 +namespace marci {
    2.51 +
    2.52 +  template <typename Graph, typename T, 
    2.53 +    typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>, 
    2.54 +    typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
    2.55 +  class preflow_hl3 {
    2.56 +    
    2.57 +    typedef typename Graph::NodeIt NodeIt;
    2.58 +    typedef typename Graph::EdgeIt EdgeIt;
    2.59 +    typedef typename Graph::EachNodeIt EachNodeIt;
    2.60 +    typedef typename Graph::OutEdgeIt OutEdgeIt;
    2.61 +    typedef typename Graph::InEdgeIt InEdgeIt;
    2.62 +    
    2.63 +    Graph& G;
    2.64 +    NodeIt s;
    2.65 +    NodeIt t;
    2.66 +    FlowMap flow;
    2.67 +    CapMap& capacity;  
    2.68 +    T value;
    2.69 +    
    2.70 +  public:
    2.71 +
    2.72 +    preflow_hl3(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
    2.73 +      G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
    2.74 +
    2.75 +
    2.76 +    void run() {
    2.77 + 
    2.78 +      bool phase=0;
    2.79 +      int n=G.nodeNum(); 
    2.80 +      int b=n-2; 
    2.81 +      /*
    2.82 +	b is a bound on the highest level of the stack. 
    2.83 +	In the beginning it is at most n-2.
    2.84 +      */
    2.85 +
    2.86 +      IntMap level(G,n);      
    2.87 +      TMap excess(G); 
    2.88 +      
    2.89 +      std::vector<int> numb(n);    
    2.90 +      /*
    2.91 +	The number of nodes on level i < n. It is
    2.92 +	initialized to n+1, because of the reverse_bfs-part.
    2.93 +	Needed only in phase 0.
    2.94 +      */
    2.95 +
    2.96 +      std::vector<std::stack<NodeIt> > stack(n);    
    2.97 +      //Stack of the active nodes in level i < n.
    2.98 +      //We use it in both phases.
    2.99 +
   2.100 +
   2.101 +      /*Reverse_bfs from t, to find the starting level.*/
   2.102 +      level.set(t,0);
   2.103 +      std::queue<NodeIt> bfs_queue;
   2.104 +      bfs_queue.push(t);
   2.105 +
   2.106 +      while (!bfs_queue.empty()) {
   2.107 +
   2.108 +	NodeIt v=bfs_queue.front();	
   2.109 +	bfs_queue.pop();
   2.110 +	int l=level.get(v)+1;
   2.111 +
   2.112 +	for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
   2.113 +	  NodeIt w=G.tail(e);
   2.114 +	  if ( level.get(w) == n ) {
   2.115 +	    bfs_queue.push(w);
   2.116 +	    ++numb[l];
   2.117 +	    level.set(w, l);
   2.118 +	  }
   2.119 +	}
   2.120 +      }
   2.121 +      
   2.122 +      level.set(s,n);
   2.123 +
   2.124 +
   2.125 +
   2.126 +      /* Starting flow. It is everywhere 0 at the moment. */     
   2.127 +      for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) 
   2.128 +	{
   2.129 +	  T c=capacity.get(e);
   2.130 +	  if ( c == 0 ) continue;
   2.131 +	  NodeIt w=G.head(e);
   2.132 +	  if ( level.get(w) < n ) {	  
   2.133 +	    if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); 
   2.134 +	    flow.set(e, c); 
   2.135 +	    excess.set(w, excess.get(w)+c);
   2.136 +	  }
   2.137 +	}
   2.138 +
   2.139 +      /* 
   2.140 +	 End of preprocessing 
   2.141 +      */
   2.142 +
   2.143 +
   2.144 +
   2.145 +      /*
   2.146 +	Push/relabel on the highest level active nodes.
   2.147 +      */	
   2.148 +      /*While there exists an active node.*/
   2.149 +      while ( true ) {
   2.150 +
   2.151 +	if ( b == 0 ) {
   2.152 +	  if ( phase ) break; 
   2.153 +	  phase=1;
   2.154 +
   2.155 +	  level.set(s,0);
   2.156 +
   2.157 +	  std::queue<NodeIt> bfs_queue;
   2.158 +	  bfs_queue.push(s);
   2.159 +	  
   2.160 +	  while (!bfs_queue.empty()) {
   2.161 +	    
   2.162 +	    NodeIt v=bfs_queue.front();	
   2.163 +	    bfs_queue.pop();
   2.164 +	    int l=level.get(v)+1;
   2.165 +
   2.166 +	    for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
   2.167 +	      if ( capacity.get(e) == flow.get(e) ) continue;
   2.168 +	      NodeIt u=G.tail(e);
   2.169 +	      if ( level.get(u) == n ) { 
   2.170 +		bfs_queue.push(u);
   2.171 +		level.set(u, l);
   2.172 +		if ( excess.get(u) > 0 ) stack[l].push(u);
   2.173 +	      }
   2.174 +	    }
   2.175 +
   2.176 +	    for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
   2.177 +	      if ( 0 == flow.get(e) ) continue;
   2.178 +	      NodeIt u=G.head(e);
   2.179 +	      if ( level.get(u) == n ) { 
   2.180 +		bfs_queue.push(u);
   2.181 +		level.set(u, l);
   2.182 +		if ( excess.get(u) > 0 ) stack[l].push(u);
   2.183 +	      }
   2.184 +	    }
   2.185 +	  }
   2.186 +
   2.187 +	  b=n-2;
   2.188 +	}
   2.189 +
   2.190 +	if ( stack[b].empty() ) --b;
   2.191 +	else {
   2.192 +	  
   2.193 +	  NodeIt w=stack[b].top();        //w is a highest label active node.
   2.194 +	  stack[b].pop();           
   2.195 +	  int lev=level.get(w);
   2.196 +	  int exc=excess.get(w);
   2.197 +	  int newlevel=n;                 //In newlevel we bound the next level of w.
   2.198 +	  
   2.199 +	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
   2.200 +	    
   2.201 +	    if ( flow.get(e) == capacity.get(e) ) continue; 
   2.202 +	    NodeIt v=G.head(e);            
   2.203 +	    //e=wv	    
   2.204 +	    
   2.205 +	    if( lev > level.get(v) ) {      
   2.206 +	      /*Push is allowed now*/
   2.207 +	      
   2.208 +	      if ( excess.get(v)==0 && v !=t && v!=s ) 
   2.209 +		stack[level.get(v)].push(v); 
   2.210 +	      /*v becomes active.*/
   2.211 +	      
   2.212 +	      int cap=capacity.get(e);
   2.213 +	      int flo=flow.get(e);
   2.214 +	      int remcap=cap-flo;
   2.215 +	      
   2.216 +	      if ( remcap >= exc ) {       
   2.217 +		/*A nonsaturating push.*/
   2.218 +		
   2.219 +		flow.set(e, flo+exc);
   2.220 +		excess.set(v, excess.get(v)+exc);
   2.221 +		exc=0;
   2.222 +		break; 
   2.223 +		
   2.224 +	      } else { 
   2.225 +		/*A saturating push.*/
   2.226 +		
   2.227 +		flow.set(e, cap );
   2.228 +		excess.set(v, excess.get(v)+remcap);
   2.229 +		exc-=remcap;
   2.230 +	      }
   2.231 +	    } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
   2.232 +	    
   2.233 +	  } //for out edges wv 
   2.234 +	
   2.235 +	
   2.236 +	if ( exc > 0 ) {	
   2.237 +	  for( InEdgeIt e=G.template first<InEdgeIt>(w); e.valid(); ++e) {
   2.238 +	    
   2.239 +	    if( flow.get(e) == 0 ) continue; 
   2.240 +	    NodeIt v=G.tail(e);  
   2.241 +	    //e=vw
   2.242 +	    
   2.243 +	    if( lev > level.get(v) ) {  
   2.244 +	      /*Push is allowed now*/
   2.245 +	      
   2.246 +	      if ( excess.get(v)==0 && v !=t && v!=s ) 
   2.247 +		stack[level.get(v)].push(v); 
   2.248 +	      /*v becomes active.*/
   2.249 +	      
   2.250 +	      int flo=flow.get(e);
   2.251 +	      
   2.252 +	      if ( flo >= exc ) { 
   2.253 +		/*A nonsaturating push.*/
   2.254 +		
   2.255 +		flow.set(e, flo-exc);
   2.256 +		excess.set(v, excess.get(v)+exc);
   2.257 +		exc=0;
   2.258 +		break; 
   2.259 +	      } else {                                               
   2.260 +		/*A saturating push.*/
   2.261 +		
   2.262 +		excess.set(v, excess.get(v)+flo);
   2.263 +		exc-=flo;
   2.264 +		flow.set(e,0);
   2.265 +	      }  
   2.266 +	    } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
   2.267 +	    
   2.268 +	  } //for in edges vw
   2.269 +	  
   2.270 +	} // if w still has excess after the out edge for cycle
   2.271 +	
   2.272 +	excess.set(w, exc);
   2.273 +	  
   2.274 +
   2.275 +	/*
   2.276 +	  Relabel
   2.277 +	*/
   2.278 +	
   2.279 +	if ( exc > 0 ) {
   2.280 +	  //now 'lev' is the old level of w
   2.281 +	
   2.282 +	  if ( phase ) {
   2.283 +	    level.set(w,++newlevel);
   2.284 +	    stack[newlevel].push(w);
   2.285 +	    b=newlevel;
   2.286 +	  } else {
   2.287 +
   2.288 +	    if ( newlevel >= n-2 || --numb[lev] == 0 ) {
   2.289 +	      
   2.290 +	      level.set(w,n);
   2.291 +	      
   2.292 +	      if ( newlevel < n ) {
   2.293 +		
   2.294 +		std::queue<NodeIt> bfs_queue;
   2.295 +		bfs_queue.push(w);
   2.296 +
   2.297 +		while (!bfs_queue.empty()) {
   2.298 +
   2.299 +		  NodeIt v=bfs_queue.front();	
   2.300 +		  bfs_queue.pop();
   2.301 +
   2.302 +		  for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
   2.303 +		    if ( capacity.get(e) == flow.get(e) ) continue;
   2.304 +		    NodeIt u=G.head(e);
   2.305 +		    if ( level.get(u) < n ) { 
   2.306 +		      bfs_queue.push(u);
   2.307 +		      --numb[level.get(u)];
   2.308 +		      level.set(u, n);
   2.309 +		    }
   2.310 +		  }
   2.311 +
   2.312 +		  for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
   2.313 +		    if ( 0 == flow.get(e) ) continue;
   2.314 +		    NodeIt u=G.tail(e);
   2.315 +		    if ( level.get(u) < n ) { 
   2.316 +		      bfs_queue.push(u);
   2.317 +		      --numb[level.get(u)];
   2.318 +		      level.set(u, n);
   2.319 +		    }
   2.320 +		  }
   2.321 +		}
   2.322 +	      }
   2.323 +	      b=n-1;
   2.324 +
   2.325 +	    } else {
   2.326 +	      level.set(w,++newlevel);
   2.327 +	      stack[newlevel].push(w);
   2.328 +	      ++numb[newlevel];
   2.329 +	      b=newlevel;
   2.330 +	    }
   2.331 +	  }
   2.332 +	}
   2.333 +
   2.334 + 
   2.335 +	
   2.336 +	} // if stack[b] is nonempty
   2.337 +	
   2.338 +      } // while(true)
   2.339 +
   2.340 +
   2.341 +      value = excess.get(t);
   2.342 +      /*Max flow value.*/
   2.343 +
   2.344 +
   2.345 +    } //void run()
   2.346 +
   2.347 +
   2.348 +
   2.349 +
   2.350 +
   2.351 +    /*
   2.352 +      Returns the maximum value of a flow.
   2.353 +     */
   2.354 +
   2.355 +    T maxflow() {
   2.356 +      return value;
   2.357 +    }
   2.358 +
   2.359 +
   2.360 +
   2.361 +    /*
   2.362 +      For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). 
   2.363 +    */
   2.364 +
   2.365 +    T flowonedge(EdgeIt e) {
   2.366 +      return flow.get(e);
   2.367 +    }
   2.368 +
   2.369 +
   2.370 +
   2.371 +    /*
   2.372 +      Returns the maximum flow x found by the algorithm.
   2.373 +    */
   2.374 +
   2.375 +    FlowMap allflow() {
   2.376 +      return flow;
   2.377 +    }
   2.378 +
   2.379 +
   2.380 +
   2.381 +
   2.382 +    /*
   2.383 +      Returns the minimum min cut, by a bfs from s in the residual graph.
   2.384 +    */
   2.385 +    
   2.386 +    template<typename CutMap>
   2.387 +    void mincut(CutMap& M) {
   2.388 +    
   2.389 +      std::queue<NodeIt> queue;
   2.390 +      
   2.391 +      M.set(s,true);      
   2.392 +      queue.push(s);
   2.393 +
   2.394 +      while (!queue.empty()) {
   2.395 +        NodeIt w=queue.front();
   2.396 +	queue.pop();
   2.397 +
   2.398 +	for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
   2.399 +	  NodeIt v=G.head(e);
   2.400 +	  if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
   2.401 +	    queue.push(v);
   2.402 +	    M.set(v, true);
   2.403 +	  }
   2.404 +	} 
   2.405 +
   2.406 +	for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
   2.407 +	  NodeIt v=G.tail(e);
   2.408 +	  if (!M.get(v) && flow.get(e) > 0 ) {
   2.409 +	    queue.push(v);
   2.410 +	    M.set(v, true);
   2.411 +	  }
   2.412 +	} 
   2.413 +
   2.414 +      }
   2.415 +
   2.416 +    }
   2.417 +
   2.418 +
   2.419 +
   2.420 +    /*
   2.421 +      Returns the maximum min cut, by a reverse bfs 
   2.422 +      from t in the residual graph.
   2.423 +    */
   2.424 +    
   2.425 +    template<typename CutMap>
   2.426 +    void max_mincut(CutMap& M) {
   2.427 +    
   2.428 +      std::queue<NodeIt> queue;
   2.429 +      
   2.430 +      M.set(t,true);        
   2.431 +      queue.push(t);
   2.432 +
   2.433 +      while (!queue.empty()) {
   2.434 +        NodeIt w=queue.front();
   2.435 +	queue.pop();
   2.436 +
   2.437 +	for(InEdgeIt e=G.template first<InEdgeIt>(w) ; e.valid(); ++e) {
   2.438 +	  NodeIt v=G.tail(e);
   2.439 +	  if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
   2.440 +	    queue.push(v);
   2.441 +	    M.set(v, true);
   2.442 +	  }
   2.443 +	}
   2.444 +
   2.445 +	for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
   2.446 +	  NodeIt v=G.head(e);
   2.447 +	  if (!M.get(v) && flow.get(e) > 0 ) {
   2.448 +	    queue.push(v);
   2.449 +	    M.set(v, true);
   2.450 +	  }
   2.451 +	}
   2.452 +      }
   2.453 +
   2.454 +      for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v) {
   2.455 +	M.set(v, !M.get(v));
   2.456 +      }
   2.457 +
   2.458 +    }
   2.459 +
   2.460 +
   2.461 +
   2.462 +    template<typename CutMap>
   2.463 +    void min_mincut(CutMap& M) {
   2.464 +      mincut(M);
   2.465 +    }
   2.466 +
   2.467 +
   2.468 +
   2.469 +  };
   2.470 +}//namespace marci
   2.471 +#endif 
   2.472 +
   2.473 +
   2.474 +
   2.475 +
     3.1 --- a/src/work/jacint/preflow_push_hl.h	Wed Feb 18 17:27:13 2004 +0000
     3.2 +++ b/src/work/jacint/preflow_push_hl.h	Wed Feb 18 21:50:45 2004 +0000
     3.3 @@ -1,4 +1,8 @@
     3.4  // -*- C++ -*-
     3.5 +
     3.6 +//kerdesek: nem tudom lehet-e a 
     3.7 +//kieleket csak a legf n szintu pontokra nezni.
     3.8 +
     3.9  /*
    3.10  preflow_push_hl.h
    3.11  by jacint. 
    3.12 @@ -116,12 +120,13 @@
    3.13        /* Starting flow. It is everywhere 0 at the moment. */     
    3.14        for(OutEdgeIt e=G.template first<OutEdgeIt>(s); e.valid(); ++e) 
    3.15  	{
    3.16 -	  if ( capacity.get(e) == 0 ) continue;
    3.17 +	  T c=capacity.get(e);
    3.18 +	  if ( c == 0 ) continue;
    3.19  	  NodeIt w=G.head(e);
    3.20  	  if ( w!=s ) {	  
    3.21  	    if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); 
    3.22 -	    flow.set(e, capacity.get(e)); 
    3.23 -	    excess.set(w, excess.get(w)+capacity.get(e));
    3.24 +	    flow.set(e, c); 
    3.25 +	    excess.set(w, excess.get(w)+c);
    3.26  	  }
    3.27  	}
    3.28        
    3.29 @@ -144,8 +149,9 @@
    3.30  	stack[b].pop();           
    3.31  	int lev=level.get(w);
    3.32  	int exc=excess.get(w);
    3.33 -	int newlevel=2*n-2;      //In newlevel we bound the next level of w.
    3.34 -	
    3.35 +	int newlevel=2*n;      //In newlevel we bound the next level of w.
    3.36 +	//vagy MAXINT	
    3.37 +
    3.38  	//  if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
    3.39  	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
    3.40