A max flow algorithm
authorjacint
Tue, 20 Jan 2004 21:22:35 +0000
changeset 20bf088f14b87a
parent 19 3151a1026db9
child 21 181b37336b29
A max flow algorithm
src/work/preflow_push_hl.hh
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/work/preflow_push_hl.hh	Tue Jan 20 21:22:35 2004 +0000
     1.3 @@ -0,0 +1,320 @@
     1.4 +/*
     1.5 +preflow_push_hl.hh
     1.6 +by jacint. 
     1.7 +Runs the highest label variant of the preflow push algorithm with 
     1.8 +running time O(n^2\sqrt(m)). 
     1.9 +
    1.10 +Member functions:
    1.11 +
    1.12 +void run() : runs the algorithm
    1.13 +
    1.14 + The following functions should be used after run() was already run.
    1.15 +
    1.16 +T maxflow() : returns the value of a maximum flow
    1.17 +
    1.18 +T flowonedge(edge_iterator e) : for a fixed maximum flow x it returns x(e) 
    1.19 +
    1.20 +edge_property_vector<graph_type, T> allflow() : returns the fixed maximum flow x
    1.21 +
    1.22 +node_property_vector<graph_type, bool> mincut() : returns a 
    1.23 +     characteristic vector of a minimum cut. (An empty level 
    1.24 +     in the algorithm gives a minimum cut.)
    1.25 +*/
    1.26 +
    1.27 +#ifndef PREFLOW_PUSH_HL_HH
    1.28 +#define PREFLOW_PUSH_HL_HH
    1.29 +
    1.30 +#include <algorithm>
    1.31 +#include <vector>
    1.32 +#include <stack>
    1.33 +
    1.34 +#include <marci_graph_traits.hh>
    1.35 +#include <marci_property_vector.hh>
    1.36 +#include <reverse_bfs.hh>
    1.37 +
    1.38 +namespace marci {
    1.39 +
    1.40 +  template <typename graph_type, typename T>
    1.41 +  class preflow_push_hl {
    1.42 +    
    1.43 +    typedef typename graph_traits<graph_type>::node_iterator node_iterator;
    1.44 +    typedef typename graph_traits<graph_type>::edge_iterator edge_iterator;
    1.45 +    typedef typename graph_traits<graph_type>::each_node_iterator each_node_iterator;
    1.46 +    typedef typename graph_traits<graph_type>::out_edge_iterator out_edge_iterator;
    1.47 +    typedef typename graph_traits<graph_type>::in_edge_iterator in_edge_iterator;
    1.48 +    typedef typename graph_traits<graph_type>::each_edge_iterator each_edge_iterator;
    1.49 +    
    1.50 +
    1.51 +    graph_type& G;
    1.52 +    node_iterator s;
    1.53 +    node_iterator t;
    1.54 +    edge_property_vector<graph_type, T> flow;
    1.55 +    edge_property_vector<graph_type, T>& capacity; 
    1.56 +    T value;
    1.57 +    node_property_vector<graph_type, bool> mincutvector;
    1.58 +
    1.59 +   
    1.60 +  public:
    1.61 +
    1.62 +    preflow_push_hl(graph_type& _G, node_iterator _s, node_iterator _t, edge_property_vector<graph_type, T>& _capacity) : G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), mincutvector(_G, true) { }
    1.63 +
    1.64 +
    1.65 +
    1.66 +
    1.67 +    /*
    1.68 +      The run() function runs the highest label preflow-push, 
    1.69 +      running time: O(n^2\sqrt(m))
    1.70 +    */
    1.71 +    void run() {
    1.72 + 
    1.73 +      node_property_vector<graph_type, int> level(G);         //level of node
    1.74 +      node_property_vector<graph_type, T> excess(G);          //excess of node
    1.75 +            
    1.76 +      int n=number_of(G.first_node());                        //number of nodes 
    1.77 +      int b=n; 
    1.78 +      /*b is a bound on the highest level of an active node. In the beginning it is at most n-2.*/
    1.79 +
    1.80 +      std::vector<std::stack<node_iterator> > stack(2*n-1);    //Stack of the active nodes in level i.
    1.81 +
    1.82 +
    1.83 +
    1.84 +
    1.85 +      /*Reverse_bfs from t, to find the starting level.*/
    1.86 +
    1.87 +      reverse_bfs<list_graph> bfs(G, t);
    1.88 +      bfs.run();
    1.89 +      for(each_node_iterator v=G.first_node(); v.is_valid(); ++v) {
    1.90 +	level.put(v, bfs.dist(v)); 
    1.91 +	//std::cout << "the level of " << v << " is " << bfs.dist(v);
    1.92 +      }
    1.93 +
    1.94 +      /*The level of s is fixed to n*/ 
    1.95 +      level.put(s,n);
    1.96 +
    1.97 +
    1.98 +
    1.99 +
   1.100 +
   1.101 +      /* Starting flow. It is everywhere 0 at the moment. */
   1.102 +     
   1.103 +      for(out_edge_iterator i=G.first_out_edge(s); i.is_valid(); ++i) 
   1.104 +	{
   1.105 +	  node_iterator w=G.head(i);
   1.106 +	  flow.put(i, capacity.get(i)); 
   1.107 +	  stack[bfs.dist(w)].push(w); 
   1.108 +	  excess.put(w, capacity.get(i));
   1.109 +	}
   1.110 +
   1.111 +
   1.112 +      /* 
   1.113 +	 End of preprocessing 
   1.114 +      */
   1.115 +
   1.116 +
   1.117 +
   1.118 +      /*
   1.119 +	Push/relabel on the highest level active nodes.
   1.120 +      */
   1.121 +	
   1.122 +      /*While there exists active node.*/
   1.123 +      while (b) { 
   1.124 +
   1.125 +	/*We decrease the bound if there is no active node of level b.*/
   1.126 +	if (stack[b].empty()) {
   1.127 +	  --b;
   1.128 +	} else {
   1.129 +
   1.130 +	  node_iterator w=stack[b].top();    //w is the highest label active node.
   1.131 +	  stack[b].pop();                    //We delete w from the stack.
   1.132 +	
   1.133 +	  int newlevel=2*n-2;                   //In newlevel we maintain the next level of w.
   1.134 +	
   1.135 +	  for(out_edge_iterator e=G.first_out_edge(w); e.is_valid(); ++e) {
   1.136 +	    node_iterator v=G.head(e);
   1.137 +	    /*e is the edge wv.*/
   1.138 +
   1.139 +	    if (flow.get(e)<capacity.get(e)) {              
   1.140 +	      /*e is an edge of the residual graph */
   1.141 +
   1.142 +	      if(level.get(w)==level.get(v)+1) {      
   1.143 +		/*Push is allowed now*/
   1.144 +
   1.145 +		if (capacity.get(e)-flow.get(e) > excess.get(w)) {       
   1.146 +		  /*A nonsaturating push.*/
   1.147 +		  
   1.148 +		  if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v); 
   1.149 +		  /*v becomes active.*/
   1.150 +		  
   1.151 +		  flow.put(e, flow.get(e)+excess.get(w));
   1.152 +		  excess.put(v, excess.get(v)+excess.get(w));
   1.153 +		  excess.put(w,0);
   1.154 +		  //std::cout << w << " " << v <<" elore elen nonsat pump "  << std::endl;
   1.155 +		  break; 
   1.156 +		} else { 
   1.157 +		  /*A saturating push.*/
   1.158 +
   1.159 +		  if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v); 
   1.160 +		  /*v becomes active.*/
   1.161 +
   1.162 +		  excess.put(v, excess.get(v)+capacity.get(e)-flow.get(e));
   1.163 +		  excess.put(w, excess.get(w)-capacity.get(e)+flow.get(e));
   1.164 +		  flow.put(e, capacity.get(e));
   1.165 +		  //std::cout << w<<" " <<v<<" elore elen sat pump "   << std::endl;
   1.166 +		  if (excess.get(w)==0) break;
   1.167 +		  /*If w is not active any more, then we go on to the next node.*/
   1.168 +		  
   1.169 +		} // if (capacity.get(e)-flow.get(e) > excess.get(w))
   1.170 +	      } // if(level.get(w)==level.get(v)+1)
   1.171 +	    
   1.172 +	      else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
   1.173 +	    
   1.174 +	    } //if (flow.get(e)<capacity.get(e))
   1.175 +	 
   1.176 +	  } //for(out_edge_iterator e=G.first_out_edge(w); e.is_valid(); ++e) 
   1.177 +	  
   1.178 +
   1.179 +
   1.180 +	  for(in_edge_iterator e=G.first_in_edge(w); e.is_valid(); ++e) {
   1.181 +	    node_iterator v=G.tail(e);
   1.182 +	    /*e is the edge vw.*/
   1.183 +
   1.184 +	    if (excess.get(w)==0) break;
   1.185 +	    /*It may happen, that w became inactive in the first for cycle.*/		
   1.186 +	    if(flow.get(e)>0) {             
   1.187 +	      /*e is an edge of the residual graph */
   1.188 +
   1.189 +	      if(level.get(w)==level.get(v)+1) {  
   1.190 +		/*Push is allowed now*/
   1.191 +		
   1.192 +		if (flow.get(e) > excess.get(w)) { 
   1.193 +		  /*A nonsaturating push.*/
   1.194 +		  
   1.195 +		  if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v); 
   1.196 +		  /*v becomes active.*/
   1.197 +
   1.198 +		  flow.put(e, flow.get(e)-excess.get(w));
   1.199 +		  excess.put(v, excess.get(v)+excess.get(w));
   1.200 +		  excess.put(w,0);
   1.201 +		  //std::cout << v << " " << w << " vissza elen nonsat pump "     << std::endl;
   1.202 +		  break; 
   1.203 +		} else {                                               
   1.204 +		  /*A saturating push.*/
   1.205 +		  
   1.206 +		  if (excess.get(v)==0 && v != s) stack[level.get(v)].push(v); 
   1.207 +		  /*v becomes active.*/
   1.208 +		  
   1.209 +		  excess.put(v, excess.get(v)+flow.get(e));
   1.210 +		  excess.put(w, excess.get(w)-flow.get(e));
   1.211 +		  flow.put(e,0);
   1.212 +		  //std::cout << v <<" " << w << " vissza elen sat pump "     << std::endl;
   1.213 +		  if (excess.get(w)==0) { break;}
   1.214 +		} //if (flow.get(e) > excess.get(v)) 
   1.215 +	      } //if(level.get(w)==level.get(v)+1)
   1.216 +	      
   1.217 +	      else {newlevel = newlevel < level.get(v) ? newlevel : level.get(v);}
   1.218 +	      
   1.219 +
   1.220 +	    } //if (flow.get(e)>0)
   1.221 +
   1.222 +	  } //for
   1.223 +
   1.224 +
   1.225 +	  if (excess.get(w)>0) {
   1.226 +	    level.put(w,++newlevel);
   1.227 +	    stack[newlevel].push(w);
   1.228 +	    b=newlevel;
   1.229 +	    //std::cout << "The new level of " << w << " is "<< newlevel <<std::endl; 
   1.230 +	  }
   1.231 +
   1.232 +
   1.233 +	} //else
   1.234 +       
   1.235 +      } //while(b)
   1.236 +
   1.237 +      value = excess.get(t);
   1.238 +      /*Max flow value.*/
   1.239 +
   1.240 +
   1.241 +
   1.242 +
   1.243 +    } //void run()
   1.244 +
   1.245 +
   1.246 +
   1.247 +
   1.248 +
   1.249 +    /*
   1.250 +      Returns the maximum value of a flow.
   1.251 +     */
   1.252 +
   1.253 +    T maxflow() {
   1.254 +      return value;
   1.255 +    }
   1.256 +
   1.257 +
   1.258 +
   1.259 +    /*
   1.260 +      For the maximum flow x found by the algorithm, it returns the flow value on edge e, i.e. x(e). 
   1.261 +    */
   1.262 +
   1.263 +    T flowonedge(edge_iterator e) {
   1.264 +      return flow.get(e);
   1.265 +    }
   1.266 +
   1.267 +
   1.268 +
   1.269 +    /*
   1.270 +      Returns the maximum flow x found by the algorithm.
   1.271 +    */
   1.272 +
   1.273 +    edge_property_vector<graph_type, T> allflow() {
   1.274 +      return flow;
   1.275 +    }
   1.276 +
   1.277 +
   1.278 +
   1.279 +    /*
   1.280 +      Returns a minimum cut by using a reverse bfs from t in the residual graph.
   1.281 +    */
   1.282 +    
   1.283 +    node_property_vector<graph_type, bool> mincut() {
   1.284 +    
   1.285 +      std::queue<node_iterator> queue;
   1.286 +      
   1.287 +      mincutvector.put(t,false);      
   1.288 +      queue.push(t);
   1.289 +
   1.290 +      while (!queue.empty()) {
   1.291 +        node_iterator w=queue.front();
   1.292 +	queue.pop();
   1.293 +
   1.294 +	for(in_edge_iterator e=G.first_in_edge(w) ; e.is_valid(); ++e) {
   1.295 +	  node_iterator v=G.tail(e);
   1.296 +	  if (mincutvector.get(v) && flow.get(e) < capacity.get(e) ) {
   1.297 +	    queue.push(v);
   1.298 +	    mincutvector.put(v, false);
   1.299 +	  }
   1.300 +	} // for
   1.301 +
   1.302 +	for(out_edge_iterator e=G.first_out_edge(w) ; e.is_valid(); ++e) {
   1.303 +	  node_iterator v=G.head(e);
   1.304 +	  if (mincutvector.get(v) && flow.get(e) > 0 ) {
   1.305 +	    queue.push(v);
   1.306 +	    mincutvector.put(v, false);
   1.307 +	  }
   1.308 +	} // for
   1.309 +
   1.310 +      }
   1.311 +
   1.312 +      return mincutvector;
   1.313 +    
   1.314 +    }
   1.315 +
   1.316 +
   1.317 +  };
   1.318 +}//namespace marci
   1.319 +#endif 
   1.320 +
   1.321 +
   1.322 +
   1.323 +