# HG changeset patch
# User jacint
# Date 1083106755 0
# Node ID 6b36be4cffa40d1b74d6560e10b2155b7aac5d27
# Parent  5caac2f7829b90163fbc268dad187530747bdbde
Changes in the interface and new test program added.

diff -r 5caac2f7829b -r 6b36be4cffa4 src/work/jacint/preflow.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/jacint/preflow.cc	Tue Apr 27 22:59:15 2004 +0000
@@ -0,0 +1,239 @@
+#include <iostream>
+
+#include <smart_graph.h>
+#include <dimacs.h>
+#include <preflow.h>
+#include <time_measure.h>
+
+using namespace hugo;
+
+int main(int, char **) {
+ 
+  typedef SmartGraph Graph;
+  
+  typedef Graph::Node Node;
+  typedef Graph::EdgeIt EdgeIt;
+
+  Graph G;
+  Node s, t;
+  Graph::EdgeMap<int> cap(G);
+  readDimacsMaxFlow(std::cin, G, s, t, cap);
+  Timer ts;
+  
+  std::cout <<
+    "\n  Testing preflow.h on a graph with " << 
+    G.nodeNum() << " nodes and " << G.edgeNum() << " edges..."
+	   << std::endl;
+
+
+  Graph::EdgeMap<int> flow(G,0);
+  Preflow<Graph, int> preflow_test(G, s, t, cap, flow);
+  std::cout << "\nCalling run() (flow must be constant zero)..."<<std::endl;
+  ts.reset();
+  preflow_test.run();
+  std::cout << "Elapsed time: " << ts << std::endl;
+
+  Graph::NodeMap<bool> mincut(G);
+  preflow_test.minMinCut(mincut); 
+  int min_min_cut_value=0;
+  Graph::NodeMap<bool> cut(G);
+  preflow_test.minCut(cut); 
+  int min_cut_value=0;
+  Graph::NodeMap<bool> maxcut(G);
+  preflow_test.maxMinCut(maxcut); 
+  int max_min_cut_value=0;
+  EdgeIt e;
+  for(G.first(e); G.valid(e); G.next(e)) {
+    int c=cap[e];
+    if (mincut[G.tail(e)] && !mincut[G.head(e)]) min_min_cut_value+=c;
+    if (cut[G.tail(e)] && !cut[G.head(e)]) min_cut_value+=c; 
+    if (maxcut[G.tail(e)] && !maxcut[G.head(e)]) max_min_cut_value+=c;
+  }
+
+  std::cout << "\nChecking the result: " <<std::endl;  
+  std::cout << "Flow value: "<< preflow_test.flowValue() << std::endl;
+  std::cout << "Min cut value: "<< min_cut_value << std::endl;
+  std::cout << "Min min cut value: "<< min_min_cut_value << std::endl;
+  std::cout << "Max min cut value: "<< max_min_cut_value << 
+    std::endl;
+
+  if ( preflow_test.flowValue() == min_cut_value &&
+       min_cut_value == min_min_cut_value &&
+       min_min_cut_value == max_min_cut_value )
+    std::cout << "They are equal. " <<std::endl;  
+
+
+
+
+
+  Preflow<Graph, int> preflow_test2(G, s, t, cap, flow);
+  std::cout << "\n\nCalling preflow(GEN_FLOW) with the given maximum flow..."<<std::endl;
+  ts.reset();
+  preflow_test2.preflow(preflow_test2.GEN_FLOW);
+  std::cout << "Elapsed time: " << ts << std::endl;
+
+  Graph::NodeMap<bool> mincut2(G);
+  preflow_test.minMinCut(mincut2); 
+  int min_min_cut2_value=0;
+  Graph::NodeMap<bool> cut2(G);
+  preflow_test.minCut(cut2); 
+  int min_cut2_value=0;
+  Graph::NodeMap<bool> maxcut2(G);
+  preflow_test.maxMinCut(maxcut2); 
+  int max_min_cut2_value=0;
+  for(G.first(e); G.valid(e); G.next(e)) {
+    int c=cap[e];
+    if (mincut2[G.tail(e)] && !mincut2[G.head(e)]) min_min_cut2_value+=c;
+    if (cut2[G.tail(e)] && !cut2[G.head(e)]) min_cut2_value+=c; 
+    if (maxcut2[G.tail(e)] && !maxcut2[G.head(e)]) max_min_cut2_value+=c;
+  }
+
+  std::cout << "\nThe given flow value is "
+	    << preflow_test2.flowValue();
+
+  if ( preflow_test2.flowValue() == min_cut2_value &&
+       min_cut2_value == min_min_cut2_value &&
+       min_min_cut2_value == max_min_cut2_value )
+    std::cout <<", which is equal to all three min cut values." 
+	      <<std::endl;  
+
+
+
+
+
+  Graph::EdgeMap<int> flow3(G,0);
+  Preflow<Graph, int> preflow_test3(G, s, t, cap, flow3);
+  std::cout << "\n\nCalling preflowPhase0(PREFLOW) on the constant zero flow..."<<std::endl;
+  ts.reset();
+  preflow_test3.preflowPhase0(preflow_test3.PREFLOW);
+  std::cout << "Elapsed time: " << ts << std::endl;
+  Graph::NodeMap<bool> actcut3(G);
+  std::cout << "\nCalling actMinCut()..."<<std::endl;
+  preflow_test3.actMinCut(actcut3); 
+  std::cout << "Calling preflowPhase1() on the given flow..."<<std::endl;
+  ts.reset();
+  preflow_test3.preflowPhase1();
+  std::cout << "Elapsed time: " << ts << std::endl;
+  
+  int act_min_cut3_value=0;
+  
+  Graph::NodeMap<bool> mincut3(G);
+  preflow_test.minMinCut(mincut3); 
+  int min_min_cut3_value=0;
+  
+  Graph::NodeMap<bool> cut3(G);
+  preflow_test.minCut(cut3); 
+  int min_cut3_value=0;
+  
+  Graph::NodeMap<bool> maxcut3(G);
+  preflow_test.maxMinCut(maxcut3); 
+  int max_min_cut3_value=0;
+  
+  for(G.first(e); G.valid(e); G.next(e)) {
+    int c=cap[e];
+    if (mincut3[G.tail(e)] && !mincut3[G.head(e)]) min_min_cut3_value+=c;
+    if (cut3[G.tail(e)] && !cut3[G.head(e)]) min_cut3_value+=c; 
+    if (maxcut3[G.tail(e)] && !maxcut3[G.head(e)]) max_min_cut3_value+=c;
+    if (actcut3[G.tail(e)] && !actcut3[G.head(e)]) act_min_cut3_value+=c;
+  }
+
+ std::cout << "\nThe min cut value given by actMinCut() after phase 0 is "<<
+   act_min_cut3_value;
+
+  if ( preflow_test3.flowValue() == min_cut3_value &&
+       min_cut3_value == min_min_cut3_value &&
+       min_min_cut3_value == max_min_cut3_value &&
+       max_min_cut3_value == act_min_cut3_value ) {
+    std::cout << 
+      ", which is equal to the given flow value and to all three min cut values after phase 1." 
+	      <<std::endl;  
+  }
+
+
+
+
+
+  Graph::EdgeMap<int> flow4(G,0);
+  Preflow<Graph, int> preflow_test4(G, s, t, cap, flow4);
+  std::cout << 
+    "\n\nCalling preflow(PREFLOW) with the constant 0 flow, the result is f..."
+	    <<std::endl;
+  preflow_test4.preflow(preflow_test4.PREFLOW);
+
+  std::cout << "Swapping the source and the target, "<<std::endl;
+  std::cout << "by calling resetSource(t) and resetTarget(s)..."
+	    <<std::endl;
+  preflow_test4.resetSource(t);
+  preflow_test4.resetTarget(s);
+
+  std::cout << 
+    "Calling preflow(PREFLOW) to find a maximum t-s flow starting with flow f..."
+	    <<std::endl;
+  preflow_test4.preflow(preflow_test4.PREFLOW);
+
+  Graph::NodeMap<bool> mincut4(G);
+  preflow_test4.minMinCut(mincut4); 
+  int min_min_cut4_value=0;
+  Graph::NodeMap<bool> cut4(G);
+  preflow_test4.minCut(cut4); 
+  int min_cut4_value=0;
+  Graph::NodeMap<bool> maxcut4(G);
+  preflow_test4.maxMinCut(maxcut4); 
+  int max_min_cut4_value=0;
+  for(G.first(e); G.valid(e); G.next(e)) {
+    int c=cap[e];
+    if (mincut4[G.tail(e)] && !mincut4[G.head(e)]) min_min_cut4_value+=c;
+    if (cut4[G.tail(e)] && !cut4[G.head(e)]) min_cut4_value+=c; 
+    if (maxcut4[G.tail(e)] && !maxcut4[G.head(e)]) max_min_cut4_value+=c;
+  }
+
+  std::cout << "\nThe given flow value is "
+	    << preflow_test4.flowValue();
+  
+  if ( preflow_test4.flowValue() == min_cut4_value &&
+       min_cut4_value == min_min_cut4_value &&
+       min_min_cut4_value == max_min_cut4_value )
+    std::cout <<", which is equal to all three min cut values." 
+	      <<std::endl;  
+
+
+
+
+  Graph::EdgeMap<int> flow5(G,0);
+  std::cout << "Resetting the stored flow to constant zero, by calling resetFlow..."
+	    <<std::endl;
+  preflow_test4.resetFlow(flow5);
+  std::cout << 
+    "Calling preflow(GEN_FLOW) to find a maximum t-s flow "<<std::endl;
+  std::cout << 
+    "starting with this constant zero flow..." <<std::endl;
+  preflow_test4.preflow(preflow_test4.GEN_FLOW);
+
+  Graph::NodeMap<bool> mincut5(G);
+  preflow_test4.minMinCut(mincut5); 
+  int min_min_cut5_value=0;
+  Graph::NodeMap<bool> cut5(G);
+  preflow_test4.minCut(cut5); 
+  int min_cut5_value=0;
+  Graph::NodeMap<bool> maxcut5(G);
+  preflow_test4.maxMinCut(maxcut5); 
+  int max_min_cut5_value=0;
+  for(G.first(e); G.valid(e); G.next(e)) {
+    int c=cap[e];
+    if (mincut5[G.tail(e)] && !mincut5[G.head(e)]) min_min_cut5_value+=c;
+    if (cut5[G.tail(e)] && !cut5[G.head(e)]) min_cut5_value+=c; 
+    if (maxcut5[G.tail(e)] && !maxcut5[G.head(e)]) max_min_cut5_value+=c;
+  }
+
+  std::cout << "\nThe given flow value is "
+	    << preflow_test4.flowValue();
+  
+  if ( preflow_test4.flowValue() == min_cut5_value &&
+       min_cut5_value == min_min_cut5_value &&
+       min_min_cut5_value == max_min_cut5_value )
+    std::cout <<", which is equal to all three min cut values." 
+	      <<std::endl<<std::endl;  
+
+
+  return 0;
+}
diff -r 5caac2f7829b -r 6b36be4cffa4 src/work/jacint/preflow.h
--- a/src/work/jacint/preflow.h	Tue Apr 27 22:29:11 2004 +0000
+++ b/src/work/jacint/preflow.h	Tue Apr 27 22:59:15 2004 +0000
@@ -1,20 +1,15 @@
 // -*- C++ -*-
 
-//run gyorsan tudna adni a minmincutot a 2 fazis elejen , ne vegyuk be konstruktorba egy cutmapet?
-//constzero jo igy?
-
-//majd marci megmondja betegyem-e bfs-t meg resgraphot
-
 /*
 Heuristics: 
  2 phase
  gap
  list 'level_list' on the nodes on level i implemented by hand
- stack 'active' on the active nodes on level i implemented by hand
+ stack 'active' on the active nodes on level i
  runs heuristic 'highest label' for H1*n relabels
  runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
  
-Parameters H0 and H1 are initialized to 20 and 10.
+Parameters H0 and H1 are initialized to 20 and 1.
 
 Constructors:
 
@@ -28,7 +23,7 @@
 T flowValue() : returns the value of a maximum flow
 
 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.
+     minimum min cut. M should be a map of bools initialized to false. ??Is it OK?
 
 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.
@@ -36,8 +31,6 @@
 void minCut(CutMap& M) : sets M to the characteristic vector of 
      a min cut. M should be a map of bools initialized to false.
 
-FIXME reset
-
 */
 
 #ifndef HUGO_PREFLOW_H
@@ -48,6 +41,7 @@
 
 #include <vector>
 #include <queue>
+#include <stack>
 
 namespace hugo {
 
@@ -57,494 +51,223 @@
   class Preflow {
     
     typedef typename Graph::Node Node;
-    typedef typename Graph::Edge Edge;
     typedef typename Graph::NodeIt NodeIt;
     typedef typename Graph::OutEdgeIt OutEdgeIt;
     typedef typename Graph::InEdgeIt InEdgeIt;
-    
+
+    typedef typename std::vector<std::stack<Node> > VecStack;
+    typedef typename Graph::template NodeMap<Node> NNMap;
+    typedef typename std::vector<Node> VecNode;
+
     const Graph& G;
     Node s;
     Node t;
-    const CapMap& capacity;  
-    FlowMap& flow;
-    T value;
-    bool constzero;
+    CapMap* capacity;  
+    FlowMap* flow;
+    int n;      //the number of nodes of G
+    typename Graph::template NodeMap<int> level;      
+    typename Graph::template NodeMap<T> excess; 
+
 
   public:
+ 
+    enum flowEnum{
+      ZERO_FLOW=0,
+      GEN_FLOW=1,
+      PREFLOW=2
+    };
+
     Preflow(Graph& _G, Node _s, Node _t, CapMap& _capacity, 
-	    FlowMap& _flow, bool _constzero ) :
-      G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow), constzero(_constzero) {}
+	    FlowMap& _flow) :
+      G(_G), s(_s), t(_t), capacity(&_capacity), 
+      flow(&_flow), n(_G.nodeNum()), level(_G), excess(_G,0) {}
+
+    void run() {
+      preflow( ZERO_FLOW );
+    }
     
-    
-    void run() {
+    void preflow( flowEnum fe ) {
+      preflowPhase0(fe);
+      preflowPhase1();
+    }
+
+    void preflowPhase0( flowEnum fe ) {
       
-      value=0;                //for the subsequent runs
-
-      bool phase=0;        //phase 0 is the 1st phase, phase 1 is the 2nd
-      int n=G.nodeNum(); 
       int heur0=(int)(H0*n);  //time while running 'bound decrease' 
       int heur1=(int)(H1*n);  //time while running 'highest label'
       int heur=heur1;         //starting time interval (#of relabels)
+      int numrelabel=0;
+     
       bool what_heur=1;       
-      /*
-	what_heur is 0 in case 'bound decrease' 
-	and 1 in case 'highest label'
-      */
+      //It is 0 in case 'bound decrease' and 1 in case 'highest label'
+
       bool end=false;     
-      /*
-	Needed for 'bound decrease', 'true'
-	means no active nodes are above bound b.
-      */
-      int relabel=0;
+      //Needed for 'bound decrease', true means no active nodes are above bound b.
+
       int k=n-2;  //bound on the highest level under n containing a node
       int b=k;    //bound on the highest level under n of an active node
       
-      typename Graph::template NodeMap<int> level(G,n);      
-      typename Graph::template NodeMap<T> excess(G); 
+      VecStack active(n);
+      
+      NNMap left(G,INVALID);
+      NNMap right(G,INVALID);
+      VecNode level_list(n,INVALID);
+      //List of the nodes in level i<n, set to n.
 
-      std::vector<Node> active(n-1,INVALID);
-      typename Graph::template NodeMap<Node> next(G,INVALID);
-      //Stack of the active nodes in level i < n.
-      //We use it in both phases.
-
-      typename Graph::template NodeMap<Node> left(G,INVALID);
-      typename Graph::template NodeMap<Node> right(G,INVALID);
-      std::vector<Node> level_list(n,INVALID);
-      /*
-	List of the nodes in level i<n.
-      */
-
-
-      if ( constzero ) {
-     
-	/*Reverse_bfs from t, to find the starting level.*/
-	level.set(t,0);
-	std::queue<Node> bfs_queue;
-	bfs_queue.push(t);
-	
-	while (!bfs_queue.empty()) {
+      NodeIt v;
+      for(G.first(v); G.valid(v); G.next(v)) level.set(v,n);
+      //setting each node to level n
+      
+      switch ( fe ) {
+      case PREFLOW:
+	{
+	  //counting the excess
+	  NodeIt v;
+	  for(G.first(v); G.valid(v); G.next(v)) {
+	    T exc=0;
 	  
-	  Node v=bfs_queue.front();	
-	  bfs_queue.pop();
-	  int l=level[v]+1;
+	    InEdgeIt e;
+	    for(G.first(e,v); G.valid(e); G.next(e)) exc+=(*flow)[e];
+	    OutEdgeIt f;
+	    for(G.first(f,v); G.valid(f); G.next(f)) exc-=(*flow)[f];
+	    
+	    excess.set(v,exc);	  
+	    
+	    //putting the active nodes into the stack
+	    int lev=level[v];
+	    if ( exc > 0 && lev < n && v != t ) active[lev].push(v);
+	  }
+	  break;
+	}
+      case GEN_FLOW:
+	{
+	  //Counting the excess of t
+	  T exc=0;
 	  
 	  InEdgeIt e;
-	  for(G.first(e,v); G.valid(e); G.next(e)) {
-	    Node w=G.tail(e);
-	    if ( level[w] == n && w != s ) {
-	      bfs_queue.push(w);
-	      Node first=level_list[l];
-	      if ( G.valid(first) ) left.set(first,w);
-	      right.set(w,first);
-	      level_list[l]=w;
-	      level.set(w, l);
+	  for(G.first(e,t); G.valid(e); G.next(e)) exc+=(*flow)[e];
+	  OutEdgeIt f;
+	  for(G.first(f,t); G.valid(f); G.next(f)) exc-=(*flow)[f];
+	  
+	  excess.set(t,exc);	
+	  
+	  break;
+	}
+      }
+      
+      preflowPreproc( fe, active, level_list, left, right );
+      //End of preprocessing 
+      
+      
+      //Push/relabel on the highest level active nodes.
+      while ( true ) {
+	if ( b == 0 ) {
+	  if ( !what_heur && !end && k > 0 ) {
+	    b=k;
+	    end=true;
+	  } else break;
+	}
+	
+	if ( active[b].empty() ) --b; 
+	else {
+	  end=false;  
+	  Node w=active[b].top();
+	  active[b].pop();
+	  int newlevel=push(w,active);
+	  if ( excess[w] > 0 ) relabel(w, newlevel, active, level_list, 
+				       left, right, b, k, what_heur);
+	  
+	  ++numrelabel; 
+	  if ( numrelabel >= heur ) {
+	    numrelabel=0;
+	    if ( what_heur ) {
+	      what_heur=0;
+	      heur=heur0;
+	      end=false;
+	    } else {
+	      what_heur=1;
+	      heur=heur1;
+	      b=k; 
 	    }
 	  }
+	} 
+      } 
+    }
+
+
+    void preflowPhase1() {
+      
+      int k=n-2;  //bound on the highest level under n containing a node
+      int b=k;    //bound on the highest level under n of an active node
+      
+      VecStack active(n);
+      level.set(s,0);
+      std::queue<Node> bfs_queue;
+      bfs_queue.push(s);
+	    
+      while (!bfs_queue.empty()) {
+	
+	Node v=bfs_queue.front();	
+	bfs_queue.pop();
+	int l=level[v]+1;
+	      
+	InEdgeIt e;
+	for(G.first(e,v); G.valid(e); G.next(e)) {
+	  if ( (*capacity)[e] == (*flow)[e] ) continue;
+	  Node u=G.tail(e);
+	  if ( level[u] >= n ) { 
+	    bfs_queue.push(u);
+	    level.set(u, l);
+	    if ( excess[u] > 0 ) active[l].push(u);
+	  }
 	}
-
-	//the starting flow
-	OutEdgeIt e;
-	for(G.first(e,s); G.valid(e); G.next(e)) 
-	{
-	  T c=capacity[e];
-	  if ( c == 0 ) continue;
-	  Node w=G.head(e);
-	  if ( level[w] < n ) {	  
-	    if ( excess[w] == 0 && w!=t ) {
-	      next.set(w,active[level[w]]);
-	      active[level[w]]=w;
-	    }
-	    flow.set(e, c); 
-	    excess.set(w, excess[w]+c);
+	
+	OutEdgeIt f;
+	for(G.first(f,v); G.valid(f); G.next(f)) {
+	  if ( 0 == (*flow)[f] ) continue;
+	  Node u=G.head(f);
+	  if ( level[u] >= n ) { 
+	    bfs_queue.push(u);
+	    level.set(u, l);
+	    if ( excess[u] > 0 ) active[l].push(u);
 	  }
 	}
       }
-      else 
-      {
-	
-	/*
-	  Reverse_bfs from t in the residual graph, 
-	  to find the starting level.
-	*/
-	level.set(t,0);
-	std::queue<Node> bfs_queue;
-	bfs_queue.push(t);
-	
-	while (!bfs_queue.empty()) {
-	  
-	  Node v=bfs_queue.front();	
-	  bfs_queue.pop();
-	  int l=level[v]+1;
-	  
-	  InEdgeIt e;
-	  for(G.first(e,v); G.valid(e); G.next(e)) {
-	    if ( capacity[e] == flow[e] ) continue;
-	    Node w=G.tail(e);
-	    if ( level[w] == n && w != s ) {
-	      bfs_queue.push(w);
-	      Node first=level_list[l];
-	      if ( G.valid(first) ) left.set(first,w);
-	      right.set(w,first);
-	      level_list[l]=w;
-	      level.set(w, l);
-	    }
-	  }
-	    
-	  OutEdgeIt f;
-	  for(G.first(f,v); G.valid(f); G.next(f)) {
-	    if ( 0 == flow[f] ) continue;
-	    Node w=G.head(f);
-	    if ( level[w] == n && w != s ) {
-	      bfs_queue.push(w);
-	      Node first=level_list[l];
-	      if ( G.valid(first) ) left.set(first,w);
-	      right.set(w,first);
-	      level_list[l]=w;
-	      level.set(w, l);
-	    }
-	  }
-	}
-      
-	
-	/*
-	  Counting the excess
-	*/
-	NodeIt v;
-	for(G.first(v); G.valid(v); G.next(v)) {
-	  T exc=0;
+      b=n-2;
 
-	  InEdgeIt e;
-	  for(G.first(e,v); G.valid(e); G.next(e)) exc+=flow[e];
-	  OutEdgeIt f;
-	  for(G.first(f,v); G.valid(f); G.next(f)) exc-=flow[e];
-
-	  excess.set(v,exc);	  
-
-	  //putting the active nodes into the stack
-	  int lev=level[v];
-	  if ( exc > 0 && lev < n ) {
-	    next.set(v,active[lev]);
-	    active[lev]=v;
-	  }
-	}
-
-
-	//the starting flow
-	OutEdgeIt e;
-	for(G.first(e,s); G.valid(e); G.next(e)) 
-	{
-	  T rem=capacity[e]-flow[e];
-	  if ( rem == 0 ) continue;
-	  Node w=G.head(e);
-	  if ( level[w] < n ) {	  
-	    if ( excess[w] == 0 && w!=t ) {
-	      next.set(w,active[level[w]]);
-	      active[level[w]]=w;
-	    }
-	    flow.set(e, capacity[e]); 
-	    excess.set(w, excess[w]+rem);
-	  }
-	}
-	
-	InEdgeIt f;
-	for(G.first(f,s); G.valid(f); G.next(f)) 
-	{
-	  if ( flow[f] == 0 ) continue;
-	  Node w=G.head(f);
-	  if ( level[w] < n ) {	  
-	    if ( excess[w] == 0 && w!=t ) {
-	      next.set(w,active[level[w]]);
-	      active[level[w]]=w;
-	    }
-	    excess.set(w, excess[w]+flow[f]);
-	    flow.set(f, 0); 
-	  }
-	}
-      }
-
-
-
-
-      /* 
-	 End of preprocessing 
-      */
-
-
-
-      /*
-	Push/relabel on the highest level active nodes.
-      */	
       while ( true ) {
 	
-	if ( b == 0 ) {
-	  if ( phase ) break;
-	  
-	  if ( !what_heur && !end && k > 0 ) {
-	    b=k;
-	    end=true;
-	  } else {
-	    phase=1;
-	    level.set(s,0);
-	    std::queue<Node> bfs_queue;
-	    bfs_queue.push(s);
-	    
-	    while (!bfs_queue.empty()) {
-	      
-	      Node v=bfs_queue.front();	
-	      bfs_queue.pop();
-	      int l=level[v]+1;
-	      
-	      InEdgeIt e;
-	      for(G.first(e,v); G.valid(e); G.next(e)) {
-		if ( capacity[e] == flow[e] ) continue;
-		Node u=G.tail(e);
-		if ( level[u] >= n ) { 
-		  bfs_queue.push(u);
-		  level.set(u, l);
-		  if ( excess[u] > 0 ) {
-		    next.set(u,active[l]);
-		    active[l]=u;
-		  }
-		}
-	      }
-	    
-	      OutEdgeIt f;
-	      for(G.first(f,v); G.valid(f); G.next(f)) {
-		if ( 0 == flow[f] ) continue;
-		Node u=G.head(f);
-		if ( level[u] >= n ) { 
-		  bfs_queue.push(u);
-		  level.set(u, l);
-		  if ( excess[u] > 0 ) {
-		    next.set(u,active[l]);
-		    active[l]=u;
-		  }
-		}
-	      }
-	    }
-	    b=n-2;
-	    }
-	    
-	}
-	  
-	  
-	if ( !G.valid(active[b]) ) --b; 
+	if ( b == 0 ) break;
+
+	if ( active[b].empty() ) --b; 
 	else {
-	  end=false;  
+	  Node w=active[b].top();
+	  active[b].pop();
+	  int newlevel=push(w,active);	  
 
-	  Node w=active[b];
-	  active[b]=next[w];
-	  int lev=level[w];
-	  T exc=excess[w];
-	  int newlevel=n;       //bound on the next level of w
-	  
-	  OutEdgeIt e;
-	  for(G.first(e,w); G.valid(e); G.next(e)) {
-	    
-	    if ( flow[e] == capacity[e] ) continue; 
-	    Node v=G.head(e);            
-	    //e=wv	    
-	    
-	    if( lev > level[v] ) {      
-	      /*Push is allowed now*/
-	      
-	      if ( excess[v]==0 && v!=t && v!=s ) {
-		int lev_v=level[v];
-		next.set(v,active[lev_v]);
-		active[lev_v]=v;
-	      }
-	      
-	      T cap=capacity[e];
-	      T flo=flow[e];
-	      T remcap=cap-flo;
-	      
-	      if ( remcap >= exc ) {       
-		/*A nonsaturating push.*/
-		
-		flow.set(e, flo+exc);
-		excess.set(v, excess[v]+exc);
-		exc=0;
-		break; 
-		
-	      } else { 
-		/*A saturating push.*/
-		
-		flow.set(e, cap);
-		excess.set(v, excess[v]+remcap);
-		exc-=remcap;
-	      }
-	    } else if ( newlevel > level[v] ){
-	      newlevel = level[v];
-	    }	    
-	    
-	  } //for out edges wv 
-	
-	
-	if ( exc > 0 ) {	
-	  InEdgeIt e;
-	  for(G.first(e,w); G.valid(e); G.next(e)) {
-	    
-	    if( flow[e] == 0 ) continue; 
-	    Node v=G.tail(e);  
-	    //e=vw
-	    
-	    if( lev > level[v] ) {  
-	      /*Push is allowed now*/
-	      
-	      if ( excess[v]==0 && v!=t && v!=s ) {
-		int lev_v=level[v];
-		next.set(v,active[lev_v]);
-		active[lev_v]=v;
-	      }
-	      
-	      T flo=flow[e];
-	      
-	      if ( flo >= exc ) { 
-		/*A nonsaturating push.*/
-		
-		flow.set(e, flo-exc);
-		excess.set(v, excess[v]+exc);
-		exc=0;
-		break; 
-	      } else {                                               
-		/*A saturating push.*/
-		
-		excess.set(v, excess[v]+flo);
-		exc-=flo;
-		flow.set(e,0);
-	      }  
-	    } else if ( newlevel > level[v] ) {
-	      newlevel = level[v];
-	    }	    
-	  } //for in edges vw
-	  
-	} // if w still has excess after the out edge for cycle
-	
-	excess.set(w, exc);
-	 
-	/*
-	  Relabel
-	*/
-	
-
-	if ( exc > 0 ) {
-	  //now 'lev' is the old level of w
-	
-	  if ( phase ) {
+	  //relabel
+	  if ( excess[w] > 0 ) {
 	    level.set(w,++newlevel);
-	    next.set(w,active[newlevel]);
-	    active[newlevel]=w;
+	    active[newlevel].push(w);
 	    b=newlevel;
-	  } else {
-	    //unlacing starts
-	    Node right_n=right[w];
-	    Node left_n=left[w];
-
-	    if ( G.valid(right_n) ) {
-	      if ( G.valid(left_n) ) {
-		right.set(left_n, right_n);
-		left.set(right_n, left_n);
-	      } else {
-		level_list[lev]=right_n;   
-		left.set(right_n, INVALID);
-	      } 
-	    } else {
-	      if ( G.valid(left_n) ) {
-		right.set(left_n, INVALID);
-	      } else { 
-		level_list[lev]=INVALID;   
-	      } 
-	    } 
-	    //unlacing ends
-		
-	    if ( !G.valid(level_list[lev]) ) {
-	      
-	       //gapping starts
-	      for (int i=lev; i!=k ; ) {
-		Node v=level_list[++i];
-		while ( G.valid(v) ) {
-		  level.set(v,n);
-		  v=right[v];
-		}
-		level_list[i]=INVALID;
-		if ( !what_heur ) active[i]=INVALID;
-	      }	     
-
-	      level.set(w,n);
-	      b=lev-1;
-	      k=b;
-	      //gapping ends
-	    
-	    } else {
-	      
-	      if ( newlevel == n ) level.set(w,n); 
-	      else {
-		level.set(w,++newlevel);
-		next.set(w,active[newlevel]);
-		active[newlevel]=w;
-		if ( what_heur ) b=newlevel;
-		if ( k < newlevel ) ++k;      //now k=newlevel
-		Node first=level_list[newlevel];
-		if ( G.valid(first) ) left.set(first,w);
-		right.set(w,first);
-		left.set(w,INVALID);
-		level_list[newlevel]=w;
-	      }
-	    }
-
-
-	    ++relabel; 
-	    if ( relabel >= heur ) {
-	      relabel=0;
-	      if ( what_heur ) {
-		what_heur=0;
-		heur=heur0;
-		end=false;
-	      } else {
-		what_heur=1;
-		heur=heur1;
-		b=k; 
-	      }
-	    }
-	  } //phase 0
-	  
-	  
-	} // if ( exc > 0 )
-	  
-	
+	  }
 	}  // if stack[b] is nonempty
-	
       } // while(true)
-
-
-      value = excess[t];
-      /*Max flow value.*/
-     
-    } //void run()
-
-
-
-
-
-    /*
-      Returns the maximum value of a flow.
-     */
-
-    T flowValue() {
-      return value;
     }
 
 
-    FlowMap Flow() {
-      return flow;
-      }
+    //Returns the maximum value of a flow.
+    T flowValue() {
+      return excess[t];
+    }
 
-
-    
-    void Flow(FlowMap& _flow ) {
+    //should be used only between preflowPhase0 and preflowPhase1
+    template<typename _CutMap>
+    void actMinCut(_CutMap& M) {
       NodeIt v;
-      for(G.first(v) ; G.valid(v); G.next(v))
-	_flow.set(v,flow[v]);
+      for(G.first(v); G.valid(v); G.next(v)) 
+	if ( level[v] < n ) M.set(v,false);
+	else M.set(v,true);
     }
 
 
@@ -552,7 +275,6 @@
     /*
       Returns the minimum min cut, by a bfs from s in the residual graph.
     */
-   
     template<typename _CutMap>
     void minMinCut(_CutMap& M) {
     
@@ -568,7 +290,7 @@
 	OutEdgeIt e;
 	for(G.first(e,w) ; G.valid(e); G.next(e)) {
 	  Node v=G.head(e);
-	  if (!M[v] && flow[e] < capacity[e] ) {
+	  if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
 	    queue.push(v);
 	    M.set(v, true);
 	  }
@@ -577,7 +299,7 @@
 	InEdgeIt f;
 	for(G.first(f,w) ; G.valid(f); G.next(f)) {
 	  Node v=G.tail(f);
-	  if (!M[v] && flow[f] > 0 ) {
+	  if (!M[v] && (*flow)[f] > 0 ) {
 	    queue.push(v);
 	    M.set(v, true);
 	  }
@@ -594,10 +316,15 @@
     
     template<typename _CutMap>
     void maxMinCut(_CutMap& M) {
-    
+
+      NodeIt v;
+      for(G.first(v) ; G.valid(v); G.next(v)) {
+	M.set(v, true);
+      }
+
       std::queue<Node> queue;
       
-      M.set(t,true);        
+      M.set(t,false);        
       queue.push(t);
 
       while (!queue.empty()) {
@@ -608,56 +335,319 @@
 	InEdgeIt e;
 	for(G.first(e,w) ; G.valid(e); G.next(e)) {
 	  Node v=G.tail(e);
-	  if (!M[v] && flow[e] < capacity[e] ) {
+	  if (M[v] && (*flow)[e] < (*capacity)[e] ) {
 	    queue.push(v);
-	    M.set(v, true);
+	    M.set(v, false);
 	  }
 	}
 	
 	OutEdgeIt f;
 	for(G.first(f,w) ; G.valid(f); G.next(f)) {
 	  Node v=G.head(f);
-	  if (!M[v] && flow[f] > 0 ) {
+	  if (M[v] && (*flow)[f] > 0 ) {
 	    queue.push(v);
-	    M.set(v, true);
+	    M.set(v, false);
 	  }
 	}
       }
-
-      NodeIt v;
-      for(G.first(v) ; G.valid(v); G.next(v)) {
-	M.set(v, !M[v]);
-      }
-
     }
 
 
-
     template<typename CutMap>
     void minCut(CutMap& M) {
       minMinCut(M);
     }
 
     
-    void reset_target (Node _t) {t=_t;}
-    void reset_source (Node _s) {s=_s;}
+    void resetTarget (const Node _t) {t=_t;}
+
+    void resetSource (const Node _s) {s=_s;}
    
-    template<typename _CapMap>   
-    void reset_cap (_CapMap _cap) {capacity=_cap;}
-
-    template<typename _FlowMap>   
-    void reset_cap (_FlowMap _flow, bool _constzero) {
-      flow=_flow;
-      constzero=_constzero;
+    void resetCap (const CapMap& _cap) {
+      capacity=&_cap;
+    }
+    
+    void resetFlow (FlowMap& _flow) {
+      flow=&_flow;
     }
 
 
+  private:
+
+    int push(const Node w, VecStack& active) {
+      
+      int lev=level[w];
+      T exc=excess[w];
+      int newlevel=n;       //bound on the next level of w
+	  
+      OutEdgeIt e;
+      for(G.first(e,w); G.valid(e); G.next(e)) {
+	    
+	if ( (*flow)[e] == (*capacity)[e] ) continue; 
+	Node v=G.head(e);            
+	    
+	if( lev > level[v] ) { //Push is allowed now
+	  
+	  if ( excess[v]==0 && v!=t && v!=s ) {
+	    int lev_v=level[v];
+	    active[lev_v].push(v);
+	  }
+	  
+	  T cap=(*capacity)[e];
+	  T flo=(*flow)[e];
+	  T remcap=cap-flo;
+	  
+	  if ( remcap >= exc ) { //A nonsaturating push.
+	    
+	    flow->set(e, flo+exc);
+	    excess.set(v, excess[v]+exc);
+	    exc=0;
+	    break; 
+	    
+	  } else { //A saturating push.
+	    flow->set(e, cap);
+	    excess.set(v, excess[v]+remcap);
+	    exc-=remcap;
+	  }
+	} else if ( newlevel > level[v] ) newlevel = level[v];
+      } //for out edges wv 
+      
+      if ( exc > 0 ) {	
+	InEdgeIt e;
+	for(G.first(e,w); G.valid(e); G.next(e)) {
+	  
+	  if( (*flow)[e] == 0 ) continue; 
+	  Node v=G.tail(e); 
+	  
+	  if( lev > level[v] ) { //Push is allowed now
+	    
+	    if ( excess[v]==0 && v!=t && v!=s ) {
+	      int lev_v=level[v];
+	      active[lev_v].push(v);
+	    }
+	    
+	    T flo=(*flow)[e];
+	    
+	    if ( flo >= exc ) { //A nonsaturating push.
+	      
+	      flow->set(e, flo-exc);
+	      excess.set(v, excess[v]+exc);
+	      exc=0;
+	      break; 
+	    } else {  //A saturating push.
+	      
+	      excess.set(v, excess[v]+flo);
+	      exc-=flo;
+	      flow->set(e,0);
+	    }  
+	  } else if ( newlevel > level[v] ) newlevel = level[v];
+	} //for in edges vw
+	
+      } // if w still has excess after the out edge for cycle
+      
+      excess.set(w, exc);
+      
+      return newlevel;
+     }
+
+
+    void preflowPreproc ( flowEnum fe, VecStack& active, 
+			  VecNode& level_list, NNMap& left, NNMap& right ) {
+
+      std::queue<Node> bfs_queue;
+      
+      switch ( fe ) {
+      case ZERO_FLOW: 
+	{
+	  //Reverse_bfs from t, to find the starting level.
+	  level.set(t,0);
+	  bfs_queue.push(t);
+	
+	  while (!bfs_queue.empty()) {
+	    
+	    Node v=bfs_queue.front();	
+	    bfs_queue.pop();
+	    int l=level[v]+1;
+	    
+	    InEdgeIt e;
+	    for(G.first(e,v); G.valid(e); G.next(e)) {
+	      Node w=G.tail(e);
+	      if ( level[w] == n && w != s ) {
+		bfs_queue.push(w);
+		Node first=level_list[l];
+		if ( G.valid(first) ) left.set(first,w);
+		right.set(w,first);
+		level_list[l]=w;
+		level.set(w, l);
+	      }
+	    }
+	  }
+	  
+	  //the starting flow
+	  OutEdgeIt e;
+	  for(G.first(e,s); G.valid(e); G.next(e)) 
+	    {
+	      T c=(*capacity)[e];
+	      if ( c == 0 ) continue;
+	      Node w=G.head(e);
+	      if ( level[w] < n ) {	  
+		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
+		flow->set(e, c); 
+		excess.set(w, excess[w]+c);
+	      }
+	    }
+	  break;
+	}
+	
+      case GEN_FLOW:
+      case PREFLOW:
+	{
+	  //Reverse_bfs from t in the residual graph, 
+	  //to find the starting level.
+	  level.set(t,0);
+	  bfs_queue.push(t);
+	  
+	  while (!bfs_queue.empty()) {
+	    
+	    Node v=bfs_queue.front();	
+	    bfs_queue.pop();
+	    int l=level[v]+1;
+	    
+	    InEdgeIt e;
+	    for(G.first(e,v); G.valid(e); G.next(e)) {
+	      if ( (*capacity)[e] == (*flow)[e] ) continue;
+	      Node w=G.tail(e);
+	      if ( level[w] == n && w != s ) {
+		bfs_queue.push(w);
+		Node first=level_list[l];
+		if ( G.valid(first) ) left.set(first,w);
+		right.set(w,first);
+		level_list[l]=w;
+		level.set(w, l);
+	      }
+	    }
+	    
+	    OutEdgeIt f;
+	    for(G.first(f,v); G.valid(f); G.next(f)) {
+	      if ( 0 == (*flow)[f] ) continue;
+	      Node w=G.head(f);
+	      if ( level[w] == n && w != s ) {
+		bfs_queue.push(w);
+		Node first=level_list[l];
+		if ( G.valid(first) ) left.set(first,w);
+		right.set(w,first);
+		level_list[l]=w;
+		level.set(w, l);
+	      }
+	    }
+	  }
+	  
+	  
+	  //the starting flow
+	  OutEdgeIt e;
+	  for(G.first(e,s); G.valid(e); G.next(e)) 
+	    {
+	      T rem=(*capacity)[e]-(*flow)[e];
+	      if ( rem == 0 ) continue;
+	      Node w=G.head(e);
+	      if ( level[w] < n ) {	  
+		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
+		flow->set(e, (*capacity)[e]); 
+		excess.set(w, excess[w]+rem);
+	      }
+	    }
+	  
+	  InEdgeIt f;
+	  for(G.first(f,s); G.valid(f); G.next(f)) 
+	    {
+	      if ( (*flow)[f] == 0 ) continue;
+	      Node w=G.tail(f);
+	      if ( level[w] < n ) {	  
+		if ( excess[w] == 0 && w!=t ) active[level[w]].push(w);
+		excess.set(w, excess[w]+(*flow)[f]);
+		flow->set(f, 0); 
+	      }
+	    }  
+	  break;
+	} //case PREFLOW
+      }
+    } //preflowPreproc
+
+
+
+    void relabel( const Node w, int newlevel, VecStack& active,  
+		  VecNode& level_list, NNMap& left, 
+		  NNMap& right, int& b, int& k, const bool what_heur ) {
+
+      T lev=level[w];	
+      
+      Node right_n=right[w];
+      Node left_n=left[w];
+      
+      //unlacing starts
+      if ( G.valid(right_n) ) {
+	if ( G.valid(left_n) ) {
+	  right.set(left_n, right_n);
+	  left.set(right_n, left_n);
+	} else {
+	  level_list[lev]=right_n;   
+	  left.set(right_n, INVALID);
+	} 
+      } else {
+	if ( G.valid(left_n) ) {
+	  right.set(left_n, INVALID);
+	} else { 
+	  level_list[lev]=INVALID;   
+	} 
+      } 
+      //unlacing ends
+		
+      if ( !G.valid(level_list[lev]) ) {
+	      
+	//gapping starts
+	for (int i=lev; i!=k ; ) {
+	  Node v=level_list[++i];
+	  while ( G.valid(v) ) {
+	    level.set(v,n);
+	    v=right[v];
+	  }
+	  level_list[i]=INVALID;
+	  if ( !what_heur ) {
+	    while ( !active[i].empty() ) {
+	      active[i].pop();    //FIXME: ezt szebben kene
+	    }
+	  }	     
+	}
+	
+	level.set(w,n);
+	b=lev-1;
+	k=b;
+	//gapping ends
+	
+      } else {
+	
+	if ( newlevel == n ) level.set(w,n); 
+	else {
+	  level.set(w,++newlevel);
+	  active[newlevel].push(w);
+	  if ( what_heur ) b=newlevel;
+	  if ( k < newlevel ) ++k;      //now k=newlevel
+	  Node first=level_list[newlevel];
+	  if ( G.valid(first) ) left.set(first,w);
+	  right.set(w,first);
+	  left.set(w,INVALID);
+	  level_list[newlevel]=w;
+	}
+      }
+      
+    } //relabel
+    
 
   };
 
 } //namespace hugo
 
-#endif //PREFLOW_H
+#endif //HUGO_PREFLOW_H
 
 
 
diff -r 5caac2f7829b -r 6b36be4cffa4 src/work/jacint/preflow_res_comp.cc
--- a/src/work/jacint/preflow_res_comp.cc	Tue Apr 27 22:29:11 2004 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,125 +0,0 @@
-/*
-The only difference between preflow.h and preflow_res.h is that the latter
-uses the ResGraphWrapper, while the first does not. (Bfs is implemented by
-hand in both.) This test program runs Preflow and PreflowRes on the same
-graph, tests the result of these implementations and writes the running time
-of them.  */
-#include <iostream>
-
-#include <smart_graph.h>
-#include <dimacs.h>
-#include <preflow.h>
-#include <preflow_res.h>
-#include <time_measure.h>
-
-using namespace hugo;
-
-int main(int, char **) {
- 
-  typedef SmartGraph Graph;
-  
-  typedef Graph::Node Node;
-  typedef Graph::EdgeIt EdgeIt;
-
-  Graph G;
-  Node s, t;
-  Graph::EdgeMap<int> cap(G);
-  readDimacsMaxFlow(std::cin, G, s, t, cap);
-  Timer ts;
-  
-  std::cout <<
-    "\n  In which way are we faster: using ResGraphWrapper or not?"
-	    <<std::endl;
-  std::cout <<
-    "\n  Running preflow.h on a graph with " << 
-    G.nodeNum() << " nodes and " << G.edgeNum() << " edges..."
-	   << std::endl<<std::endl;
-
-
-  
-  Graph::EdgeMap<int> flow(G,0);
-  Preflow<Graph, int> max_flow_test(G, s, t, cap, flow, 1);
-  ts.reset();
-  max_flow_test.run();
-  std::cout << "Elapsed time NOT using the ResGraphWrapper: " << std::endl 
-	    <<ts << std::endl;
-  
-  Graph::NodeMap<bool> mincut(G);
-  max_flow_test.minMinCut(mincut); 
-  int min_min_cut_value=0;
-  EdgeIt e;
-  for(G.first(e); G.valid(e); G.next(e)) {
-    if (mincut[G.tail(e)] && !mincut[G.head(e)]) min_min_cut_value+=cap[e];
-  }
-
-  Graph::NodeMap<bool> cut(G);
-  max_flow_test.minCut(cut); 
-  int min_cut_value=0;
-  for(G.first(e); G.valid(e); G.next(e)) {
-    if (cut[G.tail(e)] && !cut[G.head(e)]) 
-      min_cut_value+=cap[e];
-  }
-
-  Graph::NodeMap<bool> maxcut(G);
-  max_flow_test.maxMinCut(maxcut); 
-  int max_min_cut_value=0;
-  for(G.first(e); G.valid(e); G.next(e)) {
-    if (maxcut[G.tail(e)] && !maxcut[G.head(e)]) 
-      max_min_cut_value+=cap[e];
-      }
-
-  std::cout << "\n Checking the result: " <<std::endl;  
-  std::cout << "Flow value: "<< max_flow_test.flowValue() << std::endl;
-  std::cout << "Min cut value: "<< min_cut_value << std::endl;
-  std::cout << "Min min cut value: "<< min_min_cut_value << std::endl;
-  std::cout << "Max min cut value: "<< max_min_cut_value << 
-    std::endl;
-
-  if ( max_flow_test.flowValue() == min_cut_value &&
-       min_cut_value == min_min_cut_value &&
-       min_min_cut_value == max_min_cut_value )
-    std::cout << "They are equal! " <<std::endl<< std::endl<<"\n";  
-  
-  Graph::EdgeMap<int> flow2(G,0);
-  PreflowRes<Graph, int> max_flow_test2(G, s, t, cap, flow2, 1);
-  ts.reset();
-  max_flow_test2.run();
-  std::cout << "Elapsed time using the ResGraphWrapper: " << std::endl 
-	    << ts << std::endl;
-  
-  Graph::NodeMap<bool> mincut2(G);
-  max_flow_test2.minMinCut(mincut2); 
-  int min_min_cut_value2=0;
-    for(G.first(e); G.valid(e); G.next(e)) {
-    if (mincut2[G.tail(e)] && !mincut2[G.head(e)]) min_min_cut_value2+=cap[e];
-  }
-
-  Graph::NodeMap<bool> cut2(G);
-  max_flow_test2.minCut(cut2); 
-  int min_cut_value2=0;
-  for(G.first(e); G.valid(e); G.next(e)) {
-    if (cut2[G.tail(e)] && !cut2[G.head(e)]) 
-      min_cut_value2+=cap[e];
-  }
-
-  Graph::NodeMap<bool> maxcut2(G);
-  max_flow_test2.maxMinCut(maxcut2); 
-  int max_min_cut_value2=0;
-  for(G.first(e); G.valid(e); G.next(e)) {
-    if (maxcut2[G.tail(e)] && !maxcut2[G.head(e)]) 
-      max_min_cut_value2+=cap[e];
-      }
-  
-  std::cout << "\n Checking the result: " <<std::endl;  
-  std::cout << "Flow value: "<< max_flow_test2.flowValue() << std::endl;
-  std::cout << "Min cut value: "<< min_cut_value2 << std::endl;
-  std::cout << "Min min cut value: "<< min_min_cut_value2 << std::endl;
-  std::cout << "Max min cut value: "<< max_min_cut_value2 << 
-    std::endl;  
-  if ( max_flow_test.flowValue() == min_cut_value &&
-       min_cut_value == min_min_cut_value &&
-       min_min_cut_value == max_min_cut_value )
-    std::cout << "They are equal! " <<std::endl<<"/n";  
-  
-  return 0;
-}