src/work/jacint/preflow_hl2.h
changeset 110 3c53698842dd
parent 105 a3c73e9b9b2e
equal deleted inserted replaced
2:709167d7212b 3:1469a2bdf3e1
     1 // -*- C++ -*-
     1 // -*- C++ -*-
     2 /*
     2 /*
     3 preflow_hl2.h
     3 preflow_hl2.h
     4 by jacint. 
     4 by jacint. 
     5 Runs the highest label variant of the preflow push algorithm with 
     5 Runs the highest label variant of the preflow push algorithm with 
     6 running time O(n^2\sqrt(m)), with the 'empty level' and with the 
     6 running time O(n^2\sqrt(m)). 
     7 heuristic that the bound b on the active nodes is not increased 
     7 
     8 only when b=0, when we put b=2*n-2.
     8 Heuristics:
     9 
     9 
    10 'A' is a parameter for the empty_level heuristic
    10   gap: we iterate through the nodes for finding the nodes above 
    11 
    11        the gap and under level n. So it is quite slow.
    12 Member functions:
    12   numb: we maintain the number of nodes in level i.
    13 
    13   highest label
    14 void run() : runs the algorithm
    14   
    15 
    15 'A' is a parameter for the gap, we only upgrade the nodes to level n,
    16  The following functions should be used after run() was already run.
    16   if the gap is under A*n.
    17 
    17 
    18 T maxflow() : returns the value of a maximum flow
    18 The constructor runs the algorithm.
    19 
    19 
    20 T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e) 
    20 Members:
    21 
    21 
    22 FlowMap allflow() : returns the fixed maximum flow x
    22 T maxFlow() : returns the value of a maximum flow
    23 
    23 
    24 void mincut(CutMap& M) : sets M to the characteristic vector of a 
    24 T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e) 
       
    25 
       
    26 FlowMap Flow() : returns the fixed maximum flow x
       
    27 
       
    28 void minCut(CutMap& M) : sets M to the characteristic vector of a 
    25      minimum cut. M should be a map of bools initialized to false.
    29      minimum cut. M should be a map of bools initialized to false.
    26 
    30 
    27 void min_mincut(CutMap& M) : sets M to the characteristic vector of the 
    31 void minMinCut(CutMap& M) : sets M to the characteristic vector of the 
    28      minimum min cut. M should be a map of bools initialized to false.
    32      minimum min cut. M should be a map of bools initialized to false.
    29 
    33 
    30 void max_mincut(CutMap& M) : sets M to the characteristic vector of the 
    34 void maxMinCut(CutMap& M) : sets M to the characteristic vector of the 
    31      maximum min cut. M should be a map of bools initialized to false.
    35      maximum min cut. M should be a map of bools initialized to false.
    32 
    36 
    33 */
    37 */
    34 
    38 
    35 #ifndef PREFLOW_HL2_H
    39 #ifndef PREFLOW_HL2_H
    36 #define PREFLOW_HL2_H
    40 #define PREFLOW_HL2_H
    37 
    41 
    38 #define A 1
    42 #define A .9
    39 
    43 
    40 #include <vector>
    44 #include <vector>
    41 #include <stack>
    45 #include <stack>
    42 #include <queue>
    46 #include <queue>
    43 
    47 
    44 namespace hugo {
    48 namespace hugo {
    45 
    49 
    46   template <typename Graph, typename T, 
    50   template <typename Graph, typename T, 
    47     typename FlowMap=typename Graph::EdgeMap<T>, typename CapMap=typename Graph::EdgeMap<T>, 
    51     typename FlowMap=typename Graph::EdgeMap<T>, 
    48     typename IntMap=typename Graph::NodeMap<int>, typename TMap=typename Graph::NodeMap<T> >
    52     typename CapMap=typename Graph::EdgeMap<T> >
    49   class preflow_hl2 {
    53   class preflow_hl2 {
    50     
    54     
    51     typedef typename Graph::NodeIt NodeIt;
    55     typedef typename Graph::NodeIt NodeIt;
    52     typedef typename Graph::EdgeIt EdgeIt;
    56     typedef typename Graph::EdgeIt EdgeIt;
    53     typedef typename Graph::EachNodeIt EachNodeIt;
    57     typedef typename Graph::EachNodeIt EachNodeIt;
    62     T value;
    66     T value;
    63     
    67     
    64   public:
    68   public:
    65 
    69 
    66     preflow_hl2(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
    70     preflow_hl2(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
    67       G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
    71       G(_G), s(_s), t(_t), flow(_G), capacity(_capacity) { 
    68 
    72 
    69 
       
    70     void run() {
       
    71  
       
    72       bool no_end=true;
       
    73       int n=G.nodeNum(); 
    73       int n=G.nodeNum(); 
    74       int b=n-2; 
    74       int b=n-2; 
    75       /*
    75       /*
    76 	b is a bound on the highest level of an active node. 
    76 	b is a bound on the highest level of an active node. 
    77 	In the beginning it is at most n-2.
       
    78       */
    77       */
    79 
    78 
    80       IntMap level(G,n);      
    79       typename Graph::NodeMap<int> level(G,n);      
    81       TMap excess(G); 
    80       typename Graph::NodeMap<T> excess(G); 
    82       
    81 
    83       std::vector<int> numb(n);    
    82       std::vector<int> numb(n);    
    84       /*
    83       /*
    85 	The number of nodes on level i < n. It is
    84 	The number of nodes on level i < n. It is
    86 	initialized to n+1, because of the reverse_bfs-part.
    85 	initialized to n+1, because of the reverse_bfs-part.
    87       */
    86       */
   108 	    ++numb[l];
   107 	    ++numb[l];
   109 	    level.set(w, l);
   108 	    level.set(w, l);
   110 	  }
   109 	  }
   111 	}
   110 	}
   112       }
   111       }
   113       
   112 	
   114       level.set(s,n);
   113       level.set(s,n);
   115 
   114 
   116 
   115 
   117 
   116 
   118       /* Starting flow. It is everywhere 0 at the moment. */     
   117       /* Starting flow. It is everywhere 0 at the moment. */     
   125 	    if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); 
   124 	    if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); 
   126 	    flow.set(e, c); 
   125 	    flow.set(e, c); 
   127 	    excess.set(w, excess.get(w)+c);
   126 	    excess.set(w, excess.get(w)+c);
   128 	  }
   127 	  }
   129 	}
   128 	}
   130 
   129       
   131       /* 
   130       /* 
   132 	 End of preprocessing 
   131 	 End of preprocessing 
   133       */
   132       */
   134 
   133 
   135 
   134 
   136 
       
   137       /*
   135       /*
   138 	Push/relabel on the highest level active nodes.
   136 	Push/relabel on the highest level active nodes.
   139       */	
   137       */
   140       /*While there exists an active node.*/
   138       /*While there exists an active node.*/
   141       while (b) { 
   139       while (b) { 
   142 	if ( stack[b].empty() ) {
   140 	if ( stack[b].empty() ) { 
   143 	  if ( b==1 ) {
       
   144 	    if ( !no_end ) break; 
       
   145 	    else {
       
   146 	      b=2*n-2;
       
   147 	      no_end=false;
       
   148 	    }
       
   149 	  } 
       
   150 	  --b;
   141 	  --b;
   151 	} else {
   142 	  continue;
   152 	  
   143 	} 
   153 	  no_end=true;
   144 	
   154 	  
   145 	NodeIt w=stack[b].top();   
   155 	  NodeIt w=stack[b].top();        //w is a highest label active node.
   146 	stack[b].pop();           
   156 	  stack[b].pop();           
   147 	int lev=level.get(w);
   157 	  int lev=level.get(w);
   148 	T exc=excess.get(w);
   158 	  int exc=excess.get(w);
   149 	int newlevel=2*n;      //In newlevel we bound the next level of w.
   159 	  int newlevel=2*n;      //In newlevel we bound the next level of w.
   150 	
   160 	  
       
   161 	  //  if ( level.get(w) < n ) { //Nem tudom ez mukodik-e
       
   162 	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
   151 	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
   163 	    
   152 	    
   164 	    if ( flow.get(e) == capacity.get(e) ) continue; 
   153 	    if ( flow.get(e) == capacity.get(e) ) continue; 
   165 	    NodeIt v=G.head(e);            
   154 	    NodeIt v=G.head(e);            
   166 	    //e=wv	    
   155 	    //e=wv	    
   170 	      
   159 	      
   171 	      if ( excess.get(v)==0 && v != s && v !=t ) 
   160 	      if ( excess.get(v)==0 && v != s && v !=t ) 
   172 		stack[level.get(v)].push(v); 
   161 		stack[level.get(v)].push(v); 
   173 	      /*v becomes active.*/
   162 	      /*v becomes active.*/
   174 	      
   163 	      
   175 	      int cap=capacity.get(e);
   164 	      T cap=capacity.get(e);
   176 	      int flo=flow.get(e);
   165 	      T flo=flow.get(e);
   177 	      int remcap=cap-flo;
   166 	      T remcap=cap-flo;
   178 	      
   167 	      
   179 	      if ( remcap >= exc ) {       
   168 	      if ( remcap >= exc ) {       
   180 		/*A nonsaturating push.*/
   169 		/*A nonsaturating push.*/
   181 		
   170 		
   182 		flow.set(e, flo+exc);
   171 		flow.set(e, flo+exc);
   208 	      
   197 	      
   209 	      if ( excess.get(v)==0 && v != s && v !=t) 
   198 	      if ( excess.get(v)==0 && v != s && v !=t) 
   210 		stack[level.get(v)].push(v); 
   199 		stack[level.get(v)].push(v); 
   211 	      /*v becomes active.*/
   200 	      /*v becomes active.*/
   212 	      
   201 	      
   213 	      int flo=flow.get(e);
   202 	      T flo=flow.get(e);
   214 	      
   203 	      
   215 	      if ( flo >= exc ) { 
   204 	      if ( flo >= exc ) { 
   216 		/*A nonsaturating push.*/
   205 		/*A nonsaturating push.*/
   217 		
   206 		
   218 		flow.set(e, flo-exc);
   207 		flow.set(e, flo-exc);
   229 	    } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
   218 	    } else if ( newlevel > level.get(v) ) newlevel = level.get(v);
   230 	    
   219 	    
   231 	  } //for in edges vw
   220 	  } //for in edges vw
   232 	  
   221 	  
   233 	} // if w still has excess after the out edge for cycle
   222 	} // if w still has excess after the out edge for cycle
   234 	 
   223 	
   235 	  excess.set(w, exc);
   224 	excess.set(w, exc);
       
   225 	
       
   226 
       
   227 	
       
   228 
       
   229 	/*
       
   230 	  Relabel
       
   231 	*/
       
   232 	
       
   233 	if ( exc > 0 ) {
       
   234 	  //now 'lev' is the old level of w
       
   235 	  level.set(w,++newlevel);
   236 	  
   236 	  
   237 
   237 	  if ( lev < n ) {
   238 	  /*
   238 	    --numb[lev];
   239 	    Relabel
   239 	    
   240 	  */
   240 	    if ( !numb[lev] && lev < A*n ) {  //If the level of w gets empty. 
       
   241 	      
       
   242 	      for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
       
   243 		if (level.get(v) > lev && level.get(v) < n ) level.set(v,n);  
       
   244 	      }
       
   245 	      for (int i=lev+1 ; i!=n ; ++i) numb[i]=0; 
       
   246 	      if ( newlevel < n ) newlevel=n; 
       
   247 	    } else { 
       
   248 	      if ( newlevel < n ) ++numb[newlevel]; 
       
   249 	    }
       
   250 	  } 
   241 	  
   251 	  
   242 	  if ( exc > 0 ) {
   252 	  stack[newlevel].push(w);
   243 	    //now 'lev' is the old level of w
   253 	  b=newlevel;
   244 	    level.set(w,++newlevel);
   254 	  
   245 	    
   255 	}
   246 	    if ( lev < n ) {
   256 	
   247 	      --numb[lev];
       
   248 
       
   249 	      if ( !numb[lev] && lev < A*n ) {  //If the level of w gets empty. 
       
   250 		
       
   251 		for (EachNodeIt v=G.template first<EachNodeIt>(); v.valid() ; ++v) {
       
   252 		  if (level.get(v) > lev && level.get(v) < n ) level.set(v,n);  
       
   253 		}
       
   254 		for (int i=lev+1 ; i!=n ; ++i) numb[i]=0; 
       
   255 		if ( newlevel < n ) newlevel=n; 
       
   256 	      } else { 
       
   257 		if ( newlevel < n ) ++numb[newlevel]; 
       
   258 	      }
       
   259 	    } 
       
   260 	    
       
   261 	    stack[newlevel].push(w);
       
   262 
       
   263 	  }
       
   264 
       
   265 	} // if stack[b] is nonempty
       
   266 
       
   267       } // while(b)
   257       } // while(b)
   268 
   258       
   269 
   259       
   270       value = excess.get(t);
   260       value = excess.get(t);
   271       /*Max flow value.*/
   261       /*Max flow value.*/
   272 
   262 
   273 
   263 
   274     } //void run()
   264     } //void run()
   279 
   269 
   280     /*
   270     /*
   281       Returns the maximum value of a flow.
   271       Returns the maximum value of a flow.
   282      */
   272      */
   283 
   273 
   284     T maxflow() {
   274     T maxFlow() {
   285       return value;
   275       return value;
   286     }
   276     }
   287 
   277 
   288 
   278 
   289 
   279 
   290     /*
   280     /*
   291       For the maximum flow x found by the algorithm, it returns the flow value on Edge e, i.e. x(e). 
   281       For the maximum flow x found by the algorithm, 
       
   282       it returns the flow value on edge e, i.e. x(e). 
   292     */
   283     */
   293 
   284 
   294     T flowonedge(EdgeIt e) {
   285     T flowOnEdge(const EdgeIt e) {
   295       return flow.get(e);
   286       return flow.get(e);
   296     }
   287     }
   297 
   288 
   298 
   289 
   299 
   290 
   300     /*
   291     /*
   301       Returns the maximum flow x found by the algorithm.
   292       Returns the maximum flow x found by the algorithm.
   302     */
   293     */
   303 
   294 
   304     FlowMap allflow() {
   295     FlowMap Flow() {
   305       return flow;
   296       return flow;
   306     }
   297     }
   307 
   298 
   308 
   299 
   309 
   300 
   311     /*
   302     /*
   312       Returns the minimum min cut, by a bfs from s in the residual graph.
   303       Returns the minimum min cut, by a bfs from s in the residual graph.
   313     */
   304     */
   314     
   305     
   315     template<typename CutMap>
   306     template<typename CutMap>
   316     void mincut(CutMap& M) {
   307     void minCut(CutMap& M) {
   317     
   308     
   318       std::queue<NodeIt> queue;
   309       std::queue<NodeIt> queue;
   319       
   310       
   320       M.set(s,true);      
   311       M.set(s,true);      
   321       queue.push(s);
   312       queue.push(s);
   322 
   313 
   323       while (!queue.empty()) {
   314       while (!queue.empty()) {
   324         NodeIt w=queue.front();
   315         NodeIt w=queue.front();
   325 	queue.pop();
   316 	queue.pop();
   326 
   317 	
   327 	for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
   318 	for(OutEdgeIt e=G.template first<OutEdgeIt>(w) ; e.valid(); ++e) {
   328 	  NodeIt v=G.head(e);
   319 	  NodeIt v=G.head(e);
   329 	  if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
   320 	  if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
   330 	    queue.push(v);
   321 	    queue.push(v);
   331 	    M.set(v, true);
   322 	    M.set(v, true);
   336 	  NodeIt v=G.tail(e);
   327 	  NodeIt v=G.tail(e);
   337 	  if (!M.get(v) && flow.get(e) > 0 ) {
   328 	  if (!M.get(v) && flow.get(e) > 0 ) {
   338 	    queue.push(v);
   329 	    queue.push(v);
   339 	    M.set(v, true);
   330 	    M.set(v, true);
   340 	  }
   331 	  }
   341 	} 
   332 	}
   342 
   333 
   343       }
   334       }
   344 
       
   345     }
   335     }
   346 
   336 
   347 
   337 
   348 
   338 
   349     /*
   339     /*
   350       Returns the maximum min cut, by a reverse bfs 
   340       Returns the maximum min cut, by a reverse bfs 
   351       from t in the residual graph.
   341       from t in the residual graph.
   352     */
   342     */
   353     
   343     
   354     template<typename CutMap>
   344     template<typename CutMap>
   355     void max_mincut(CutMap& M) {
   345     void maxMinCut(CutMap& M) {
   356     
   346     
   357       std::queue<NodeIt> queue;
   347       std::queue<NodeIt> queue;
   358       
   348       
   359       M.set(t,true);        
   349       M.set(t,true);        
   360       queue.push(t);
   350       queue.push(t);
   387     }
   377     }
   388 
   378 
   389 
   379 
   390 
   380 
   391     template<typename CutMap>
   381     template<typename CutMap>
   392     void min_mincut(CutMap& M) {
   382     void minMinCut(CutMap& M) {
   393       mincut(M);
   383       minCut(M);
   394     }
   384     }
   395 
   385 
   396 
   386 
   397 
   387 
   398   };
   388   };
   399 }//namespace hugo
   389 }//namespace marci
   400 #endif 
   390 #endif 
   401 
   391 
   402 
   392 
   403 
   393 
   404 
   394