Changeset 372:e6a156fc186d in lemon-0.x for src

Timestamp:

04/22/04 16:11:28 (20 years ago)

Author:

jacint

Branch:

default

Phase:

public

Convert:

svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@502

Message:

 

Location:

src/work/jacint

Files:

• dijkstra.cc (modified) (1 diff)
• dijkstra.h (modified) (3 diffs)
• f9 (added)
• makefile (modified) (1 diff)
• preflow.cc (modified) (4 diffs)
• preflow.h (modified) (12 diffs)
• preflowproba.h (modified) (3 diffs)

src/work/jacint/dijkstra.cc

@@ -2 +2 @@
 #include <fstream>
 
-#include <smart_graph.h>
-#include <list_graph.h>
-#include <dimacs.h>
+#include <list_graph.hh>
+#include <dimacs.hh>
 #include <dijkstra.h>
 #include <time_measure.h>
-
-#include <bin_heap.h>
-#include <fib_heap.h>
+#include <bin_heap.hh>
 
 using namespace hugo;
 
 int main(int, char **) {
-  typedef SmartGraph::Node Node;
-  typedef SmartGraph::NodeIt NodeIt;
-  typedef SmartGraph::InEdgeIt InEdgeIt; 
+  
+  typedef ListGraph Graph;
 
-  SmartGraph G;
+  typedef Graph::Node Node;
+  typedef Graph::EdgeIt EdgeIt;
+
+  Graph G;
   Node s, t;
-  SmartGraph::EdgeMap<int> cap(G);
+  Graph::EdgeMap<int> cap(G);
   readDimacsMaxFlow(std::cin, G, s, t, cap);
-
-  std::cout << "dijkstra demo ..." << std::endl;
+  Timer ts;
   
-  double pre_time=currTime();
-  Dijkstra<SmartGraph, int, FibHeap<SmartGraph::Node, int, 
-    SmartGraph::NodeMap<int> > > dijkstra_test(G, cap); 
-  dijkstra_test.run(s);
+  std::cout << "Testing dijkstra.h with Fibonacci-heap 
+implementation fib_heap.h ..." << std::endl;
+  
+  Dijkstra<Graph, int
+      //      , BinHeap<ListGraph::NodeIt, int, ListGraph::NodeMap<int> > 
+> dijkstra_test(G, s, cap);
+   ts.reset();
+    dijkstra_test.run();
+    std::cout << "elapsed time: " << ts << std::endl;
   double post_time=currTime();
     
-    std::cout << "running time with fib_heap: " 
-            << post_time-pre_time << " sec"<< std::endl; 
+  std::cout << "running time: " << post_time-pre_time << " sec"<< std::endl; 
  
-  pre_time=currTime();
-  Dijkstra<SmartGraph, int, BinHeap<SmartGraph::Node, int, 
-    SmartGraph::NodeMap<int> > > dijkstra_test2(G, cap);
-  dijkstra_test2.run(s);
-  post_time=currTime();
-  
-  std::cout << "running time with bin_heap: " 
-            << post_time-pre_time << " sec"<< std::endl; 
-  
+  EachEdgeIt e;
 
-  int hiba_fib=0;
-  int hiba_bin=0;
-  NodeIt u;
-  for ( G.first(u) ; G.valid(u); G.next(u) ) {
-    InEdgeIt e;
-    for ( G.first(e,u); G.valid(e); G.next(e) ) {
-      Node v=G.tail(e);
-      if ( dijkstra_test.dist(u) - dijkstra_test.dist(v) > cap[e] )
-        {
-          std::cout<<"Hibas el a fibonaccis Dijkstraban: " 
-                   << dijkstra_test.dist(u) - dijkstra_test.dist(v) - 
-            cap[e]<<std::endl;
-          ++hiba_fib;
-        }
-      if ( dijkstra_test2.dist(u) - dijkstra_test2.dist(v) > cap[e] )
-        {
-          std::cout<<"Hibas el a binarisos Dijkstraban: " 
-                   << dijkstra_test2.dist(u) - dijkstra_test2.dist(v) - 
-            cap[e]<<std::endl;
-          ++hiba_bin;
-        }
-      if ( e==dijkstra_test.pred(u) && 
-           dijkstra_test.dist(u) - dijkstra_test.dist(v) != cap[e] )
-        {
-          std::cout<<"Hibas fael a fibonaccis Dijkstraban: "<<
-            dijkstra_test.dist(u) - dijkstra_test.dist(v)- cap[e]<<std::endl;
-          ++hiba_fib;
-        }
-      if ( e==dijkstra_test2.pred(u) && 
-           dijkstra_test2.dist(u) - dijkstra_test2.dist(v) != cap[e] )
-        {
-          std::cout<<"Hibas fael a binarisos Dijkstraban: "<<
-            dijkstra_test2.dist(u) - dijkstra_test2.dist(v)- cap[e]<<std::endl;
-          ++hiba_bin;
-        }
+  int hiba=0;
+
+  int edge=0;
+
+  for ( G.getFirst(e) ; G.valid(e); G.next(e) ) {
+    NodeIt u=G.tail(e);
+    NodeIt v=G.head(e);
+    ++edge;
+    if ( dijkstra_test.dist(v) - dijkstra_test.dist(u) > cap.get(e) ) {
+      std::cout<<"Hiba: "<<edge<<": "<<dijkstra_test.dist(v) - dijkstra_test.dist(u) - cap.get(e)<<std::endl;
+      ++hiba;
     }
- 
-    if ( dijkstra_test.dist(u) != dijkstra_test2.dist(u) ) 
-      std::cout << "Nem egyezik meg a tavolsag!"<<std::endl;
+  }
 
-
- }
-
-  std::cout << "Hibas elek szama a fibonaccis Dijkstraban: " 
-            << hiba_fib << " a " << G.edgeNum() <<"-bol."<< std::endl;
-  
-  std::cout << "Hibas elek szama a binarisos Dijkstraban: " 
-            << hiba_bin << " a " << G.edgeNum() <<"-bol."<< std::endl;
-  
-
-
+  std::cout << "Osszhibas el: " << hiba << " osszel: " << G.edgeNum() << std::endl;
 
   return 0;

src/work/jacint/dijkstra.h

@@ -1 +1 @@
 // -*- C++ -*-
+
 /* 
  *template <Graph, T, Heap=FibHeap, LengthMap=Graph::EdgeMap<T> >
@@ -27 +28 @@
 #define HUGO_DIJKSTRA_H
 
+///\file
+///\brief Dijkstra algorithm.
+
 #include <fib_heap.h>
+#include <bin_heap.h>
 #include <invalid.h>
 
 namespace hugo {
   
-  template <typename Graph, typename T, 
-    typename Heap=FibHeap<typename Graph::Node, T, 
-    typename Graph::NodeMap<int> >, 
-    typename LengthMap=typename Graph::EdgeMap<T> >
+  //Alpar: Changed the order of the parameters
+  
+  ///%Dijkstra algorithm class.
+
+  ///This class provides an efficient implementation of %Dijkstra algorithm.
+  ///The edge lengths are passed to the algorithm using a
+  ///\ref ReadMapSkeleton "readable map",
+  ///so it is easy to change it to any kind of length.
+  ///
+  ///The type of the length is determined by the \c ValueType of the length map.
+  ///
+  ///It is also possible to change the underlying priority heap.
+  ///
+  ///\param Graph The graph type the algorithm runs on.
+  ///\param LengthMap This read-only
+  ///EdgeMap
+  ///determines the
+  ///lengths of the edges. It is read once for each edge, so the map
+  ///may involve in relatively time consuming process to compute the edge
+  ///length if it is necessary. The default map type is
+  ///\ref GraphSkeleton::EdgeMap "Graph::EdgeMap<int>"
+  ///\param Heap The heap type used by the %Dijkstra
+  ///algorithm. The default
+  ///is using \ref BinHeap "binary heap".
+  
+#ifdef DOXYGEN
+  template <typename Graph,
+            typename LengthMap,
+            typename Heap>
+#else
+  template <typename Graph,
+            typename LengthMap=typename Graph::EdgeMap<int>,
+            template <class,class,class> class Heap = BinHeap >
+#endif
   class Dijkstra{
+  public:
     typedef typename Graph::Node Node;
     typedef typename Graph::NodeIt NodeIt;
@@ -42 +78 @@
     typedef typename Graph::OutEdgeIt OutEdgeIt;
     
+    typedef typename LengthMap::ValueType ValueType;
+    typedef typename Graph::NodeMap<Edge> PredMap;
+    typedef typename Graph::NodeMap<Node> PredNodeMap;
+    typedef typename Graph::NodeMap<ValueType> DistMap;
+
+  private:
     const Graph& G;
     const LengthMap& length;
-    typename Graph::NodeMap<Edge> predecessor;
-    typename Graph::NodeMap<T> distance;
-    //FIXME:
-    typename Graph::NodeMap<bool> reach;
-    //typename Graph::NodeMap<int> reach;
+    PredMap predecessor;
+    PredNodeMap pred_node;
+    DistMap distance;
     
   public :
     
-    /*
-      The distance of the nodes is 0.
-    */
-    Dijkstra(Graph& _G, LengthMap& _length) : G(_G), 
-      length(_length), predecessor(_G), distance(_G), reach(_G) { }
-    
-
-    void run(Node s) {
-      
-      NodeIt u;
-      for ( G.first(u) ; G.valid(u) ; G.next(u) ) {
-        predecessor.set(u,INVALID);
-        distance.set(u,0);
-        reach.set(u,false);
-      }
-     
-      //FIXME:
-      typename Graph::NodeMap<bool> scanned(G,false);
-      //typename Graph::NodeMap<int> scanned(G,false);
-      typename Graph::NodeMap<int> heap_map(G,-1);
-      
-      Heap heap(heap_map);
-
-      heap.push(s,0); 
-      reach.set(s, true);
-
+    Dijkstra(Graph& _G, LengthMap& _length) :
+      G(_G), length(_length), predecessor(_G), pred_node(_G), distance(_G) { }
+    
+    void run(Node s);
+    
+    ///The distance of a node from the source.
+
+    ///Returns the distance of a node from the source.
+    ///\pre \ref run() must be called before using this function.
+    ///\warning If node \c v in unreachable from the source the return value
+    ///of this funcion is undefined.
+    ValueType dist(Node v) const { return distance[v]; }
+    ///Returns the edges of the shortest path tree.
+
+    ///For a node \c v it returns the last edge of the shortest path
+    ///from the source to \c v or INVALID if \c v is unreachable
+    ///from the source.
+    ///\pre \ref run() must be called before using this function.
+    Edge pred(Node v) const { return predecessor[v]; }
+    ///Returns the nodes of the shortest paths.
+
+    ///For a node \c v it returns the last but one node of the shortest path
+    ///from the source to \c v or INVALID if \c v is unreachable
+    ///from the source.
+    ///\pre \ref run() must be called before using this function.
+    Node predNode(Node v) const { return pred_node[v]; }
+    
+    ///Returns a reference to the NodeMap of distances.
+
+    ///\pre \ref run() must be called before using this function.
+    ///
+    const DistMap &distMap() const { return distance;}
+    ///Returns a reference to the shortest path tree map.
+
+    ///Returns a reference to the NodeMap of the edges of the
+    ///shortest path tree.
+    ///\pre \ref run() must be called before using this function.
+    const PredMap &predMap() const { return predecessor;}
+    ///Returns a reference to the map of nodes of  shortest paths.
+
+    ///Returns a reference to the NodeMap of the last but one nodes of the
+    ///shortest paths.
+    ///\pre \ref run() must be called before using this function.
+    const PredNodeMap &predNodeMap() const { return pred_node;}
+
+    ///Checks if a node is reachable from the source.
+
+    ///Returns \c true if \c v is reachable from the source.
+    ///\warning the source node is reported to be unreached!
+    ///\todo Is this what we want?
+    ///\pre \ref run() must be called before using this function.
+    ///
+    bool reached(Node v) { return G.valid(predecessor[v]); }
+    
+  };
+  
+
+  // **********************************************************************
+  //  IMPLEMENTATIONS
+  // **********************************************************************
+
+  ///Runs %Dijkstra algorithm from node the source.
+
+  ///This method runs the %Dijkstra algorithm from a source node \c s
+  ///in order to
+  ///compute the
+  ///shortest path to each node. The algorithm computes
+  ///- The shortest path tree.
+  ///- The distance of each node from the source.
+  template <typename Graph, typename LengthMap,
+            template<class,class,class> class Heap >
+  void Dijkstra<Graph,LengthMap,Heap>::run(Node s) {
+    
+    NodeIt u;
+    for ( G.first(u) ; G.valid(u) ; G.next(u) ) {
+      predecessor.set(u,INVALID);
+      pred_node.set(u,INVALID);
+      // If a node is unreacheable, then why should be the dist=0?
+      // distance.set(u,0);
+      //      reach.set(u,false);
+    }
+    
+    typename Graph::NodeMap<int> heap_map(G,-1);
+    
+    Heap<Node,ValueType,typename Graph::NodeMap<int> > heap(heap_map);
+    
+    heap.push(s,0); 
+    
       while ( !heap.empty() ) {
         
         Node v=heap.top(); 
-        T oldvalue=heap.get(v);
+        ValueType oldvalue=heap[v];
         heap.pop();
         distance.set(v, oldvalue);
-        scanned.set(v,true);
-
-        OutEdgeIt e;
-        for( G.first(e,v); G.valid(e); G.next(e)) {
+        
+        { //FIXME this bracket is for e to be local
+          OutEdgeIt e;
+        for(G.first(e, v);
+            G.valid(e); G.next(e)) {
           Node w=G.head(e); 
-            
-          if ( !scanned[w] ) {
-            if ( !reach[w] ) {
-              reach.set(w,true);
-              heap.push(w,oldvalue+length[e]); 
+          
+          switch(heap.state(w)) {
+          case heap.PRE_HEAP:
+            heap.push(w,oldvalue+length[e]); 
+            predecessor.set(w,e);
+            pred_node.set(w,v);
+            break;
+          case heap.IN_HEAP:
+            if ( oldvalue+length[e] < heap[w] ) {
+              heap.decrease(w, oldvalue+length[e]); 
               predecessor.set(w,e);
-            } else if ( oldvalue+length[e] < heap.get(w) ) {
-              predecessor.set(w,e);
-              heap.decrease(w, oldvalue+length[e]); 
+              pred_node.set(w,v);
             }
+            break;
+          case heap.POST_HEAP:
+            break;
           }
         }
+      } //FIXME tis bracket
       }
-    }
-    
-    T dist(Node v) {
-      return distance[v];
-    }
-
-    Edge pred(Node v) {
-      return predecessor[v];
-    }
-     
-    bool reached(Node v) {
-      return reach[v];
-    }
-    
-  };
-  
-}
+  }
+  
+} //END OF NAMESPACE HUGO
 
 #endif

src/work/jacint/makefile

@@ -4 +4 @@
 CC=$(CXX)
 INCLUDEDIRS ?= -I../../include -I.. -I../{marci,jacint,alpar,klao,akos,athos} 
-CXXFLAGS = -W -Wall -ansi -pedantic -O3 $(INCLUDEDIRS) 
+CXXFLAGS = -W -Wall -ansi -pedantic -O3 $(INCLUDEDIRS)
 LEDAROOT ?= /ledasrc/LEDA-4.1
 
-BINARIES = preflow #dijkstra prim
+BINARIES = preflow #dijkstra_bin_fib_test #preflow # prim
 
 all: $(BINARIES)

src/work/jacint/preflow.cc

@@ -1 +1 @@
 #include <iostream>
-#include <fstream>
 
 #include <smart_graph.h>
 #include <list_graph.h>
 #include <dimacs.h>
-#include <preflow.h>
+#include <preflowproba.h>
 #include <time_measure.h>
 
 using namespace hugo;
 
-// Use a DIMACS max flow file as stdin.
-// read_dimacs_demo < dimacs_max_flow_file
 int main(int, char **) {
-  typedef SmartGraph::Node Node;
-  typedef SmartGraph::EdgeIt EdgeIt;
+ 
+  typedef SmartGraph Graph;
+  
+  typedef Graph::Node Node;
+  typedef Graph::EdgeIt EdgeIt;
 
-  SmartGraph G;
+  Graph G;
   Node s, t;
-  SmartGraph::EdgeMap<int> cap(G);
+  Graph::EdgeMap<int> cap(G);
   readDimacsMaxFlow(std::cin, G, s, t, cap);
-
+  Timer ts;
+  
   std::cout << "preflow demo ..." << std::endl;
   
-  double mintime=1000000;
-
-  for ( int i=1; i!=11; ++i ) {
-    SmartGraph::EdgeMap<int> flow(G);
-    double pre_time=currTime();
-    Preflow<SmartGraph, int> max_flow_test(G, s, t, cap, flow);
-    max_flow_test.run();
-    double post_time=currTime();
-    if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
-  }
-
-  SmartGraph::EdgeMap<int> flow(G);
-  Preflow<SmartGraph, int> max_flow_test(G, s, t, cap, flow);
+  Graph::EdgeMap<int> flow(G);
+  Preflow<Graph, int> max_flow_test(G, s, t, cap, flow, 1);
+  ts.reset();
   max_flow_test.run();
+  std::cout << "elapsed time: " << ts << std::endl;
   
-  SmartGraph::NodeMap<bool> cut(G);
+  Graph::NodeMap<bool> cut(G);
   max_flow_test.minCut(cut); 
   int min_cut_value=0;
@@ -46 +38 @@
   }
 
-  SmartGraph::NodeMap<bool> cut1(G);
+  Graph::NodeMap<bool> cut1(G);
   max_flow_test.minMinCut(cut1); 
   int min_min_cut_value=0;
@@ -54 +46 @@
   }
 
-  SmartGraph::NodeMap<bool> cut2(G);
+  Graph::NodeMap<bool> cut2(G);
   max_flow_test.maxMinCut(cut2); 
   int max_min_cut_value=0;
@@ -62 +54 @@
       }
   
-  std::cout << "min time of 10 runs: " << mintime << " sec"<< std::endl; 
   std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
   std::cout << "min cut value: "<< min_cut_value << std::endl;

src/work/jacint/preflow.h

@@ -1 +1 @@
 // -*- 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: 
@@ -13 +19 @@
 Constructors:
 
-Preflow(Graph, Node, Node, CapMap, FlowMap)
+Preflow(Graph, Node, Node, CapMap, FlowMap, bool) : bool must be false if 
+     FlowMap is not constant zero, and should be true if it is
 
 Members:
@@ -30 +37 @@
      a min cut. M should be a map of bools initialized to false.
 
+FIXME reset
+
 */
 
@@ -45 +54 @@
   template <typename Graph, typename T, 
             typename CapMap=typename Graph::EdgeMap<T>, 
-            typename FlowMap=typename Graph::EdgeMap<T> >
+            typename FlowMap=typename Graph::EdgeMap<T> >
   class Preflow {
     
@@ -60 +69 @@
     FlowMap& flow;
     T value;
+    bool constzero;
 
   public:
-    Preflow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity, 
-            FlowMap& _flow ) :
-      G(_G), s(_s), t(_t), capacity(_capacity), flow(_flow) {}
-
+    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) {}
+    
+    
     void run() {
+      
+      value=0;                //for the subsequent runs
 
       bool phase=0;        //phase 0 is the 1st phase, phase 1 is the 2nd
@@ -90 +103 @@
       typename Graph::NodeMap<T> excess(G); 
 
-      std::vector<Node> active(n,INVALID);
+      std::vector<Node> active(n-1,INVALID);
       typename Graph::NodeMap<Node> next(G,INVALID);
       //Stack of the active nodes in level i < n.
@@ -102 +115 @@
       */
 
-      /*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()) {
-
-        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);
-          }
-        }
-      }
-      
-      level.set(s,n);
-      
-
-      /* Starting flow. It is everywhere 0 at the moment. */     
-      OutEdgeIt e;
-      for(G.first(e,s); G.valid(e); G.next(e)) 
+
+      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()) {
+          
+          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];
@@ -146 +159 @@
           }
         }
+      }
+      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;
+
+          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); 
+          }
+        }
+      }
+
+
+
 
       /* 
@@ -339 +458 @@
             //unlacing ends
                 
-            //gapping starts
             if ( !G.valid(level_list[lev]) ) {
               
+               //gapping starts
               for (int i=lev; i!=k ; ) {
                 Node v=level_list[++i];
@@ -356 +475 @@
               k=b;
               //gapping ends
+            
             } else {
               
@@ -364 +484 @@
                 active[newlevel]=w;
                 if ( what_heur ) b=newlevel;
-                if ( k < newlevel ) ++k;
+                if ( k < newlevel ) ++k;      //now k=newlevel
                 Node first=level_list[newlevel];
                 if ( G.valid(first) ) left.set(first,w);
@@ -519 +639 @@
     }
 
+    
+    void reset_target (Node _t) {t=_t;}
+    void reset_source (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;
+    }
+
+
 
   };

src/work/jacint/preflowproba.h

@@ -303 +303 @@
               int l=level[v]+1;
               
-              ResInEdgeIt e;
+              /*              ResInEdgeIt e;
               for(res_graph.first(e,s); res_graph.valid(e); 
                   res_graph.next(e)) {
@@ -315 +315 @@
                   }
                 }
-              }
-              /*              InEdgeIt e;
+                }*/
+                    InEdgeIt e;
               for(G.first(e,v); G.valid(e); G.next(e)) {
                 if ( capacity[e] == flow[e] ) continue;
@@ -342 +342 @@
                   }
                 }
-                }*/
+                }
             }
             b=n-2;