src/work/jacint/preflow_hl4.h
changeset 138 c6297c121409
parent 105 a3c73e9b9b2e
equal deleted inserted replaced
1:723ab01137f3 2:435e37fdebf6
     1 // -*- C++ -*-
     1 // -*- C++ -*-
     2 /*
     2 /*
     3 preflow_hl4.h
     3 preflow_h5.h
     4 by jacint. 
     4 by jacint. 
     5 Runs the two phase highest label preflow push algorithm. In phase 0
     5 Heuristics: 
     6 we maintain in a list the nodes in level i < n, and we maintain a 
     6  2 phase
     7 bound k on the max level i < n containing a node, so we can do
     7  gap
     8 the gap heuristic fast. Phase 1 is the same. (The algorithm is the 
     8  list 'level_list' on the nodes on level i implemented by hand
     9 same as preflow.hl3, the only diff is that here we use the gap
     9  highest label
    10 heuristic with the list of the nodes on level i, and not a bfs form the
    10  relevel: in phase 0, after BFS*n relabels, it runs a reverse 
    11 upgraded node.)
    11    bfs from t in the res graph to relevel the nodes reachable from t.
    12 
    12    BFS is initialized to 20
    13 In phase 1 we shift everything downwards by n.
    13 
    14 
    14 Due to the last heuristic, this algorithm is quite fast on very 
    15 Member functions:
    15 sparse graphs, but relatively bad on even the dense graphs.
    16 
    16 
    17 void run() : runs the algorithm
    17 'NodeMap<bool> cut' is a member, in this way we can count it fast, after
    18 
    18 the algorithm was run.
    19  The following functions should be used after run() was already run.
    19 
    20 
    20 The constructor runs the algorithm.
    21 T maxflow() : returns the value of a maximum flow
    21 
    22 
    22 Members:
    23 T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e) 
    23 
    24 
    24 T maxFlow() : returns the value of a maximum flow
    25 FlowMap allflow() : returns a maximum flow
    25 
    26 
    26 T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e) 
    27 void allflow(FlowMap& _flow ) : returns a maximum flow
    27 
    28 
    28 FlowMap Flow() : returns the fixed maximum flow x
    29 void mincut(CutMap& M) : sets M to the characteristic vector of a 
    29 
    30      minimum cut. M should be a map of bools initialized to false.
    30 void Flow(FlowMap& _flow ) : returns the fixed maximum flow x
    31 
    31 
    32 void min_mincut(CutMap& M) : sets M to the characteristic vector of the 
    32 void minMinCut(CutMap& M) : sets M to the characteristic vector of the 
    33      minimum min cut. M should be a map of bools initialized to false.
    33      minimum min cut. M should be a map of bools initialized to false.
    34 
    34 
    35 void max_mincut(CutMap& M) : sets M to the characteristic vector of the 
    35 void maxMinCut(CutMap& M) : sets M to the characteristic vector of the 
    36      maximum min cut. M should be a map of bools initialized to false.
    36      maximum min cut. M should be a map of bools initialized to false.
       
    37 
       
    38 void minCut(CutMap& M) : fast function, sets M to the characteristic 
       
    39      vector of a minimum cut. 
       
    40 
       
    41 Different member from the other preflow_hl-s (here we have a member 
       
    42 'NodeMap<bool> cut').
       
    43 
       
    44 CutMap minCut() : fast function, giving the characteristic 
       
    45      vector of a minimum cut.
       
    46 
    37 
    47 
    38 */
    48 */
    39 
    49 
    40 #ifndef PREFLOW_HL4_H
    50 #ifndef PREFLOW_HL4_H
    41 #define PREFLOW_HL4_H
    51 #define PREFLOW_HL4_H
    42 
    52 
       
    53 #define BFS 20
       
    54 
    43 #include <vector>
    55 #include <vector>
    44 #include <stack>
       
    45 #include <queue>
    56 #include <queue>
       
    57 
       
    58 #include <time_measure.h>  //for test
    46 
    59 
    47 namespace hugo {
    60 namespace hugo {
    48 
    61 
    49   template <typename Graph, typename T, 
    62   template <typename Graph, typename T, 
    50     typename FlowMap=typename Graph::EdgeMap<T>, 
    63     typename FlowMap=typename Graph::EdgeMap<T>, 
       
    64     typename CutMap=typename Graph::NodeMap<bool>,
    51     typename CapMap=typename Graph::EdgeMap<T> >
    65     typename CapMap=typename Graph::EdgeMap<T> >
    52   class preflow_hl4 {
    66   class preflow_hl4 {
    53     
    67     
    54     typedef typename Graph::NodeIt NodeIt;
    68     typedef typename Graph::NodeIt NodeIt;
    55     typedef typename Graph::EdgeIt EdgeIt;
    69     typedef typename Graph::EdgeIt EdgeIt;
    60     Graph& G;
    74     Graph& G;
    61     NodeIt s;
    75     NodeIt s;
    62     NodeIt t;
    76     NodeIt t;
    63     FlowMap flow;
    77     FlowMap flow;
    64     CapMap& capacity;  
    78     CapMap& capacity;  
       
    79     CutMap cut;
    65     T value;
    80     T value;
    66     
    81     
    67   public:
    82   public:
    68 
    83 
       
    84     double time;
       
    85 
    69     preflow_hl4(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
    86     preflow_hl4(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
    70       G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
    87       G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), 
    71 
    88       cut(G, false) {
    72 
    89 
    73     void run() {
    90       bool phase=0;   //phase 0 is the 1st phase, phase 1 is the 2nd
    74  
       
    75       bool phase=0;
       
    76       int n=G.nodeNum(); 
    91       int n=G.nodeNum(); 
       
    92       int relabel=0;
       
    93       int heur=(int)BFS*n;
    77       int k=n-2;
    94       int k=n-2;
    78       int b=k;
    95       int b=k;
    79       /*
    96       /*
    80 	b is a bound on the highest level of the stack. 
    97 	b is a bound on the highest level of the stack. 
    81 	k is a bound on the highest nonempty level i < n.
    98 	k is a bound on the highest nonempty level i < n.
    82       */
    99       */
    83 
   100 
    84       typename Graph::NodeMap<int> level(G,n);      
   101       typename Graph::NodeMap<int> level(G,n);      
    85       typename Graph::NodeMap<T> excess(G); 
   102       typename Graph::NodeMap<T> excess(G); 
    86       std::vector<std::stack<NodeIt> > stack(n);    
   103       
       
   104       std::vector<NodeIt> active(n);
       
   105       typename Graph::NodeMap<NodeIt> next(G);
    87       //Stack of the active nodes in level i < n.
   106       //Stack of the active nodes in level i < n.
    88       //We use it in both phases.
   107       //We use it in both phases.
    89 
   108 
    90       typename Graph::NodeMap<NodeIt> left(G);
   109       typename Graph::NodeMap<NodeIt> left(G);
    91       typename Graph::NodeMap<NodeIt> right(G);
   110       typename Graph::NodeMap<NodeIt> right(G);
   105 	bfs_queue.pop();
   124 	bfs_queue.pop();
   106 	int l=level.get(v)+1;
   125 	int l=level.get(v)+1;
   107 
   126 
   108 	for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
   127 	for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
   109 	  NodeIt w=G.tail(e);
   128 	  NodeIt w=G.tail(e);
   110 	  if ( level.get(w) == n ) {
   129 	  if ( level.get(w) == n && w !=s ) {
   111 	    bfs_queue.push(w);
   130 	    bfs_queue.push(w);
   112 	    NodeIt first=level_list[l];
   131 	    NodeIt first=level_list[l];
   113 	    if ( first != 0 ) left.set(first,w);
   132 	    if ( first != 0 ) left.set(first,w);
   114 	    right.set(w,first);
   133 	    right.set(w,first);
   115 	    level_list[l]=w;
   134 	    level_list[l]=w;
   126 	{
   145 	{
   127 	  T c=capacity.get(e);
   146 	  T c=capacity.get(e);
   128 	  if ( c == 0 ) continue;
   147 	  if ( c == 0 ) continue;
   129 	  NodeIt w=G.head(e);
   148 	  NodeIt w=G.head(e);
   130 	  if ( level.get(w) < n ) {	  
   149 	  if ( level.get(w) < n ) {	  
   131 	    if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); 
   150 	    if ( excess.get(w) == 0 && w!=t ) {
       
   151 	      next.set(w,active[level.get(w)]);
       
   152 	      active[level.get(w)]=w;
       
   153 	    }
   132 	    flow.set(e, c); 
   154 	    flow.set(e, c); 
   133 	    excess.set(w, excess.get(w)+c);
   155 	    excess.set(w, excess.get(w)+c);
   134 	  }
   156 	  }
   135 	}
   157 	}
   136       /* 
   158       /* 
   149 	  /*
   171 	  /*
   150 	    In the end of phase 0 we apply a bfs from s in
   172 	    In the end of phase 0 we apply a bfs from s in
   151 	    the residual graph.
   173 	    the residual graph.
   152 	  */
   174 	  */
   153 	  phase=1;
   175 	  phase=1;
       
   176 	  
       
   177 	  //Now have a min cut.
       
   178 	  for( EachNodeIt v=G.template first<EachNodeIt>(); 
       
   179 	       v.valid(); ++v) 
       
   180 	    if (level.get(v) >= n ) cut.set(v,true);
       
   181 	  
       
   182 	  time=currTime();
   154 	  level.set(s,0);
   183 	  level.set(s,0);
   155 	  std::queue<NodeIt> bfs_queue;
   184 	  std::queue<NodeIt> bfs_queue;
   156 	  bfs_queue.push(s);
   185 	  bfs_queue.push(s);
   157 	  
   186 	  
   158 	  while (!bfs_queue.empty()) {
   187 	  while (!bfs_queue.empty()) {
   165 	      if ( capacity.get(e) == flow.get(e) ) continue;
   194 	      if ( capacity.get(e) == flow.get(e) ) continue;
   166 	      NodeIt u=G.tail(e);
   195 	      NodeIt u=G.tail(e);
   167 	      if ( level.get(u) >= n ) { 
   196 	      if ( level.get(u) >= n ) { 
   168 		bfs_queue.push(u);
   197 		bfs_queue.push(u);
   169 		level.set(u, l);
   198 		level.set(u, l);
   170 		if ( excess.get(u) > 0 ) stack[l].push(u);
   199 		if ( excess.get(u) > 0 ) {
       
   200 		    next.set(u,active[l]);
       
   201 		    active[l]=u;
       
   202 		}
   171 	      }
   203 	      }
   172 	    }
   204 	    }
   173 
   205 
   174 	    for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
   206 	    for(OutEdgeIt e=G.template first<OutEdgeIt>(v); e.valid(); ++e) {
   175 	      if ( 0 == flow.get(e) ) continue;
   207 	      if ( 0 == flow.get(e) ) continue;
   176 	      NodeIt u=G.head(e);
   208 	      NodeIt u=G.head(e);
   177 	      if ( level.get(u) >= n ) { 
   209 	      if ( level.get(u) >= n ) { 
   178 		bfs_queue.push(u);
   210 		bfs_queue.push(u);
   179 		level.set(u, l);
   211 		level.set(u, l);
   180 		if ( excess.get(u) > 0 ) stack[l].push(u);
   212 		if ( excess.get(u) > 0 ) {
       
   213 		    next.set(u,active[l]);
       
   214 		    active[l]=u;
       
   215 		}
   181 	      }
   216 	      }
   182 	    }
   217 	    }
   183 	  }
   218 	  }
   184 	  b=n-2;
   219 	  b=n-2;
   185 	}
   220 	}
   186 
   221 
   187 
   222 
   188 	if ( stack[b].empty() ) --b;
   223 	if ( active[b] == 0 ) --b;
   189 	else {
   224 	else {
   190 	  
   225 	  
   191 	  NodeIt w=stack[b].top();        //w is a highest label active node.
   226 	  NodeIt w=active[b];
   192 	  stack[b].pop();           
   227 	  active[b]=next.get(w);
   193 	  int lev=level.get(w);
   228 	  int lev=level.get(w);
   194 	  T exc=excess.get(w);
   229 	  T exc=excess.get(w);
   195 	  int newlevel=n;          //In newlevel we bound the next level of w.
   230 	  int newlevel=n;          //bound on the next level of w.
   196 	  
   231 	  
   197 	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
   232 	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
   198 	    
   233 	    
   199 	    if ( flow.get(e) == capacity.get(e) ) continue; 
   234 	    if ( flow.get(e) == capacity.get(e) ) continue; 
   200 	    NodeIt v=G.head(e);            
   235 	    NodeIt v=G.head(e);            
   201 	    //e=wv	    
   236 	    //e=wv	    
   202 	    
   237 	    
   203 	    if( lev > level.get(v) ) {      
   238 	    if( lev > level.get(v) ) {      
   204 	      /*Push is allowed now*/
   239 	      /*Push is allowed now*/
   205 	      
   240 	      
   206 	      if ( excess.get(v)==0 && v!=t && v!=s ) 
   241 	      if ( excess.get(v)==0 && v!=t && v!=s )  {
   207 		stack[level.get(v)].push(v); 
   242 		int lev_v=level.get(v);
   208 	      /*v becomes active.*/
   243 		next.set(v,active[lev_v]);
       
   244 		active[lev_v]=v;
       
   245 	      }
   209 	      
   246 	      
   210 	      T cap=capacity.get(e);
   247 	      T cap=capacity.get(e);
   211 	      T flo=flow.get(e);
   248 	      T flo=flow.get(e);
   212 	      T remcap=cap-flo;
   249 	      T remcap=cap-flo;
   213 	      
   250 	      
   241 	    //e=vw
   278 	    //e=vw
   242 	    
   279 	    
   243 	    if( lev > level.get(v) ) {  
   280 	    if( lev > level.get(v) ) {  
   244 	      /*Push is allowed now*/
   281 	      /*Push is allowed now*/
   245 	      
   282 	      
   246 	      if ( excess.get(v)==0 && v!=t && v!=s ) 
   283 	      if ( excess.get(v)==0 && v!=t && v!=s ) {
   247 		stack[level.get(v)].push(v); 
   284 		int lev_v=level.get(v);
   248 	      /*v becomes active.*/
   285 		next.set(v,active[lev_v]);
       
   286 		active[lev_v]=v;
       
   287 	      }
   249 	      
   288 	      
   250 	      T flo=flow.get(e);
   289 	      T flo=flow.get(e);
   251 	      
   290 	      
   252 	      if ( flo >= exc ) { 
   291 	      if ( flo >= exc ) { 
   253 		/*A nonsaturating push.*/
   292 		/*A nonsaturating push.*/
   279 	if ( exc > 0 ) {
   318 	if ( exc > 0 ) {
   280 	  //now 'lev' is the old level of w
   319 	  //now 'lev' is the old level of w
   281 	
   320 	
   282 	  if ( phase ) {
   321 	  if ( phase ) {
   283 	    level.set(w,++newlevel);
   322 	    level.set(w,++newlevel);
   284 	    stack[newlevel].push(w);
   323 	    next.set(w,active[newlevel]);
       
   324 	    active[newlevel]=w;
   285 	    b=newlevel;
   325 	    b=newlevel;
   286 	  } else {
   326 	  } else {
   287 	    //unlacing
   327 	    //unlacing
   288 	    NodeIt right_n=right.get(w);
   328 	    NodeIt right_n=right.get(w);
   289 	    NodeIt left_n=left.get(w);
   329 	    NodeIt left_n=left.get(w);
   301 		right.set(left_n, 0);
   341 		right.set(left_n, 0);
   302 	      } else { 
   342 	      } else { 
   303 		level_list[lev]=0;
   343 		level_list[lev]=0;
   304 	      } 
   344 	      } 
   305 	    }
   345 	    }
       
   346 	    //unlacing ends
   306 		
   347 		
   307 
   348 
   308 	    if ( level_list[lev]==0 ) {
   349 	    if ( level_list[lev]==0 ) {
   309 
   350 
   310 	      for (int i=lev; i!=k ; ) {
   351 	      for (int i=lev; i!=k ; ) {
   326 	      if ( newlevel == n ) {
   367 	      if ( newlevel == n ) {
   327 		level.set(w,n);
   368 		level.set(w,n);
   328 	      } else {
   369 	      } else {
   329 		
   370 		
   330 		level.set(w,++newlevel);
   371 		level.set(w,++newlevel);
   331 		stack[newlevel].push(w);
   372 		next.set(w,active[newlevel]);
       
   373 		active[newlevel]=w;
   332 		b=newlevel;
   374 		b=newlevel;
   333 		if ( k < newlevel ) ++k;
   375 		if ( k < newlevel ) ++k;
   334 		NodeIt first=level_list[newlevel];
   376 		NodeIt first=level_list[newlevel];
   335 		if ( first != 0 ) left.set(first,w);
   377 		if ( first != 0 ) left.set(first,w);
   336 		right.set(w,first);
   378 		right.set(w,first);
   337 		left.set(w,0);
   379 		left.set(w,0);
   338 		level_list[newlevel]=w;
   380 		level_list[newlevel]=w;
   339 	      }
   381 	      }
   340 	    }
   382 	    }
       
   383 
       
   384 	    ++relabel;
       
   385 	    if ( relabel >= heur ) {
       
   386 	      relabel=0;
       
   387 	      b=n-2;
       
   388 	      k=b;
       
   389 		
       
   390 	      for ( int i=1; i!=n; ++i ) { 
       
   391 		active[i]=0;
       
   392 		level_list[i]=0;
       
   393 	      }
       
   394 
       
   395 	      //bfs from t in the res graph to relevel the nodes
       
   396 	      for( EachNodeIt v=G.template first<EachNodeIt>(); 
       
   397 		   v.valid(); ++v) level.set(v,n);
       
   398 
       
   399 	      level.set(t,0);
       
   400 	      std::queue<NodeIt> bfs_queue;
       
   401 	      bfs_queue.push(t);
       
   402 	      
       
   403 	      while (!bfs_queue.empty()) {
       
   404 		
       
   405 		NodeIt v=bfs_queue.front();	
       
   406 		bfs_queue.pop();
       
   407 		int l=level.get(v)+1;
       
   408 		
       
   409 		for(InEdgeIt e=G.template first<InEdgeIt>(v); 
       
   410 		    e.valid(); ++e) {
       
   411 		  if ( capacity.get(e) == flow.get(e) ) continue;
       
   412 		  NodeIt u=G.tail(e);
       
   413 		  if ( level.get(u) == n ) { 
       
   414 		    bfs_queue.push(u);
       
   415 		    level.set(u, l);
       
   416 		    if ( excess.get(u) > 0 ) {
       
   417 		      next.set(u,active[l]);
       
   418 		      active[l]=u;
       
   419 		    }
       
   420 		    NodeIt first=level_list[l];
       
   421 		    if ( first != 0 ) left.set(first,w);
       
   422 		    right.set(w,first);
       
   423 		    left.set(w,0);
       
   424 		    level_list[l]=w;
       
   425 		  }
       
   426 		}
       
   427 		
       
   428 		
       
   429 		for(OutEdgeIt e=G.template first<OutEdgeIt>(v); 
       
   430 		    e.valid(); ++e) {
       
   431 		  if ( 0 == flow.get(e) ) continue;
       
   432 		  NodeIt u=G.head(e);
       
   433 		  if ( level.get(u) == n ) { 
       
   434 		    bfs_queue.push(u);
       
   435 		    level.set(u, l);
       
   436 		    if ( excess.get(u) > 0 ) {
       
   437 		      next.set(u,active[l]);
       
   438 		      active[l]=u;
       
   439 		    }
       
   440 		    NodeIt first=level_list[l];
       
   441 		    if ( first != 0 ) left.set(first,w);
       
   442 		    right.set(w,first);
       
   443 		    left.set(w,0);
       
   444 		    level_list[l]=w;
       
   445 		  }
       
   446 		}
       
   447 	      }
       
   448 	      
       
   449 	      level.set(s,n);
       
   450 	    }
       
   451 	  
   341 	  } //phase 0
   452 	  } //phase 0
   342 	} // if ( exc > 0 )
   453 	} // if ( exc > 0 )
   343  
   454 	
   344 	
   455 	
   345 	} // if stack[b] is nonempty
   456 	} // if stack[b] is nonempty
   346 	
   457 	
   347       } // while(true)
   458       } // while(true)
   348 
   459 
   359 
   470 
   360     /*
   471     /*
   361       Returns the maximum value of a flow.
   472       Returns the maximum value of a flow.
   362      */
   473      */
   363 
   474 
   364     T maxflow() {
   475     T maxFlow() {
   365       return value;
   476       return value;
   366     }
   477     }
   367 
   478 
   368 
   479 
   369 
   480 
   370     /*
   481     /*
   371       For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). 
   482       For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). 
   372     */
   483     */
   373 
   484 
   374     T flowonedge(EdgeIt e) {
   485     T flowOnEdge(EdgeIt e) {
   375       return flow.get(e);
   486       return flow.get(e);
   376     }
   487     }
   377 
   488 
   378 
   489 
   379 
   490 
   380     FlowMap allflow() {
   491     FlowMap Flow() {
   381       return flow;
   492       return flow;
   382     }
   493     }
   383 
   494 
   384 
   495 
   385     
   496     
   386     void allflow(FlowMap& _flow ) {
   497     void Flow(FlowMap& _flow ) {
   387       for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
   498       for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
   388 	_flow.set(v,flow.get(v));
   499 	_flow.set(v,flow.get(v));
   389     }
   500     }
   390 
   501 
   391 
   502 
   392 
   503 
   393     /*
   504     /*
   394       Returns the minimum min cut, by a bfs from s in the residual graph.
   505       Returns the minimum min cut, by a bfs from s in the residual graph.
   395     */
   506     */
   396     
   507     
   397     template<typename CutMap>
   508     template<typename _CutMap>
   398     void mincut(CutMap& M) {
   509     void minMinCut(_CutMap& M) {
   399     
   510     
   400       std::queue<NodeIt> queue;
   511       std::queue<NodeIt> queue;
   401       
   512       
   402       M.set(s,true);      
   513       M.set(s,true);      
   403       queue.push(s);
   514       queue.push(s);
   431     /*
   542     /*
   432       Returns the maximum min cut, by a reverse bfs 
   543       Returns the maximum min cut, by a reverse bfs 
   433       from t in the residual graph.
   544       from t in the residual graph.
   434     */
   545     */
   435     
   546     
   436     template<typename CutMap>
   547     template<typename _CutMap>
   437     void max_mincut(CutMap& M) {
   548     void maxMinCut(_CutMap& M) {
   438     
   549     
   439       std::queue<NodeIt> queue;
   550       std::queue<NodeIt> queue;
   440       
   551       
   441       M.set(t,true);        
   552       M.set(t,true);        
   442       queue.push(t);
   553       queue.push(t);
   467       }
   578       }
   468 
   579 
   469     }
   580     }
   470 
   581 
   471 
   582 
   472 
   583     template<typename _CutMap>
   473     template<typename CutMap>
   584     void minCut(_CutMap& M) {
   474     void min_mincut(CutMap& M) {
   585       for( EachNodeIt v=G.template first<EachNodeIt>(); 
   475       mincut(M);
   586 	   v.valid(); ++v) 
   476     }
   587 	M.set(v, cut.get(v));
   477 
   588     }
       
   589 
       
   590    
       
   591     CutMap minCut() {
       
   592       return cut;
       
   593     }
   478 
   594 
   479 
   595 
   480   };
   596   };
   481 }//namespace hugo
   597 }//namespace marci
   482 #endif 
   598 #endif 
   483 
   599 
   484 
   600 
   485 
   601 
   486 
   602