Index: src/work/jacint/preflow_hl4.h
===================================================================
--- src/work/jacint/preflow_hl4.h	(revision 105)
+++ src/work/jacint/preflow_hl4.h	(revision 109)
@@ -1,38 +1,48 @@
 // -*- C++ -*-
 /*
-preflow_hl4.h
+preflow_h5.h
 by jacint. 
-Runs the two phase highest label preflow push algorithm. In phase 0
-we maintain in a list the nodes in level i < n, and we maintain a 
-bound k on the max level i < n containing a node, so we can do
-the gap heuristic fast. Phase 1 is the same. (The algorithm is the 
-same as preflow.hl3, the only diff is that here we use the gap
-heuristic with the list of the nodes on level i, and not a bfs form the
-upgraded node.)
-
-In phase 1 we shift everything downwards by n.
-
-Member functions:
-
-void run() : runs the algorithm
-
- The following functions should be used after run() was already run.
-
-T maxflow() : returns the value of a maximum flow
-
-T flowonedge(EdgeIt e) : for a fixed maximum flow x it returns x(e) 
-
-FlowMap allflow() : returns a maximum flow
-
-void allflow(FlowMap& _flow ) : returns a maximum flow
-
-void mincut(CutMap& M) : sets M to the characteristic vector of a 
-     minimum cut. M should be a map of bools initialized to false.
-
-void min_mincut(CutMap& M) : sets M to the characteristic vector of the 
+Heuristics: 
+ 2 phase
+ gap
+ list 'level_list' on the nodes on level i implemented by hand
+ highest label
+ relevel: in phase 0, after BFS*n relabels, it runs a reverse 
+   bfs from t in the res graph to relevel the nodes reachable from t.
+   BFS is initialized to 20
+
+Due to the last heuristic, this algorithm is quite fast on very 
+sparse graphs, but relatively bad on even the dense graphs.
+
+'NodeMap<bool> cut' is a member, in this way we can count it fast, after
+the algorithm was run.
+
+The constructor runs the algorithm.
+
+Members:
+
+T maxFlow() : returns the value of a maximum flow
+
+T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e) 
+
+FlowMap Flow() : returns the fixed maximum flow x
+
+void Flow(FlowMap& _flow ) : returns the fixed maximum flow x
+
+void minMinCut(CutMap& M) : sets M to the characteristic vector of the 
      minimum min cut. M should be a map of bools initialized to false.
 
-void max_mincut(CutMap& M) : sets M to the characteristic vector of the 
+void maxMinCut(CutMap& M) : sets M to the characteristic vector of the 
      maximum min cut. M should be a map of bools initialized to false.
+
+void minCut(CutMap& M) : fast function, sets M to the characteristic 
+     vector of a minimum cut. 
+
+Different member from the other preflow_hl-s (here we have a member 
+'NodeMap<bool> cut').
+
+CutMap minCut() : fast function, giving the characteristic 
+     vector of a minimum cut.
+
 
 */
@@ -41,7 +51,10 @@
 #define PREFLOW_HL4_H
 
+#define BFS 20
+
 #include <vector>
-#include <stack>
 #include <queue>
+
+#include <time_measure.h>  //for test
 
 namespace hugo {
@@ -49,4 +62,5 @@
   template <typename Graph, typename T, 
     typename FlowMap=typename Graph::EdgeMap<T>, 
+    typename CutMap=typename Graph::NodeMap<bool>,
     typename CapMap=typename Graph::EdgeMap<T> >
   class preflow_hl4 {
@@ -63,16 +77,19 @@
     FlowMap flow;
     CapMap& capacity;  
+    CutMap cut;
     T value;
     
   public:
 
+    double time;
+
     preflow_hl4(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity) :
-      G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity) { }
-
-
-    void run() {
- 
-      bool phase=0;
+      G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity), 
+      cut(G, false) {
+
+      bool phase=0;   //phase 0 is the 1st phase, phase 1 is the 2nd
       int n=G.nodeNum(); 
+      int relabel=0;
+      int heur=(int)BFS*n;
       int k=n-2;
       int b=k;
@@ -84,5 +101,7 @@
       typename Graph::NodeMap<int> level(G,n);      
       typename Graph::NodeMap<T> excess(G); 
-      std::vector<std::stack<NodeIt> > stack(n);    
+      
+      std::vector<NodeIt> active(n);
+      typename Graph::NodeMap<NodeIt> next(G);
       //Stack of the active nodes in level i < n.
       //We use it in both phases.
@@ -108,5 +127,5 @@
 	for(InEdgeIt e=G.template first<InEdgeIt>(v); e.valid(); ++e) {
 	  NodeIt w=G.tail(e);
-	  if ( level.get(w) == n ) {
+	  if ( level.get(w) == n && w !=s ) {
 	    bfs_queue.push(w);
 	    NodeIt first=level_list[l];
@@ -129,5 +148,8 @@
 	  NodeIt w=G.head(e);
 	  if ( level.get(w) < n ) {	  
-	    if ( excess.get(w) == 0 && w!=t ) stack[level.get(w)].push(w); 
+	    if ( excess.get(w) == 0 && w!=t ) {
+	      next.set(w,active[level.get(w)]);
+	      active[level.get(w)]=w;
+	    }
 	    flow.set(e, c); 
 	    excess.set(w, excess.get(w)+c);
@@ -152,4 +174,11 @@
 	  */
 	  phase=1;
+	  
+	  //Now have a min cut.
+	  for( EachNodeIt v=G.template first<EachNodeIt>(); 
+	       v.valid(); ++v) 
+	    if (level.get(v) >= n ) cut.set(v,true);
+	  
+	  time=currTime();
 	  level.set(s,0);
 	  std::queue<NodeIt> bfs_queue;
@@ -168,5 +197,8 @@
 		bfs_queue.push(u);
 		level.set(u, l);
-		if ( excess.get(u) > 0 ) stack[l].push(u);
+		if ( excess.get(u) > 0 ) {
+		    next.set(u,active[l]);
+		    active[l]=u;
+		}
 	      }
 	    }
@@ -178,5 +210,8 @@
 		bfs_queue.push(u);
 		level.set(u, l);
-		if ( excess.get(u) > 0 ) stack[l].push(u);
+		if ( excess.get(u) > 0 ) {
+		    next.set(u,active[l]);
+		    active[l]=u;
+		}
 	      }
 	    }
@@ -186,12 +221,12 @@
 
 
-	if ( stack[b].empty() ) --b;
+	if ( active[b] == 0 ) --b;
 	else {
 	  
-	  NodeIt w=stack[b].top();        //w is a highest label active node.
-	  stack[b].pop();           
+	  NodeIt w=active[b];
+	  active[b]=next.get(w);
 	  int lev=level.get(w);
 	  T exc=excess.get(w);
-	  int newlevel=n;          //In newlevel we bound the next level of w.
+	  int newlevel=n;          //bound on the next level of w.
 	  
 	  for(OutEdgeIt e=G.template first<OutEdgeIt>(w); e.valid(); ++e) {
@@ -204,7 +239,9 @@
 	      /*Push is allowed now*/
 	      
-	      if ( excess.get(v)==0 && v!=t && v!=s ) 
-		stack[level.get(v)].push(v); 
-	      /*v becomes active.*/
+	      if ( excess.get(v)==0 && v!=t && v!=s )  {
+		int lev_v=level.get(v);
+		next.set(v,active[lev_v]);
+		active[lev_v]=v;
+	      }
 	      
 	      T cap=capacity.get(e);
@@ -244,7 +281,9 @@
 	      /*Push is allowed now*/
 	      
-	      if ( excess.get(v)==0 && v!=t && v!=s ) 
-		stack[level.get(v)].push(v); 
-	      /*v becomes active.*/
+	      if ( excess.get(v)==0 && v!=t && v!=s ) {
+		int lev_v=level.get(v);
+		next.set(v,active[lev_v]);
+		active[lev_v]=v;
+	      }
 	      
 	      T flo=flow.get(e);
@@ -282,5 +321,6 @@
 	  if ( phase ) {
 	    level.set(w,++newlevel);
-	    stack[newlevel].push(w);
+	    next.set(w,active[newlevel]);
+	    active[newlevel]=w;
 	    b=newlevel;
 	  } else {
@@ -304,4 +344,5 @@
 	      } 
 	    }
+	    //unlacing ends
 		
 
@@ -329,5 +370,6 @@
 		
 		level.set(w,++newlevel);
-		stack[newlevel].push(w);
+		next.set(w,active[newlevel]);
+		active[newlevel]=w;
 		b=newlevel;
 		if ( k < newlevel ) ++k;
@@ -339,7 +381,76 @@
 	      }
 	    }
+
+	    ++relabel;
+	    if ( relabel >= heur ) {
+	      relabel=0;
+	      b=n-2;
+	      k=b;
+		
+	      for ( int i=1; i!=n; ++i ) { 
+		active[i]=0;
+		level_list[i]=0;
+	      }
+
+	      //bfs from t in the res graph to relevel the nodes
+	      for( EachNodeIt v=G.template first<EachNodeIt>(); 
+		   v.valid(); ++v) level.set(v,n);
+
+	      level.set(t,0);
+	      std::queue<NodeIt> bfs_queue;
+	      bfs_queue.push(t);
+	      
+	      while (!bfs_queue.empty()) {
+		
+		NodeIt v=bfs_queue.front();	
+		bfs_queue.pop();
+		int l=level.get(v)+1;
+		
+		for(InEdgeIt e=G.template first<InEdgeIt>(v); 
+		    e.valid(); ++e) {
+		  if ( capacity.get(e) == flow.get(e) ) continue;
+		  NodeIt u=G.tail(e);
+		  if ( level.get(u) == n ) { 
+		    bfs_queue.push(u);
+		    level.set(u, l);
+		    if ( excess.get(u) > 0 ) {
+		      next.set(u,active[l]);
+		      active[l]=u;
+		    }
+		    NodeIt first=level_list[l];
+		    if ( first != 0 ) left.set(first,w);
+		    right.set(w,first);
+		    left.set(w,0);
+		    level_list[l]=w;
+		  }
+		}
+		
+		
+		for(OutEdgeIt e=G.template first<OutEdgeIt>(v); 
+		    e.valid(); ++e) {
+		  if ( 0 == flow.get(e) ) continue;
+		  NodeIt u=G.head(e);
+		  if ( level.get(u) == n ) { 
+		    bfs_queue.push(u);
+		    level.set(u, l);
+		    if ( excess.get(u) > 0 ) {
+		      next.set(u,active[l]);
+		      active[l]=u;
+		    }
+		    NodeIt first=level_list[l];
+		    if ( first != 0 ) left.set(first,w);
+		    right.set(w,first);
+		    left.set(w,0);
+		    level_list[l]=w;
+		  }
+		}
+	      }
+	      
+	      level.set(s,n);
+	    }
+	  
 	  } //phase 0
 	} // if ( exc > 0 )
- 
+	
 	
 	} // if stack[b] is nonempty
@@ -362,5 +473,5 @@
      */
 
-    T maxflow() {
+    T maxFlow() {
       return value;
     }
@@ -372,5 +483,5 @@
     */
 
-    T flowonedge(EdgeIt e) {
+    T flowOnEdge(EdgeIt e) {
       return flow.get(e);
     }
@@ -378,5 +489,5 @@
 
 
-    FlowMap allflow() {
+    FlowMap Flow() {
       return flow;
     }
@@ -384,5 +495,5 @@
 
     
-    void allflow(FlowMap& _flow ) {
+    void Flow(FlowMap& _flow ) {
       for(EachNodeIt v=G.template first<EachNodeIt>() ; v.valid(); ++v)
 	_flow.set(v,flow.get(v));
@@ -395,6 +506,6 @@
     */
     
-    template<typename CutMap>
-    void mincut(CutMap& M) {
+    template<typename _CutMap>
+    void minMinCut(_CutMap& M) {
     
       std::queue<NodeIt> queue;
@@ -434,6 +545,6 @@
     */
     
-    template<typename CutMap>
-    void max_mincut(CutMap& M) {
+    template<typename _CutMap>
+    void maxMinCut(_CutMap& M) {
     
       std::queue<NodeIt> queue;
@@ -470,14 +581,19 @@
 
 
-
-    template<typename CutMap>
-    void min_mincut(CutMap& M) {
-      mincut(M);
-    }
-
+    template<typename _CutMap>
+    void minCut(_CutMap& M) {
+      for( EachNodeIt v=G.template first<EachNodeIt>(); 
+	   v.valid(); ++v) 
+	M.set(v, cut.get(v));
+    }
+
+   
+    CutMap minCut() {
+      return cut;
+    }
 
 
   };
-}//namespace hugo
+}//namespace marci
 #endif 
 
