Testing preprocess.
authorjacint
Tue, 27 Apr 2004 10:27:34 +0000
changeset 4379853b743d830
parent 436 6d632cb56ea3
child 438 a0a2709cf178
Testing preprocess.
src/work/jacint/preflow_excess.h
src/work/jacint/preflow_excess_test.cc
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/work/jacint/preflow_excess.h	Tue Apr 27 10:27:34 2004 +0000
     1.3 @@ -0,0 +1,664 @@
     1.4 +// -*- C++ -*-
     1.5 +
     1.6 +//run gyorsan tudna adni a minmincutot a 2 fazis elejen , ne vegyuk be konstruktorba egy cutmapet?
     1.7 +//constzero jo igy?
     1.8 +
     1.9 +//majd marci megmondja betegyem-e bfs-t meg resgraphot
    1.10 +
    1.11 +//constzero helyett az kell hogy flow-e vagy csak preflow, ha flow akor csak
    1.12 +//excess[t]-t kell szmaolni
    1.13 +
    1.14 +/*
    1.15 +Heuristics: 
    1.16 + 2 phase
    1.17 + gap
    1.18 + list 'level_list' on the nodes on level i implemented by hand
    1.19 + stack 'active' on the active nodes on level i implemented by hand
    1.20 + runs heuristic 'highest label' for H1*n relabels
    1.21 + runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
    1.22 + 
    1.23 +Parameters H0 and H1 are initialized to 20 and 10.
    1.24 +
    1.25 +Constructors:
    1.26 +
    1.27 +Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if 
    1.28 +     FlowMap is not constant zero, and should be true if it is
    1.29 +
    1.30 +Members:
    1.31 +
    1.32 +void run()
    1.33 +
    1.34 +T flowValue() : returns the value of a maximum flow
    1.35 +
    1.36 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the 
    1.37 +     minimum min cut. M should be a map of bools initialized to false.
    1.38 +
    1.39 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the 
    1.40 +     maximum min cut. M should be a map of bools initialized to false.
    1.41 +
    1.42 +void minCut(CutMap& M) : sets M to the characteristic vector of 
    1.43 +     a min cut. M should be a map of bools initialized to false.
    1.44 +
    1.45 +FIXME reset
    1.46 +
    1.47 +*/
    1.48 +
    1.49 +#ifndef HUGO_PREFLOW_H
    1.50 +#define HUGO_PREFLOW_H
    1.51 +
    1.52 +#define H0 20
    1.53 +#define H1 1
    1.54 +
    1.55 +#include <vector>
    1.56 +#include <queue>
    1.57 +#include <stack>
    1.58 +
    1.59 +namespace hugo {
    1.60 +
    1.61 +  template <typename Graph, typename T, 
    1.62 +	    typename CapMap=typename Graph::template EdgeMap<T>, 
    1.63 +            typename FlowMap=typename Graph::template EdgeMap<T> >
    1.64 +  class Preflow {
    1.65 +    
    1.66 +    typedef typename Graph::Node Node;
    1.67 +    typedef typename Graph::Edge Edge;
    1.68 +    typedef typename Graph::NodeIt NodeIt;
    1.69 +    typedef typename Graph::OutEdgeIt OutEdgeIt;
    1.70 +    typedef typename Graph::InEdgeIt InEdgeIt;
    1.71 +    
    1.72 +    const Graph& G;
    1.73 +    Node s;
    1.74 +    Node t;
    1.75 +    const CapMap& capacity;  
    1.76 +    FlowMap& flow;
    1.77 +    T value;
    1.78 +    bool constzero;
    1.79 +    bool isflow;
    1.80 +
    1.81 +  public:
    1.82 +    Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity, 
    1.83 +	    FlowMap& _flow, bool _constzero, bool _isflow ) :
    1.84 +      G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero), isflow(_isflow) {}
    1.85 +    
    1.86 +    
    1.87 +    void run() {
    1.88 +      
    1.89 +      value=0;                //for the subsequent runs
    1.90 +
    1.91 +      bool phase=0;        //phase 0 is the 1st phase, phase 1 is the 2nd
    1.92 +      int n=G.nodeNum(); 
    1.93 +      int heur0=(int)(H0*n);  //time while running 'bound decrease' 
    1.94 +      int heur1=(int)(H1*n);  //time while running 'highest label'
    1.95 +      int heur=heur1;         //starting time interval (#of relabels)
    1.96 +      bool what_heur=1;       
    1.97 +      /*
    1.98 +	what_heur is 0 in case 'bound decrease' 
    1.99 +	and 1 in case 'highest label'
   1.100 +      */
   1.101 +      bool end=false;     
   1.102 +      /*
   1.103 +	Needed for 'bound decrease', 'true'
   1.104 +	means no active nodes are above bound b.
   1.105 +      */
   1.106 +      int relabel=0;
   1.107 +      int k=n-2;  //bound on the highest level under n containing a node
   1.108 +      int b=k;    //bound on the highest level under n of an active node
   1.109 +      
   1.110 +      typename Graph::template NodeMap<int> level(G,n);      
   1.111 +      typename Graph::template NodeMap<T> excess(G); 
   1.112 +
   1.113 +      std::vector<std::stack<Node> > active(n);
   1.114 +      /*      std::vector<Node> active(n-1,INVALID);
   1.115 +      typename Graph::template NodeMap<Node> next(G,INVALID);
   1.116 +      //Stack of the active nodes in level i < n.
   1.117 +      //We use it in both phases.*/
   1.118 +
   1.119 +      typename Graph::template NodeMap<Node> left(G,INVALID);
   1.120 +      typename Graph::template NodeMap<Node> right(G,INVALID);
   1.121 +      std::vector<Node> level_list(n,INVALID);
   1.122 +      /*
   1.123 +	List of the nodes in level i<n.
   1.124 +      */
   1.125 +
   1.126 +
   1.127 +      if ( constzero ) {
   1.128 +     
   1.129 +	/*Reverse_bfs from t, to find the starting level.*/
   1.130 +	level.set(t,0);
   1.131 +	std::queue<Node> bfs_queue;
   1.132 +	bfs_queue.push(t);
   1.133 +	
   1.134 +	while (!bfs_queue.empty()) {
   1.135 +	  
   1.136 +	  Node v=bfs_queue.front();	
   1.137 +	  bfs_queue.pop();
   1.138 +	  int l=level[v]+1;
   1.139 +	  
   1.140 +	  InEdgeIt e;
   1.141 +	  for(G.first(e,v); G.valid(e); G.next(e)) {
   1.142 +	    Node w=G.tail(e);
   1.143 +	    if ( level[w] == n && w != s ) {
   1.144 +	      bfs_queue.push(w);
   1.145 +	      Node first=level_list[l];
   1.146 +	      if ( G.valid(first) ) left.set(first,w);
   1.147 +	      right.set(w,first);
   1.148 +	      level_list[l]=w;
   1.149 +	      level.set(w, l);
   1.150 +	    }
   1.151 +	  }
   1.152 +	}
   1.153 +
   1.154 +	//the starting flow
   1.155 +	OutEdgeIt e;
   1.156 +	for(G.first(e,s); G.valid(e); G.next(e)) 
   1.157 +	{
   1.158 +	  T c=capacity[e];
   1.159 +	  if ( c == 0 ) continue;
   1.160 +	  Node w=G.head(e);
   1.161 +	  if ( level[w] < n ) {	  
   1.162 +	    if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
   1.163 +	    flow.set(e, c); 
   1.164 +	    excess.set(w, excess[w]+c);
   1.165 +	  }
   1.166 +	}
   1.167 +      }
   1.168 +      else 
   1.169 +      {
   1.170 +	
   1.171 +	/*
   1.172 +	  Reverse_bfs from t in the residual graph, 
   1.173 +	  to find the starting level.
   1.174 +	*/
   1.175 +	level.set(t,0);
   1.176 +	std::queue<Node> bfs_queue;
   1.177 +	bfs_queue.push(t);
   1.178 +	
   1.179 +	while (!bfs_queue.empty()) {
   1.180 +	  
   1.181 +	  Node v=bfs_queue.front();	
   1.182 +	  bfs_queue.pop();
   1.183 +	  int l=level[v]+1;
   1.184 +	  
   1.185 +	  InEdgeIt e;
   1.186 +	  for(G.first(e,v); G.valid(e); G.next(e)) {
   1.187 +	    if ( capacity[e] == flow[e] ) continue;
   1.188 +	    Node w=G.tail(e);
   1.189 +	    if ( level[w] == n && w != s ) {
   1.190 +	      bfs_queue.push(w);
   1.191 +	      Node first=level_list[l];
   1.192 +	      if ( G.valid(first) ) left.set(first,w);
   1.193 +	      right.set(w,first);
   1.194 +	      level_list[l]=w;
   1.195 +	      level.set(w, l);
   1.196 +	    }
   1.197 +	  }
   1.198 +	    
   1.199 +	  OutEdgeIt f;
   1.200 +	  for(G.first(f,v); G.valid(f); G.next(f)) {
   1.201 +	    if ( 0 == flow[f] ) continue;
   1.202 +	    Node w=G.head(f);
   1.203 +	    if ( level[w] == n && w != s ) {
   1.204 +	      bfs_queue.push(w);
   1.205 +	      Node first=level_list[l];
   1.206 +	      if ( G.valid(first) ) left.set(first,w);
   1.207 +	      right.set(w,first);
   1.208 +	      level_list[l]=w;
   1.209 +	      level.set(w, l);
   1.210 +	    }
   1.211 +	  }
   1.212 +	}
   1.213 +      
   1.214 +	
   1.215 +	/*
   1.216 +	  Counting the excess
   1.217 +	*/
   1.218 +
   1.219 +	if ( !isflow ) {
   1.220 +	  NodeIt v;
   1.221 +	  for(G.first(v); G.valid(v); G.next(v)) {
   1.222 +	    T exc=0;
   1.223 +	    
   1.224 +	    InEdgeIt e;
   1.225 +	    for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e];
   1.226 +	    OutEdgeIt f;
   1.227 +	    for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[f];
   1.228 +	    
   1.229 +	    excess.set(v,exc);	  
   1.230 +	    
   1.231 +	    //putting the active nodes into the stack
   1.232 +	    int lev=level[v];
   1.233 +	    if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
   1.234 +	  }
   1.235 +	} else {
   1.236 +	  T exc=0;
   1.237 +	    
   1.238 +	  InEdgeIt e;
   1.239 +	  for(G.first(e,t); G.valid(e); G.next(e)) exc+=flow[e];
   1.240 +	  OutEdgeIt f;
   1.241 +	  for(G.first(f,t); G.valid(f); G.next(f)) exc-=flow[f];
   1.242 +
   1.243 +	  excess.set(t,exc);	  
   1.244 +	}
   1.245 +
   1.246 +
   1.247 +	//the starting flow
   1.248 +	OutEdgeIt e;
   1.249 +	for(G.first(e,s); G.valid(e); G.next(e)) 
   1.250 +	{
   1.251 +	  T rem=capacity[e]-flow[e];
   1.252 +	  if ( rem == 0 ) continue;
   1.253 +	  Node w=G.head(e);
   1.254 +	  if ( level[w] < n ) {	  
   1.255 +	    if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
   1.256 +	    flow.set(e, capacity[e]); 
   1.257 +	    excess.set(w, excess[w]+rem);
   1.258 +	  }
   1.259 +	}
   1.260 +	
   1.261 +	InEdgeIt f;
   1.262 +	for(G.first(f,s); G.valid(f); G.next(f)) 
   1.263 +	{
   1.264 +	  if ( flow[f] == 0 ) continue;
   1.265 +	  Node w=G.tail(f);
   1.266 +	  if ( level[w] < n ) {	  
   1.267 +	    if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
   1.268 +	    excess.set(w, excess[w]+flow[f]);
   1.269 +	    flow.set(f, 0); 
   1.270 +	  }
   1.271 +	}
   1.272 +      }
   1.273 +
   1.274 +
   1.275 +
   1.276 +
   1.277 +      /* 
   1.278 +	 End of preprocessing 
   1.279 +      */
   1.280 +
   1.281 +
   1.282 +
   1.283 +      /*
   1.284 +	Push/relabel on the highest level active nodes.
   1.285 +      */	
   1.286 +      while ( true ) {
   1.287 +	
   1.288 +	if ( b == 0 ) {
   1.289 +	  if ( phase ) break;
   1.290 +	  
   1.291 +	  if ( !what_heur && !end && k > 0 ) {
   1.292 +	    b=k;
   1.293 +	    end=true;
   1.294 +	  } else {
   1.295 +	    phase=1;
   1.296 +	    level.set(s,0);
   1.297 +	    std::queue<Node> bfs_queue;
   1.298 +	    bfs_queue.push(s);
   1.299 +	    
   1.300 +	    while (!bfs_queue.empty()) {
   1.301 +	      
   1.302 +	      Node v=bfs_queue.front();	
   1.303 +	      bfs_queue.pop();
   1.304 +	      int l=level[v]+1;
   1.305 +	      
   1.306 +	      InEdgeIt e;
   1.307 +	      for(G.first(e,v); G.valid(e); G.next(e)) {
   1.308 +		if ( capacity[e] == flow[e] ) continue;
   1.309 +		Node u=G.tail(e);
   1.310 +		if ( level[u] >= n ) { 
   1.311 +		  bfs_queue.push(u);
   1.312 +		  level.set(u, l);
   1.313 +		  if ( excess[u] > 0 ) active[l].push(u);
   1.314 +		}
   1.315 +	      }
   1.316 +	    
   1.317 +	      OutEdgeIt f;
   1.318 +	      for(G.first(f,v); G.valid(f); G.next(f)) {
   1.319 +		if ( 0 == flow[f] ) continue;
   1.320 +		Node u=G.head(f);
   1.321 +		if ( level[u] >= n ) { 
   1.322 +		  bfs_queue.push(u);
   1.323 +		  level.set(u, l);
   1.324 +		  if ( excess[u] > 0 ) active[l].push(u);
   1.325 +		}
   1.326 +	      }
   1.327 +	    }
   1.328 +	    b=n-2;
   1.329 +	    }
   1.330 +	    
   1.331 +	}
   1.332 +	  
   1.333 +
   1.334 +	///	  
   1.335 +	if ( active[b].empty() ) --b; 
   1.336 +	else {
   1.337 +	  end=false;  
   1.338 +
   1.339 +	  Node w=active[b].top();
   1.340 +	  active[b].pop();
   1.341 +	  int lev=level[w];
   1.342 +	  T exc=excess[w];
   1.343 +	  int newlevel=n;       //bound on the next level of w
   1.344 +	  
   1.345 +	  OutEdgeIt e;
   1.346 +	  for(G.first(e,w); G.valid(e); G.next(e)) {
   1.347 +	    
   1.348 +	    if ( flow[e] == capacity[e] ) continue; 
   1.349 +	    Node v=G.head(e);            
   1.350 +	    //e=wv	    
   1.351 +	    
   1.352 +	    if( lev > level[v] ) {      
   1.353 +	      /*Push is allowed now*/
   1.354 +	      
   1.355 +	      if ( excess[v]==0 && v!=t && v!=s ) {
   1.356 +		int lev_v=level[v];
   1.357 +		active[lev_v].push(v);
   1.358 +	      }
   1.359 +	      
   1.360 +	      T cap=capacity[e];
   1.361 +	      T flo=flow[e];
   1.362 +	      T remcap=cap-flo;
   1.363 +	      
   1.364 +	      if ( remcap >= exc ) {       
   1.365 +		/*A nonsaturating push.*/
   1.366 +		
   1.367 +		flow.set(e, flo+exc);
   1.368 +		excess.set(v, excess[v]+exc);
   1.369 +		exc=0;
   1.370 +		break; 
   1.371 +		
   1.372 +	      } else { 
   1.373 +		/*A saturating push.*/
   1.374 +		
   1.375 +		flow.set(e, cap);
   1.376 +		excess.set(v, excess[v]+remcap);
   1.377 +		exc-=remcap;
   1.378 +	      }
   1.379 +	    } else if ( newlevel > level[v] ){
   1.380 +	      newlevel = level[v];
   1.381 +	    }	    
   1.382 +	    
   1.383 +	  } //for out edges wv 
   1.384 +	
   1.385 +	
   1.386 +	if ( exc > 0 ) {	
   1.387 +	  InEdgeIt e;
   1.388 +	  for(G.first(e,w); G.valid(e); G.next(e)) {
   1.389 +	    
   1.390 +	    if( flow[e] == 0 ) continue; 
   1.391 +	    Node v=G.tail(e);  
   1.392 +	    //e=vw
   1.393 +	    
   1.394 +	    if( lev > level[v] ) {  
   1.395 +	      /*Push is allowed now*/
   1.396 +	      
   1.397 +	      if ( excess[v]==0 && v!=t && v!=s ) {
   1.398 +		int lev_v=level[v];
   1.399 +		active[lev_v].push(v);
   1.400 +	      }
   1.401 +	      
   1.402 +	      T flo=flow[e];
   1.403 +	      
   1.404 +	      if ( flo >= exc ) { 
   1.405 +		/*A nonsaturating push.*/
   1.406 +		
   1.407 +		flow.set(e, flo-exc);
   1.408 +		excess.set(v, excess[v]+exc);
   1.409 +		exc=0;
   1.410 +		break; 
   1.411 +	      } else {                                               
   1.412 +		/*A saturating push.*/
   1.413 +		
   1.414 +		excess.set(v, excess[v]+flo);
   1.415 +		exc-=flo;
   1.416 +		flow.set(e,0);
   1.417 +	      }  
   1.418 +	    } else if ( newlevel > level[v] ) {
   1.419 +	      newlevel = level[v];
   1.420 +	    }	    
   1.421 +	  } //for in edges vw
   1.422 +	  
   1.423 +	} // if w still has excess after the out edge for cycle
   1.424 +	
   1.425 +	excess.set(w, exc);
   1.426 +	///	push
   1.427 +
   1.428 + 
   1.429 +	/*
   1.430 +	  Relabel
   1.431 +	*/
   1.432 +	
   1.433 +
   1.434 +	if ( exc > 0 ) {
   1.435 +	  //now 'lev' is the old level of w
   1.436 +	
   1.437 +	  if ( phase ) {
   1.438 +	    level.set(w,++newlevel);
   1.439 +	    active[newlevel].push(w);
   1.440 +	    b=newlevel;
   1.441 +	  } else {
   1.442 +	    //unlacing starts
   1.443 +	    Node right_n=right[w];
   1.444 +	    Node left_n=left[w];
   1.445 +
   1.446 +	    if ( G.valid(right_n) ) {
   1.447 +	      if ( G.valid(left_n) ) {
   1.448 +		right.set(left_n, right_n);
   1.449 +		left.set(right_n, left_n);
   1.450 +	      } else {
   1.451 +		level_list[lev]=right_n;   
   1.452 +		left.set(right_n, INVALID);
   1.453 +	      } 
   1.454 +	    } else {
   1.455 +	      if ( G.valid(left_n) ) {
   1.456 +		right.set(left_n, INVALID);
   1.457 +	      } else { 
   1.458 +		level_list[lev]=INVALID;   
   1.459 +	      } 
   1.460 +	    } 
   1.461 +	    //unlacing ends
   1.462 +		
   1.463 +	    if ( !G.valid(level_list[lev]) ) {
   1.464 +	      
   1.465 +	       //gapping starts
   1.466 +	      for (int i=lev; i!=k ; ) {
   1.467 +		Node v=level_list[++i];
   1.468 +		while ( G.valid(v) ) {
   1.469 +		  level.set(v,n);
   1.470 +		  v=right[v];
   1.471 +		}
   1.472 +		level_list[i]=INVALID;
   1.473 +		if ( !what_heur ) {
   1.474 +		  while ( !active[i].empty() ) {
   1.475 +		    active[i].pop();    //FIXME: ezt szebben kene
   1.476 +		  }
   1.477 +		}	     
   1.478 +	      }
   1.479 +
   1.480 +	      level.set(w,n);
   1.481 +	      b=lev-1;
   1.482 +	      k=b;
   1.483 +	      //gapping ends
   1.484 +	    
   1.485 +	    } else {
   1.486 +	      
   1.487 +	      if ( newlevel == n ) level.set(w,n); 
   1.488 +	      else {
   1.489 +		level.set(w,++newlevel);
   1.490 +		active[newlevel].push(w);
   1.491 +		if ( what_heur ) b=newlevel;
   1.492 +		if ( k < newlevel ) ++k;      //now k=newlevel
   1.493 +		Node first=level_list[newlevel];
   1.494 +		if ( G.valid(first) ) left.set(first,w);
   1.495 +		right.set(w,first);
   1.496 +		left.set(w,INVALID);
   1.497 +		level_list[newlevel]=w;
   1.498 +	      }
   1.499 +	    }
   1.500 +
   1.501 +
   1.502 +	    ++relabel; 
   1.503 +	    if ( relabel >= heur ) {
   1.504 +	      relabel=0;
   1.505 +	      if ( what_heur ) {
   1.506 +		what_heur=0;
   1.507 +		heur=heur0;
   1.508 +		end=false;
   1.509 +	      } else {
   1.510 +		what_heur=1;
   1.511 +		heur=heur1;
   1.512 +		b=k; 
   1.513 +	      }
   1.514 +	    }
   1.515 +	  } //phase 0
   1.516 +	  
   1.517 +	  
   1.518 +	} // if ( exc > 0 )
   1.519 +	  
   1.520 +	
   1.521 +	}  // if stack[b] is nonempty
   1.522 +	
   1.523 +      } // while(true)
   1.524 +
   1.525 +
   1.526 +      value = excess[t];
   1.527 +      /*Max flow value.*/
   1.528 +     
   1.529 +    } //void run()
   1.530 +
   1.531 +
   1.532 +
   1.533 +
   1.534 +
   1.535 +    /*
   1.536 +      Returns the maximum value of a flow.
   1.537 +     */
   1.538 +
   1.539 +    T flowValue() {
   1.540 +      return value;
   1.541 +    }
   1.542 +
   1.543 +
   1.544 +    FlowMap Flow() {
   1.545 +      return flow;
   1.546 +      }
   1.547 +
   1.548 +
   1.549 +    void Flow(FlowMap& _flow ) {
   1.550 +      NodeIt v;
   1.551 +      for(G.first(v) ; G.valid(v); G.next(v))
   1.552 +	_flow.set(v,flow[v]);
   1.553 +    }
   1.554 +
   1.555 +
   1.556 +
   1.557 +    /*
   1.558 +      Returns the minimum min cut, by a bfs from s in the residual graph.
   1.559 +    */
   1.560 +   
   1.561 +    template<typename _CutMap>
   1.562 +    void minMinCut(_CutMap& M) {
   1.563 +    
   1.564 +      std::queue<Node> queue;
   1.565 +      
   1.566 +      M.set(s,true);      
   1.567 +      queue.push(s);
   1.568 +
   1.569 +      while (!queue.empty()) {
   1.570 +        Node w=queue.front();
   1.571 +	queue.pop();
   1.572 +
   1.573 +	OutEdgeIt e;
   1.574 +	for(G.first(e,w) ; G.valid(e); G.next(e)) {
   1.575 +	  Node v=G.head(e);
   1.576 +	  if (!M[v] && flow[e] < capacity[e] ) {
   1.577 +	    queue.push(v);
   1.578 +	    M.set(v, true);
   1.579 +	  }
   1.580 +	} 
   1.581 +
   1.582 +	InEdgeIt f;
   1.583 +	for(G.first(f,w) ; G.valid(f); G.next(f)) {
   1.584 +	  Node v=G.tail(f);
   1.585 +	  if (!M[v] && flow[f] > 0 ) {
   1.586 +	    queue.push(v);
   1.587 +	    M.set(v, true);
   1.588 +	  }
   1.589 +	} 
   1.590 +      }
   1.591 +    }
   1.592 +
   1.593 +
   1.594 +  
   1.595 +    /*
   1.596 +      Returns the maximum min cut, by a reverse bfs 
   1.597 +      from t in the residual graph.
   1.598 +    */
   1.599 +    
   1.600 +    template<typename _CutMap>
   1.601 +    void maxMinCut(_CutMap& M) {
   1.602 +    
   1.603 +      std::queue<Node> queue;
   1.604 +      
   1.605 +      M.set(t,true);        
   1.606 +      queue.push(t);
   1.607 +
   1.608 +      while (!queue.empty()) {
   1.609 +        Node w=queue.front();
   1.610 +	queue.pop();
   1.611 +
   1.612 +
   1.613 +	InEdgeIt e;
   1.614 +	for(G.first(e,w) ; G.valid(e); G.next(e)) {
   1.615 +	  Node v=G.tail(e);
   1.616 +	  if (!M[v] && flow[e] < capacity[e] ) {
   1.617 +	    queue.push(v);
   1.618 +	    M.set(v, true);
   1.619 +	  }
   1.620 +	}
   1.621 +	
   1.622 +	OutEdgeIt f;
   1.623 +	for(G.first(f,w) ; G.valid(f); G.next(f)) {
   1.624 +	  Node v=G.head(f);
   1.625 +	  if (!M[v] && flow[f] > 0 ) {
   1.626 +	    queue.push(v);
   1.627 +	    M.set(v, true);
   1.628 +	  }
   1.629 +	}
   1.630 +      }
   1.631 +
   1.632 +      NodeIt v;
   1.633 +      for(G.first(v) ; G.valid(v); G.next(v)) {
   1.634 +	M.set(v, !M[v]);
   1.635 +      }
   1.636 +
   1.637 +    }
   1.638 +
   1.639 +
   1.640 +
   1.641 +    template<typename CutMap>
   1.642 +    void minCut(CutMap& M) {
   1.643 +      minMinCut(M);
   1.644 +    }
   1.645 +
   1.646 +    
   1.647 +    void resetTarget (Node _t) {t=_t;}
   1.648 +    void resetSource (Node _s) {s=_s;}
   1.649 +   
   1.650 +    void resetCap (CapMap _cap) {capacity=_cap;}
   1.651 +
   1.652 +    void resetFlow (FlowMap _flow, bool _constzero) {
   1.653 +      flow=_flow;
   1.654 +      constzero=_constzero;
   1.655 +    }
   1.656 +
   1.657 +
   1.658 +
   1.659 +  };
   1.660 +
   1.661 +} //namespace hugo
   1.662 +
   1.663 +#endif //PREFLOW_H
   1.664 +
   1.665 +
   1.666 +
   1.667 +
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/src/work/jacint/preflow_excess_test.cc	Tue Apr 27 10:27:34 2004 +0000
     2.3 @@ -0,0 +1,179 @@
     2.4 +/*
     2.5 +The only difference between preflow.h and preflow_res.h is that the latter
     2.6 +uses the ResGraphWrapper, while the first does not. (Bfs is implemented by
     2.7 +hand in both.) This test program runs Preflow and PreflowRes on the same
     2.8 +graph, tests the result of these implementations and writes the running time
     2.9 +of them.  */
    2.10 +#include <iostream>
    2.11 +
    2.12 +#include <smart_graph.h>
    2.13 +#include <dimacs.h>
    2.14 +#include <preflow_excess.h>
    2.15 +#include <time_measure.h>
    2.16 +
    2.17 +using namespace hugo;
    2.18 +
    2.19 +int main(int, char **) {
    2.20 + 
    2.21 +  typedef SmartGraph Graph;
    2.22 +  
    2.23 +  typedef Graph::Node Node;
    2.24 +  typedef Graph::EdgeIt EdgeIt;
    2.25 +
    2.26 +  Graph G;
    2.27 +  Node s, t;
    2.28 +  Graph::EdgeMap<int> cap(G);
    2.29 +  readDimacsMaxFlow(std::cin, G, s, t, cap);
    2.30 +  Timer ts;
    2.31 +  
    2.32 +  std::cout <<
    2.33 +    "\n  Are we slower?"
    2.34 +	    <<std::endl;
    2.35 +  std::cout <<
    2.36 +    "\n  Running preflow.h on a graph with " << 
    2.37 +    G.nodeNum() << " nodes and " << G.edgeNum() << " edges..."
    2.38 +	   << std::endl<<std::endl;
    2.39 +
    2.40 +
    2.41 +
    2.42 +
    2.43 +
    2.44 +
    2.45 +  
    2.46 +  Graph::EdgeMap<int> flow(G,0);
    2.47 +  Preflow<Graph, int> max_flow_test(G, s, t, cap, flow, 0 , 0);
    2.48 +  ts.reset();
    2.49 +  max_flow_test.run();
    2.50 +  std::cout << "Elapsed time from a preflow: " << std::endl 
    2.51 +	    <<ts << std::endl;
    2.52 +  
    2.53 +  Graph::NodeMap<bool> mincut(G);
    2.54 +  max_flow_test.minMinCut(mincut); 
    2.55 +  int min_min_cut_value=0;
    2.56 +  EdgeIt e;
    2.57 +  for(G.first(e); G.valid(e); G.next(e)) {
    2.58 +    if (mincut[G.tail(e)] && !mincut[G.head(e)]) min_min_cut_value+=cap[e];
    2.59 +  }
    2.60 +
    2.61 +  Graph::NodeMap<bool> cut(G);
    2.62 +  max_flow_test.minCut(cut); 
    2.63 +  int min_cut_value=0;
    2.64 +  for(G.first(e); G.valid(e); G.next(e)) {
    2.65 +    if (cut[G.tail(e)] && !cut[G.head(e)]) 
    2.66 +      min_cut_value+=cap[e];
    2.67 +  }
    2.68 +
    2.69 +  Graph::NodeMap<bool> maxcut(G);
    2.70 +  max_flow_test.maxMinCut(maxcut); 
    2.71 +  int max_min_cut_value=0;
    2.72 +  for(G.first(e); G.valid(e); G.next(e)) {
    2.73 +    if (maxcut[G.tail(e)] && !maxcut[G.head(e)]) 
    2.74 +      max_min_cut_value+=cap[e];
    2.75 +      }
    2.76 +
    2.77 +  std::cout << "\n Checking the result: " <<std::endl;  
    2.78 +  std::cout << "Flow value: "<< max_flow_test.flowValue() << std::endl;
    2.79 +  std::cout << "Min cut value: "<< min_cut_value << std::endl;
    2.80 +  std::cout << "Min min cut value: "<< min_min_cut_value << std::endl;
    2.81 +  std::cout << "Max min cut value: "<< max_min_cut_value << 
    2.82 +    std::endl;
    2.83 +
    2.84 +  if ( max_flow_test.flowValue() == min_cut_value &&
    2.85 +       min_cut_value == min_min_cut_value &&
    2.86 +       min_min_cut_value == max_min_cut_value )
    2.87 +    std::cout << "They are equal! " <<std::endl<< std::endl<<"\n";  
    2.88 +
    2.89 +
    2.90 +
    2.91 +
    2.92 +
    2.93 +  
    2.94 +  Graph::EdgeMap<int> flow2(G,0);
    2.95 +  Preflow<Graph, int> max_flow_test2(G, s, t, cap, flow2, 0 , 1);
    2.96 +  ts.reset();
    2.97 +  max_flow_test2.run();
    2.98 +  std::cout << "Elapsed time from a flow: " << std::endl 
    2.99 +	    << ts << std::endl;
   2.100 +  
   2.101 +  Graph::NodeMap<bool> mincut2(G);
   2.102 +  max_flow_test2.minMinCut(mincut2); 
   2.103 +  int min_min_cut_value2=0;
   2.104 +    for(G.first(e); G.valid(e); G.next(e)) {
   2.105 +    if (mincut2[G.tail(e)] && !mincut2[G.head(e)]) min_min_cut_value2+=cap[e];
   2.106 +  }
   2.107 +
   2.108 +  Graph::NodeMap<bool> cut2(G);
   2.109 +  max_flow_test2.minCut(cut2); 
   2.110 +  int min_cut_value2=0;
   2.111 +  for(G.first(e); G.valid(e); G.next(e)) {
   2.112 +    if (cut2[G.tail(e)] && !cut2[G.head(e)]) 
   2.113 +      min_cut_value2+=cap[e];
   2.114 +  }
   2.115 +
   2.116 +  Graph::NodeMap<bool> maxcut2(G);
   2.117 +  max_flow_test2.maxMinCut(maxcut2); 
   2.118 +  int max_min_cut_value2=0;
   2.119 +  for(G.first(e); G.valid(e); G.next(e)) {
   2.120 +    if (maxcut2[G.tail(e)] && !maxcut2[G.head(e)]) 
   2.121 +      max_min_cut_value2+=cap[e];
   2.122 +      }
   2.123 +  
   2.124 +  std::cout << "\n Checking the result: " <<std::endl;  
   2.125 +  std::cout << "Flow value: "<< max_flow_test2.flowValue() << std::endl;
   2.126 +  std::cout << "Min cut value: "<< min_cut_value2 << std::endl;
   2.127 +  std::cout << "Min min cut value: "<< min_min_cut_value2 << std::endl;
   2.128 +  std::cout << "Max min cut value: "<< max_min_cut_value2 << 
   2.129 +    std::endl;  
   2.130 +  if ( max_flow_test.flowValue() == min_cut_value &&
   2.131 +       min_cut_value == min_min_cut_value &&
   2.132 +       min_min_cut_value == max_min_cut_value )
   2.133 +    std::cout << "They are equal! " <<std::endl;  
   2.134 +
   2.135 +
   2.136 +
   2.137 +
   2.138 +
   2.139 +  Graph::EdgeMap<int> flow3(G,0);
   2.140 +  Preflow<Graph, int> max_flow_test3(G, s, t, cap, flow3, 1 , 1);
   2.141 +  ts.reset();
   2.142 +  max_flow_test3.run();
   2.143 +  std::cout << "Elapsed time from a const zero flow: " << std::endl 
   2.144 +	    <<ts << std::endl;
   2.145 +  
   2.146 +  Graph::NodeMap<bool> mincut3(G);
   2.147 +  max_flow_test3.minMinCut(mincut3); 
   2.148 +  int min_min_cut_value3=0;
   2.149 +  for(G.first(e); G.valid(e); G.next(e)) {
   2.150 +    if (mincut3[G.tail(e)] && !mincut3[G.head(e)]) min_min_cut_value3+=cap[e];
   2.151 +  }
   2.152 +
   2.153 +  Graph::NodeMap<bool> cut3(G);
   2.154 +  max_flow_test3.minCut(cut3); 
   2.155 +  int min_cut_value3=0;
   2.156 +  for(G.first(e); G.valid(e); G.next(e)) {
   2.157 +    if (cut3[G.tail(e)] && !cut3[G.head(e)]) 
   2.158 +      min_cut_value3+=cap[e];
   2.159 +  }
   2.160 +
   2.161 +  Graph::NodeMap<bool> maxcut3(G);
   2.162 +  max_flow_test3.maxMinCut(maxcut3); 
   2.163 +  int max_min_cut_value3=0;
   2.164 +  for(G.first(e); G.valid(e); G.next(e)) {
   2.165 +    if (maxcut3[G.tail(e)] && !maxcut3[G.head(e)]) 
   2.166 +      max_min_cut_value3+=cap[e];
   2.167 +      }
   2.168 +
   2.169 +  std::cout << "\n Checking the result: " <<std::endl;  
   2.170 +  std::cout << "Flow value: "<< max_flow_test3.flowValue() << std::endl;
   2.171 +  std::cout << "Min cut value: "<< min_cut_value3 << std::endl;
   2.172 +  std::cout << "Min min cut value: "<< min_min_cut_value3 << std::endl;
   2.173 +  std::cout << "Max min cut value: "<< max_min_cut_value3 << 
   2.174 +    std::endl;
   2.175 +
   2.176 +  if ( max_flow_test3.flowValue() == min_cut_value3 &&
   2.177 +       min_cut_value3 == min_min_cut_value3 &&
   2.178 +       min_min_cut_value3 == max_min_cut_value3 )
   2.179 +    std::cout << "They are equal! " <<std::endl<< std::endl<<"\n";  
   2.180 +  
   2.181 +  return 0;
   2.182 +}