Index: INSTALL
===================================================================
--- INSTALL	(revision 890)
+++ INSTALL	(revision 615)
@@ -174,24 +174,2 @@
 
    Disable COIN-OR support.
-
-
-Makefile Variables
-==================
-
-Some Makefile variables are reserved by the GNU Coding Standards for
-the use of the "user" - the person building the package. For instance,
-CXX and CXXFLAGS are such variables, and have the same meaning as
-explained in the previous section. These variables can be set on the
-command line when invoking `make' like this:
-`make [VARIABLE=VALUE]...'
-
-WARNINGCXXFLAGS is a non-standard Makefile variable introduced by us
-to hold several compiler flags related to warnings. Its default value
-can be overridden when invoking `make'. For example to disable all
-warning flags use `make WARNINGCXXFLAGS='.
-
-In order to turn off a single flag from the default set of warning
-flags, you can use the CXXFLAGS variable, since this is passed after
-WARNINGCXXFLAGS. For example to turn off `-Wold-style-cast' (which is
-used by default when g++ is detected) you can use
-`make CXXFLAGS="-g -O2 -Wno-old-style-cast"'.
Index: doc/CMakeLists.txt
===================================================================
--- doc/CMakeLists.txt	(revision 895)
+++ doc/CMakeLists.txt	(revision 791)
@@ -27,5 +27,4 @@
     COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/nodeshape_3.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_3.eps
     COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/nodeshape_4.png ${CMAKE_CURRENT_SOURCE_DIR}/images/nodeshape_4.eps
-    COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/planar.png ${CMAKE_CURRENT_SOURCE_DIR}/images/planar.eps
     COMMAND ${GHOSTSCRIPT_EXECUTABLE} ${GHOSTSCRIPT_OPTIONS} -r18 -sOutputFile=gen-images/strongly_connected_components.png ${CMAKE_CURRENT_SOURCE_DIR}/images/strongly_connected_components.eps
     COMMAND ${CMAKE_COMMAND} -E remove_directory html
Index: doc/Makefile.am
===================================================================
--- doc/Makefile.am	(revision 895)
+++ doc/Makefile.am	(revision 791)
@@ -29,5 +29,4 @@
 	edge_biconnected_components.eps \
 	node_biconnected_components.eps \
-	planar.eps \
 	strongly_connected_components.eps
 
Index: doc/groups.dox
===================================================================
--- doc/groups.dox	(revision 879)
+++ doc/groups.dox	(revision 818)
@@ -407,11 +407,13 @@
  - \ref NetworkSimplex Primal Network Simplex algorithm with various
    pivot strategies \ref dantzig63linearprog, \ref kellyoneill91netsimplex.
- - \ref CostScaling Cost Scaling algorithm based on push/augment and
-   relabel operations \ref goldberg90approximation, \ref goldberg97efficient,
+ - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
+   cost scaling \ref goldberg90approximation, \ref goldberg97efficient,
    \ref bunnagel98efficient.
- - \ref CapacityScaling Capacity Scaling algorithm based on the successive
-   shortest path method \ref edmondskarp72theoretical.
- - \ref CycleCanceling Cycle-Canceling algorithms, two of which are
-   strongly polynomial \ref klein67primal, \ref goldberg89cyclecanceling.
+ - \ref CapacityScaling Successive Shortest %Path algorithm with optional
+   capacity scaling \ref edmondskarp72theoretical.
+ - \ref CancelAndTighten The Cancel and Tighten algorithm
+   \ref goldberg89cyclecanceling.
+ - \ref CycleCanceling Cycle-Canceling algorithms
+   \ref klein67primal, \ref goldberg89cyclecanceling.
 
 In general NetworkSimplex is the most efficient implementation,
Index: c/images/planar.eps
===================================================================
--- doc/images/planar.eps	(revision 895)
+++ 	(revision )
@@ -1,181 +1,0 @@
-%!PS-Adobe-2.0 EPSF-2.0
-%%Creator: LEMON, graphToEps()
-%%CreationDate: Fri Oct 19 18:32:32 2007
-%%BoundingBox: 0 0 596 842
-%%DocumentPaperSizes: a4
-%%EndComments
-/lb { setlinewidth setrgbcolor newpath moveto
-      4 2 roll 1 index 1 index curveto stroke } bind def
-/l { setlinewidth setrgbcolor newpath moveto lineto stroke } bind def
-/c { newpath dup 3 index add 2 index moveto 0 360 arc closepath } bind def
-/sq { newpath 2 index 1 index add 2 index 2 index add moveto
-      2 index 1 index sub 2 index 2 index add lineto
-      2 index 1 index sub 2 index 2 index sub lineto
-      2 index 1 index add 2 index 2 index sub lineto
-      closepath pop pop pop} bind def
-/di { newpath 2 index 1 index add 2 index moveto
-      2 index             2 index 2 index add lineto
-      2 index 1 index sub 2 index             lineto
-      2 index             2 index 2 index sub lineto
-      closepath pop pop pop} bind def
-/nc { 0 0 0 setrgbcolor 5 index 5 index 5 index c fill
-     setrgbcolor 1.1 div c fill
-   } bind def
-/nsq { 0 0 0 setrgbcolor 5 index 5 index 5 index sq fill
-     setrgbcolor 1.1 div sq fill
-   } bind def
-/ndi { 0 0 0 setrgbcolor 5 index 5 index 5 index di fill
-     setrgbcolor 1.1 div di fill
-   } bind def
-/nfemale { 0 0 0 setrgbcolor 3 index 0.0909091 1.5 mul mul setlinewidth
-  newpath 5 index 5 index moveto 5 index 5 index 5 index 3.01 mul sub
-  lineto 5 index 4 index .7 mul sub 5 index 5 index 2.2 mul sub moveto
-  5 index 4 index .7 mul add 5 index 5 index 2.2 mul sub lineto stroke
-  5 index 5 index 5 index c fill
-  setrgbcolor 1.1 div c fill
-  } bind def
-/nmale {
-  0 0 0 setrgbcolor 3 index 0.0909091 1.5 mul mul setlinewidth
-  newpath 5 index 5 index moveto
-  5 index 4 index 1 mul 1.5 mul add
-  5 index 5 index 3 sqrt 1.5 mul mul add
-  1 index 1 index lineto
-  1 index 1 index 7 index sub moveto
-  1 index 1 index lineto
-  exch 5 index 3 sqrt .5 mul mul sub exch 5 index .5 mul sub lineto
-  stroke
-  5 index 5 index 5 index c fill
-  setrgbcolor 1.1 div c fill
-  } bind def
-/arrl 1 def
-/arrw 0.3 def
-/lrl { 2 index mul exch 2 index mul exch rlineto pop} bind def
-/arr { setrgbcolor /y1 exch def /x1 exch def /dy exch def /dx exch def
-       /w exch def /len exch def
-       newpath x1 dy w 2 div mul add y1 dx w 2 div mul sub moveto
-       len w sub arrl sub dx dy lrl
-       arrw dy dx neg lrl
-       dx arrl w add mul dy w 2 div arrw add mul sub
-       dy arrl w add mul dx w 2 div arrw add mul add rlineto
-       dx arrl w add mul neg dy w 2 div arrw add mul sub
-       dy arrl w add mul neg dx w 2 div arrw add mul add rlineto
-       arrw dy dx neg lrl
-       len w sub arrl sub neg dx dy lrl
-       closepath fill } bind def
-/cshow { 2 index 2 index moveto dup stringwidth pop
-         neg 2 div fosi .35 mul neg rmoveto show pop pop} def
-
-gsave
-15 138.307 translate
-12.7843 dup scale
-90 rotate
-0.608112 -43.6081 translate
-%Edges:
-gsave
-9 31 9.5 30.5 10 30 0 0 0 0.091217 lb
-9 31 5.5 34.5 2 38 0 0 0 0.091217 lb
-9 31 25.5 16 42 1 0 0 0 0.091217 lb
-3 40 23 20.5 43 1 0 0 0 0.091217 lb
-3 40 22.5 20.5 42 1 0 0 0 0.091217 lb
-3 40 2.5 40.5 2 41 0 0 0 0.091217 lb
-13 25 10.5 24.5 8 24 0 0 0 0.091217 lb
-13 25 12 27 11 29 0 0 0 0.091217 lb
-3 4 2.5 3 2 2 0 0 0 0.091217 lb
-3 4 4.5 3 6 2 0 0 0 0.091217 lb
-6 25 7 24.5 8 24 0 0 0 0.091217 lb
-6 25 6 24.5 6 24 0 0 0 0.091217 lb
-34 2 33.5 2 33 2 0 0 0 0.091217 lb
-34 2 35 2 36 2 0 0 0 0.091217 lb
-6 8 16 9 26 10 0 0 0 0.091217 lb
-6 8 6 10.5 6 13 0 0 0 0.091217 lb
-6 8 6 7.5 6 7 0 0 0 0.091217 lb
-26 10 27.5 8.5 29 7 0 0 0 0.091217 lb
-26 10 27.5 9 29 8 0 0 0 0.091217 lb
-10 30 10.5 29.5 11 29 0 0 0 0.091217 lb
-8 24 7 23.5 6 23 0 0 0 0.091217 lb
-8 24 8 24.5 8 25 0 0 0 0.091217 lb
-33 2 32.5 2 32 2 0 0 0 0.091217 lb
-29 7 17.5 7 6 7 0 0 0 0.091217 lb
-2 2 1.5 22 1 42 0 0 0 0.091217 lb
-2 2 3.5 2 5 2 0 0 0 0.091217 lb
-21 15 13.5 14.5 6 14 0 0 0 0.091217 lb
-21 15 21 15.5 21 16 0 0 0 0.091217 lb
-1 42 0.5 42.5 0 43 0 0 0 0.091217 lb
-1 42 1.5 41.5 2 41 0 0 0 0.091217 lb
-6 15 6 15.5 6 16 0 0 0 0.091217 lb
-6 15 6 14.5 6 14 0 0 0 0.091217 lb
-43 1 22 0.5 1 0 0 0 0 0.091217 lb
-31 2 18.5 2 6 2 0 0 0 0.091217 lb
-31 2 31.5 2 32 2 0 0 0 0.091217 lb
-6 24 6 23.5 6 23 0 0 0 0.091217 lb
-6 16 6 16.5 6 17 0 0 0 0.091217 lb
-6 23 6 20 6 17 0 0 0 0.091217 lb
-6 2 5.5 2 5 2 0 0 0 0.091217 lb
-6 2 6 4.5 6 7 0 0 0 0.091217 lb
-0 43 0.5 21.5 1 0 0 0 0 0.091217 lb
-1 1 19.5 1.5 38 2 0 0 0 0.091217 lb
-1 1 1 0.5 1 0 0 0 0 0.091217 lb
-2 38 5.5 31.5 9 25 0 0 0 0.091217 lb
-25 13 15.5 13 6 13 0 0 0 0.091217 lb
-25 13 15.5 13.5 6 14 0 0 0 0.091217 lb
-8 25 8.5 25 9 25 0 0 0 0.091217 lb
-11 29 24.5 15.5 38 2 0 0 0 0.091217 lb
-6 17 11.5 18 17 19 0 0 0 0.091217 lb
-16 23 26.5 12.5 37 2 0 0 0 0.091217 lb
-16 23 18.5 19.5 21 16 0 0 0 0.091217 lb
-36 2 36.5 2 37 2 0 0 0 0.091217 lb
-36 2 32.5 5 29 8 0 0 0 0.091217 lb
-6 13 6 13.5 6 14 0 0 0 0.091217 lb
-37 2 37.5 2 38 2 0 0 0 0.091217 lb
-21 16 19 17.5 17 19 0 0 0 0.091217 lb
-grestore
-%Nodes:
-gsave
-29 8 0.304556 1 1 1 nc
-2 41 0.304556 1 1 1 nc
-6 7 0.304556 1 1 1 nc
-5 2 0.304556 1 1 1 nc
-17 19 0.304556 1 1 1 nc
-21 16 0.304556 1 1 1 nc
-1 0 0.304556 1 1 1 nc
-9 25 0.304556 1 1 1 nc
-6 14 0.304556 1 1 1 nc
-42 1 0.304556 1 1 1 nc
-38 2 0.304556 1 1 1 nc
-37 2 0.304556 1 1 1 nc
-6 13 0.304556 1 1 1 nc
-36 2 0.304556 1 1 1 nc
-16 23 0.304556 1 1 1 nc
-6 17 0.304556 1 1 1 nc
-11 29 0.304556 1 1 1 nc
-8 25 0.304556 1 1 1 nc
-32 2 0.304556 1 1 1 nc
-25 13 0.304556 1 1 1 nc
-2 38 0.304556 1 1 1 nc
-1 1 0.304556 1 1 1 nc
-0 43 0.304556 1 1 1 nc
-6 2 0.304556 1 1 1 nc
-6 23 0.304556 1 1 1 nc
-6 16 0.304556 1 1 1 nc
-6 24 0.304556 1 1 1 nc
-31 2 0.304556 1 1 1 nc
-43 1 0.304556 1 1 1 nc
-6 15 0.304556 1 1 1 nc
-1 42 0.304556 1 1 1 nc
-21 15 0.304556 1 1 1 nc
-2 2 0.304556 1 1 1 nc
-29 7 0.304556 1 1 1 nc
-33 2 0.304556 1 1 1 nc
-8 24 0.304556 1 1 1 nc
-10 30 0.304556 1 1 1 nc
-26 10 0.304556 1 1 1 nc
-6 8 0.304556 1 1 1 nc
-34 2 0.304556 1 1 1 nc
-6 25 0.304556 1 1 1 nc
-3 4 0.304556 1 1 1 nc
-13 25 0.304556 1 1 1 nc
-3 40 0.304556 1 1 1 nc
-9 31 0.304556 1 1 1 nc
-grestore
-grestore
-showpage
Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am	(revision 888)
+++ lemon/Makefile.am	(revision 863)
@@ -63,5 +63,4 @@
 	lemon/binom_heap.h \
 	lemon/bucket_heap.h \
-	lemon/capacity_scaling.h \
 	lemon/cbc.h \
 	lemon/circulation.h \
@@ -70,9 +69,7 @@
 	lemon/concept_check.h \
 	lemon/connectivity.h \
+	lemon/counter.h \
 	lemon/core.h \
-	lemon/cost_scaling.h \
-	lemon/counter.h \
 	lemon/cplex.h \
-	lemon/cycle_canceling.h \
 	lemon/dfs.h \
 	lemon/dijkstra.h \
Index: lemon/bellman_ford.h
===================================================================
--- lemon/bellman_ford.h	(revision 891)
+++ lemon/bellman_ford.h	(revision 835)
@@ -172,9 +172,4 @@
   /// the lengths of the arcs. The default map type is
   /// \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref BellmanFordDefaultTraits
-  /// "BellmanFordDefaultTraits<GR, LEN>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
 #ifdef DOXYGEN
   template <typename GR, typename LEN, typename TR>
@@ -243,7 +238,5 @@
 	_dist = Traits::createDistMap(*_gr);
       }
-      if(!_mask) {
-        _mask = new MaskMap(*_gr);
-      }
+      _mask = new MaskMap(*_gr, false);
     }
     
@@ -411,8 +404,4 @@
 	  _process.push_back(it);
 	  _mask->set(it, true);
-	}
-      } else {
-	for (NodeIt it(*_gr); it != INVALID; ++it) {
-	  _mask->set(it, false);
 	}
       }
@@ -939,7 +928,4 @@
   /// This class should only be used through the \ref bellmanFord()
   /// function, which makes it easier to use the algorithm.
-  ///
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm.
   template<class TR>
   class BellmanFordWizard : public TR {
Index: lemon/bfs.h
===================================================================
--- lemon/bfs.h	(revision 891)
+++ lemon/bfs.h	(revision 835)
@@ -122,9 +122,4 @@
   ///\tparam GR The type of the digraph the algorithm runs on.
   ///The default type is \ref ListDigraph.
-  ///\tparam TR The traits class that defines various types used by the
-  ///algorithm. By default, it is \ref BfsDefaultTraits
-  ///"BfsDefaultTraits<GR>".
-  ///In most cases, this parameter should not be set directly,
-  ///consider to use the named template parameters instead.
 #ifdef DOXYGEN
   template <typename GR,
@@ -963,7 +958,4 @@
   /// This class should only be used through the \ref bfs() function,
   /// which makes it easier to use the algorithm.
-  ///
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm.
   template<class TR>
   class BfsWizard : public TR
@@ -1304,9 +1296,9 @@
   /// does not observe the BFS events. If you want to observe the BFS
   /// events, you should implement your own visitor class.
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref BfsVisitDefaultTraits
-  /// "BfsVisitDefaultTraits<GR>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
+  /// \tparam TR Traits class to set various data types used by the
+  /// algorithm. The default traits class is
+  /// \ref BfsVisitDefaultTraits "BfsVisitDefaultTraits<GR>".
+  /// See \ref BfsVisitDefaultTraits for the documentation of
+  /// a BFS visit traits class.
 #ifdef DOXYGEN
   template <typename GR, typename VS, typename TR>
Index: lemon/bits/map_extender.h
===================================================================
--- lemon/bits/map_extender.h	(revision 867)
+++ lemon/bits/map_extender.h	(revision 765)
@@ -85,34 +85,34 @@
       typedef typename Map::Value Value;
 
-      MapIt() : map(NULL) {}
-
-      MapIt(Invalid i) : Parent(i), map(NULL) {}
-
-      explicit MapIt(Map& _map) : map(&_map) {
-        map->notifier()->first(*this);
+      MapIt() {}
+
+      MapIt(Invalid i) : Parent(i) { }
+
+      explicit MapIt(Map& _map) : map(_map) {
+        map.notifier()->first(*this);
       }
 
       MapIt(const Map& _map, const Item& item)
-        : Parent(item), map(&_map) {}
+        : Parent(item), map(_map) {}
 
       MapIt& operator++() {
-        map->notifier()->next(*this);
+        map.notifier()->next(*this);
         return *this;
       }
 
       typename MapTraits<Map>::ConstReturnValue operator*() const {
-        return (*map)[*this];
+        return map[*this];
       }
 
       typename MapTraits<Map>::ReturnValue operator*() {
-        return (*map)[*this];
+        return map[*this];
       }
 
       void set(const Value& value) {
-        map->set(*this, value);
-      }
-
-    protected:
-      Map* map;
+        map.set(*this, value);
+      }
+
+    protected:
+      Map& map;
 
     };
@@ -125,10 +125,10 @@
       typedef typename Map::Value Value;
 
-      ConstMapIt() : map(NULL) {}
-
-      ConstMapIt(Invalid i) : Parent(i), map(NULL) {}
-
-      explicit ConstMapIt(Map& _map) : map(&_map) {
-        map->notifier()->first(*this);
+      ConstMapIt() {}
+
+      ConstMapIt(Invalid i) : Parent(i) { }
+
+      explicit ConstMapIt(Map& _map) : map(_map) {
+        map.notifier()->first(*this);
       }
 
@@ -137,5 +137,5 @@
 
       ConstMapIt& operator++() {
-        map->notifier()->next(*this);
+        map.notifier()->next(*this);
         return *this;
       }
@@ -146,5 +146,5 @@
 
     protected:
-      const Map* map;
+      const Map& map;
     };
 
@@ -153,23 +153,23 @@
 
     public:
-      ItemIt() : map(NULL) {}
-
-
-      ItemIt(Invalid i) : Parent(i), map(NULL) {}
-
-      explicit ItemIt(Map& _map) : map(&_map) {
-        map->notifier()->first(*this);
+
+      ItemIt() {}
+
+      ItemIt(Invalid i) : Parent(i) { }
+
+      explicit ItemIt(Map& _map) : map(_map) {
+        map.notifier()->first(*this);
       }
 
       ItemIt(const Map& _map, const Item& item)
-        : Parent(item), map(&_map) {}
+        : Parent(item), map(_map) {}
 
       ItemIt& operator++() {
-        map->notifier()->next(*this);
-        return *this;
-      }
-
-    protected:
-      const Map* map;
+        map.notifier()->next(*this);
+        return *this;
+      }
+
+    protected:
+      const Map& map;
 
     };
@@ -232,34 +232,34 @@
       typedef typename Map::Value Value;
 
-      MapIt() : map(NULL) {}
-
-      MapIt(Invalid i) : Parent(i), map(NULL) { }
-
-      explicit MapIt(Map& _map) : map(&_map) {
-        map->graph.first(*this);
+      MapIt() {}
+
+      MapIt(Invalid i) : Parent(i) { }
+
+      explicit MapIt(Map& _map) : map(_map) {
+        map.graph.first(*this);
       }
 
       MapIt(const Map& _map, const Item& item)
-        : Parent(item), map(&_map) {}
+        : Parent(item), map(_map) {}
 
       MapIt& operator++() {
-        map->graph.next(*this);
+        map.graph.next(*this);
         return *this;
       }
 
       typename MapTraits<Map>::ConstReturnValue operator*() const {
-        return (*map)[*this];
+        return map[*this];
       }
 
       typename MapTraits<Map>::ReturnValue operator*() {
-        return (*map)[*this];
+        return map[*this];
       }
 
       void set(const Value& value) {
-        map->set(*this, value);
-      }
-
-    protected:
-      Map* map;
+        map.set(*this, value);
+      }
+
+    protected:
+      Map& map;
 
     };
@@ -272,26 +272,26 @@
       typedef typename Map::Value Value;
 
-      ConstMapIt() : map(NULL) {}
-
-      ConstMapIt(Invalid i) : Parent(i), map(NULL) { }
-
-      explicit ConstMapIt(Map& _map) : map(&_map) {
-        map->graph.first(*this);
+      ConstMapIt() {}
+
+      ConstMapIt(Invalid i) : Parent(i) { }
+
+      explicit ConstMapIt(Map& _map) : map(_map) {
+        map.graph.first(*this);
       }
 
       ConstMapIt(const Map& _map, const Item& item)
-        : Parent(item), map(&_map) {}
+        : Parent(item), map(_map) {}
 
       ConstMapIt& operator++() {
-        map->graph.next(*this);
+        map.graph.next(*this);
         return *this;
       }
 
       typename MapTraits<Map>::ConstReturnValue operator*() const {
-        return (*map)[*this];
-      }
-
-    protected:
-      const Map* map;
+        return map[*this];
+      }
+
+    protected:
+      const Map& map;
     };
 
@@ -300,23 +300,23 @@
 
     public:
-      ItemIt() : map(NULL) {}
-
-
-      ItemIt(Invalid i) : Parent(i), map(NULL) { }
-
-      explicit ItemIt(Map& _map) : map(&_map) {
-        map->graph.first(*this);
+
+      ItemIt() {}
+
+      ItemIt(Invalid i) : Parent(i) { }
+
+      explicit ItemIt(Map& _map) : map(_map) {
+        map.graph.first(*this);
       }
 
       ItemIt(const Map& _map, const Item& item)
-        : Parent(item), map(&_map) {}
+        : Parent(item), map(_map) {}
 
       ItemIt& operator++() {
-        map->graph.next(*this);
-        return *this;
-      }
-
-    protected:
-      const Map* map;
+        map.graph.next(*this);
+        return *this;
+      }
+
+    protected:
+      const Map& map;
 
     };
Index: mon/capacity_scaling.h
===================================================================
--- lemon/capacity_scaling.h	(revision 899)
+++ 	(revision )
@@ -1,985 +1,0 @@
-/* -*- C++ -*-
- *
- * This file is a part of LEMON, a generic C++ optimization library
- *
- * Copyright (C) 2003-2008
- * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
- * (Egervary Research Group on Combinatorial Optimization, EGRES).
- *
- * Permission to use, modify and distribute this software is granted
- * provided that this copyright notice appears in all copies. For
- * precise terms see the accompanying LICENSE file.
- *
- * This software is provided "AS IS" with no warranty of any kind,
- * express or implied, and with no claim as to its suitability for any
- * purpose.
- *
- */
-
-#ifndef LEMON_CAPACITY_SCALING_H
-#define LEMON_CAPACITY_SCALING_H
-
-/// \ingroup min_cost_flow_algs
-///
-/// \file
-/// \brief Capacity Scaling algorithm for finding a minimum cost flow.
-
-#include <vector>
-#include <limits>
-#include <lemon/core.h>
-#include <lemon/bin_heap.h>
-
-namespace lemon {
-
-  /// \brief Default traits class of CapacityScaling algorithm.
-  ///
-  /// Default traits class of CapacityScaling algorithm.
-  /// \tparam GR Digraph type.
-  /// \tparam V The number type used for flow amounts, capacity bounds
-  /// and supply values. By default it is \c int.
-  /// \tparam C The number type used for costs and potentials.
-  /// By default it is the same as \c V.
-  template <typename GR, typename V = int, typename C = V>
-  struct CapacityScalingDefaultTraits
-  {
-    /// The type of the digraph
-    typedef GR Digraph;
-    /// The type of the flow amounts, capacity bounds and supply values
-    typedef V Value;
-    /// The type of the arc costs
-    typedef C Cost;
-
-    /// \brief The type of the heap used for internal Dijkstra computations.
-    ///
-    /// The type of the heap used for internal Dijkstra computations.
-    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
-    /// its priority type must be \c Cost and its cross reference type
-    /// must be \ref RangeMap "RangeMap<int>".
-    typedef BinHeap<Cost, RangeMap<int> > Heap;
-  };
-
-  /// \addtogroup min_cost_flow_algs
-  /// @{
-
-  /// \brief Implementation of the Capacity Scaling algorithm for
-  /// finding a \ref min_cost_flow "minimum cost flow".
-  ///
-  /// \ref CapacityScaling implements the capacity scaling version
-  /// of the successive shortest path algorithm for finding a
-  /// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,
-  /// \ref edmondskarp72theoretical. It is an efficient dual
-  /// solution method.
-  ///
-  /// Most of the parameters of the problem (except for the digraph)
-  /// can be given using separate functions, and the algorithm can be
-  /// executed using the \ref run() function. If some parameters are not
-  /// specified, then default values will be used.
-  ///
-  /// \tparam GR The digraph type the algorithm runs on.
-  /// \tparam V The number type used for flow amounts, capacity bounds
-  /// and supply values in the algorithm. By default, it is \c int.
-  /// \tparam C The number type used for costs and potentials in the
-  /// algorithm. By default, it is the same as \c V.
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref CapacityScalingDefaultTraits
-  /// "CapacityScalingDefaultTraits<GR, V, C>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
-  ///
-  /// \warning Both number types must be signed and all input data must
-  /// be integer.
-  /// \warning This algorithm does not support negative costs for such
-  /// arcs that have infinite upper bound.
-#ifdef DOXYGEN
-  template <typename GR, typename V, typename C, typename TR>
-#else
-  template < typename GR, typename V = int, typename C = V,
-             typename TR = CapacityScalingDefaultTraits<GR, V, C> >
-#endif
-  class CapacityScaling
-  {
-  public:
-
-    /// The type of the digraph
-    typedef typename TR::Digraph Digraph;
-    /// The type of the flow amounts, capacity bounds and supply values
-    typedef typename TR::Value Value;
-    /// The type of the arc costs
-    typedef typename TR::Cost Cost;
-
-    /// The type of the heap used for internal Dijkstra computations
-    typedef typename TR::Heap Heap;
-
-    /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
-    typedef TR Traits;
-
-  public:
-
-    /// \brief Problem type constants for the \c run() function.
-    ///
-    /// Enum type containing the problem type constants that can be
-    /// returned by the \ref run() function of the algorithm.
-    enum ProblemType {
-      /// The problem has no feasible solution (flow).
-      INFEASIBLE,
-      /// The problem has optimal solution (i.e. it is feasible and
-      /// bounded), and the algorithm has found optimal flow and node
-      /// potentials (primal and dual solutions).
-      OPTIMAL,
-      /// The digraph contains an arc of negative cost and infinite
-      /// upper bound. It means that the objective function is unbounded
-      /// on that arc, however, note that it could actually be bounded
-      /// over the feasible flows, but this algroithm cannot handle
-      /// these cases.
-      UNBOUNDED
-    };
-  
-  private:
-
-    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
-
-    typedef std::vector<int> IntVector;
-    typedef std::vector<char> BoolVector;
-    typedef std::vector<Value> ValueVector;
-    typedef std::vector<Cost> CostVector;
-
-  private:
-
-    // Data related to the underlying digraph
-    const GR &_graph;
-    int _node_num;
-    int _arc_num;
-    int _res_arc_num;
-    int _root;
-
-    // Parameters of the problem
-    bool _have_lower;
-    Value _sum_supply;
-
-    // Data structures for storing the digraph
-    IntNodeMap _node_id;
-    IntArcMap _arc_idf;
-    IntArcMap _arc_idb;
-    IntVector _first_out;
-    BoolVector _forward;
-    IntVector _source;
-    IntVector _target;
-    IntVector _reverse;
-
-    // Node and arc data
-    ValueVector _lower;
-    ValueVector _upper;
-    CostVector _cost;
-    ValueVector _supply;
-
-    ValueVector _res_cap;
-    CostVector _pi;
-    ValueVector _excess;
-    IntVector _excess_nodes;
-    IntVector _deficit_nodes;
-
-    Value _delta;
-    int _factor;
-    IntVector _pred;
-
-  public:
-  
-    /// \brief Constant for infinite upper bounds (capacities).
-    ///
-    /// Constant for infinite upper bounds (capacities).
-    /// It is \c std::numeric_limits<Value>::infinity() if available,
-    /// \c std::numeric_limits<Value>::max() otherwise.
-    const Value INF;
-
-  private:
-
-    // Special implementation of the Dijkstra algorithm for finding
-    // shortest paths in the residual network of the digraph with
-    // respect to the reduced arc costs and modifying the node
-    // potentials according to the found distance labels.
-    class ResidualDijkstra
-    {
-    private:
-
-      int _node_num;
-      bool _geq;
-      const IntVector &_first_out;
-      const IntVector &_target;
-      const CostVector &_cost;
-      const ValueVector &_res_cap;
-      const ValueVector &_excess;
-      CostVector &_pi;
-      IntVector &_pred;
-      
-      IntVector _proc_nodes;
-      CostVector _dist;
-      
-    public:
-
-      ResidualDijkstra(CapacityScaling& cs) :
-        _node_num(cs._node_num), _geq(cs._sum_supply < 0),
-        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
-        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
-        _pred(cs._pred), _dist(cs._node_num)
-      {}
-
-      int run(int s, Value delta = 1) {
-        RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
-        Heap heap(heap_cross_ref);
-        heap.push(s, 0);
-        _pred[s] = -1;
-        _proc_nodes.clear();
-
-        // Process nodes
-        while (!heap.empty() && _excess[heap.top()] > -delta) {
-          int u = heap.top(), v;
-          Cost d = heap.prio() + _pi[u], dn;
-          _dist[u] = heap.prio();
-          _proc_nodes.push_back(u);
-          heap.pop();
-
-          // Traverse outgoing residual arcs
-          int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
-          for (int a = _first_out[u]; a != last_out; ++a) {
-            if (_res_cap[a] < delta) continue;
-            v = _target[a];
-            switch (heap.state(v)) {
-              case Heap::PRE_HEAP:
-                heap.push(v, d + _cost[a] - _pi[v]);
-                _pred[v] = a;
-                break;
-              case Heap::IN_HEAP:
-                dn = d + _cost[a] - _pi[v];
-                if (dn < heap[v]) {
-                  heap.decrease(v, dn);
-                  _pred[v] = a;
-                }
-                break;
-              case Heap::POST_HEAP:
-                break;
-            }
-          }
-        }
-        if (heap.empty()) return -1;
-
-        // Update potentials of processed nodes
-        int t = heap.top();
-        Cost dt = heap.prio();
-        for (int i = 0; i < int(_proc_nodes.size()); ++i) {
-          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
-        }
-
-        return t;
-      }
-
-    }; //class ResidualDijkstra
-
-  public:
-
-    /// \name Named Template Parameters
-    /// @{
-
-    template <typename T>
-    struct SetHeapTraits : public Traits {
-      typedef T Heap;
-    };
-
-    /// \brief \ref named-templ-param "Named parameter" for setting
-    /// \c Heap type.
-    ///
-    /// \ref named-templ-param "Named parameter" for setting \c Heap
-    /// type, which is used for internal Dijkstra computations.
-    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
-    /// its priority type must be \c Cost and its cross reference type
-    /// must be \ref RangeMap "RangeMap<int>".
-    template <typename T>
-    struct SetHeap
-      : public CapacityScaling<GR, V, C, SetHeapTraits<T> > {
-      typedef  CapacityScaling<GR, V, C, SetHeapTraits<T> > Create;
-    };
-
-    /// @}
-
-  public:
-
-    /// \brief Constructor.
-    ///
-    /// The constructor of the class.
-    ///
-    /// \param graph The digraph the algorithm runs on.
-    CapacityScaling(const GR& graph) :
-      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
-      INF(std::numeric_limits<Value>::has_infinity ?
-          std::numeric_limits<Value>::infinity() :
-          std::numeric_limits<Value>::max())
-    {
-      // Check the number types
-      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
-        "The flow type of CapacityScaling must be signed");
-      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
-        "The cost type of CapacityScaling must be signed");
-
-      // Reset data structures
-      reset();
-    }
-
-    /// \name Parameters
-    /// The parameters of the algorithm can be specified using these
-    /// functions.
-
-    /// @{
-
-    /// \brief Set the lower bounds on the arcs.
-    ///
-    /// This function sets the lower bounds on the arcs.
-    /// If it is not used before calling \ref run(), the lower bounds
-    /// will be set to zero on all arcs.
-    ///
-    /// \param map An arc map storing the lower bounds.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template <typename LowerMap>
-    CapacityScaling& lowerMap(const LowerMap& map) {
-      _have_lower = true;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _lower[_arc_idf[a]] = map[a];
-        _lower[_arc_idb[a]] = map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the upper bounds (capacities) on the arcs.
-    ///
-    /// This function sets the upper bounds (capacities) on the arcs.
-    /// If it is not used before calling \ref run(), the upper bounds
-    /// will be set to \ref INF on all arcs (i.e. the flow value will be
-    /// unbounded from above).
-    ///
-    /// \param map An arc map storing the upper bounds.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename UpperMap>
-    CapacityScaling& upperMap(const UpperMap& map) {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _upper[_arc_idf[a]] = map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the costs of the arcs.
-    ///
-    /// This function sets the costs of the arcs.
-    /// If it is not used before calling \ref run(), the costs
-    /// will be set to \c 1 on all arcs.
-    ///
-    /// \param map An arc map storing the costs.
-    /// Its \c Value type must be convertible to the \c Cost type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename CostMap>
-    CapacityScaling& costMap(const CostMap& map) {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _cost[_arc_idf[a]] =  map[a];
-        _cost[_arc_idb[a]] = -map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the supply values of the nodes.
-    ///
-    /// This function sets the supply values of the nodes.
-    /// If neither this function nor \ref stSupply() is used before
-    /// calling \ref run(), the supply of each node will be set to zero.
-    ///
-    /// \param map A node map storing the supply values.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename SupplyMap>
-    CapacityScaling& supplyMap(const SupplyMap& map) {
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        _supply[_node_id[n]] = map[n];
-      }
-      return *this;
-    }
-
-    /// \brief Set single source and target nodes and a supply value.
-    ///
-    /// This function sets a single source node and a single target node
-    /// and the required flow value.
-    /// If neither this function nor \ref supplyMap() is used before
-    /// calling \ref run(), the supply of each node will be set to zero.
-    ///
-    /// Using this function has the same effect as using \ref supplyMap()
-    /// with such a map in which \c k is assigned to \c s, \c -k is
-    /// assigned to \c t and all other nodes have zero supply value.
-    ///
-    /// \param s The source node.
-    /// \param t The target node.
-    /// \param k The required amount of flow from node \c s to node \c t
-    /// (i.e. the supply of \c s and the demand of \c t).
-    ///
-    /// \return <tt>(*this)</tt>
-    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
-      for (int i = 0; i != _node_num; ++i) {
-        _supply[i] = 0;
-      }
-      _supply[_node_id[s]] =  k;
-      _supply[_node_id[t]] = -k;
-      return *this;
-    }
-    
-    /// @}
-
-    /// \name Execution control
-    /// The algorithm can be executed using \ref run().
-
-    /// @{
-
-    /// \brief Run the algorithm.
-    ///
-    /// This function runs the algorithm.
-    /// The paramters can be specified using functions \ref lowerMap(),
-    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
-    /// For example,
-    /// \code
-    ///   CapacityScaling<ListDigraph> cs(graph);
-    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
-    ///     .supplyMap(sup).run();
-    /// \endcode
-    ///
-    /// This function can be called more than once. All the given parameters
-    /// are kept for the next call, unless \ref resetParams() or \ref reset()
-    /// is used, thus only the modified parameters have to be set again.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class (or the last \ref reset() call), then the \ref reset()
-    /// function must be called.
-    ///
-    /// \param factor The capacity scaling factor. It must be larger than
-    /// one to use scaling. If it is less or equal to one, then scaling
-    /// will be disabled.
-    ///
-    /// \return \c INFEASIBLE if no feasible flow exists,
-    /// \n \c OPTIMAL if the problem has optimal solution
-    /// (i.e. it is feasible and bounded), and the algorithm has found
-    /// optimal flow and node potentials (primal and dual solutions),
-    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
-    /// and infinite upper bound. It means that the objective function
-    /// is unbounded on that arc, however, note that it could actually be
-    /// bounded over the feasible flows, but this algroithm cannot handle
-    /// these cases.
-    ///
-    /// \see ProblemType
-    /// \see resetParams(), reset()
-    ProblemType run(int factor = 4) {
-      _factor = factor;
-      ProblemType pt = init();
-      if (pt != OPTIMAL) return pt;
-      return start();
-    }
-
-    /// \brief Reset all the parameters that have been given before.
-    ///
-    /// This function resets all the paramaters that have been given
-    /// before using functions \ref lowerMap(), \ref upperMap(),
-    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
-    ///
-    /// It is useful for multiple \ref run() calls. Basically, all the given
-    /// parameters are kept for the next \ref run() call, unless
-    /// \ref resetParams() or \ref reset() is used.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class or the last \ref reset() call, then the \ref reset()
-    /// function must be used, otherwise \ref resetParams() is sufficient.
-    ///
-    /// For example,
-    /// \code
-    ///   CapacityScaling<ListDigraph> cs(graph);
-    ///
-    ///   // First run
-    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
-    ///     .supplyMap(sup).run();
-    ///
-    ///   // Run again with modified cost map (resetParams() is not called,
-    ///   // so only the cost map have to be set again)
-    ///   cost[e] += 100;
-    ///   cs.costMap(cost).run();
-    ///
-    ///   // Run again from scratch using resetParams()
-    ///   // (the lower bounds will be set to zero on all arcs)
-    ///   cs.resetParams();
-    ///   cs.upperMap(capacity).costMap(cost)
-    ///     .supplyMap(sup).run();
-    /// \endcode
-    ///
-    /// \return <tt>(*this)</tt>
-    ///
-    /// \see reset(), run()
-    CapacityScaling& resetParams() {
-      for (int i = 0; i != _node_num; ++i) {
-        _supply[i] = 0;
-      }
-      for (int j = 0; j != _res_arc_num; ++j) {
-        _lower[j] = 0;
-        _upper[j] = INF;
-        _cost[j] = _forward[j] ? 1 : -1;
-      }
-      _have_lower = false;
-      return *this;
-    }
-
-    /// \brief Reset the internal data structures and all the parameters
-    /// that have been given before.
-    ///
-    /// This function resets the internal data structures and all the
-    /// paramaters that have been given before using functions \ref lowerMap(),
-    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
-    ///
-    /// It is useful for multiple \ref run() calls. Basically, all the given
-    /// parameters are kept for the next \ref run() call, unless
-    /// \ref resetParams() or \ref reset() is used.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class or the last \ref reset() call, then the \ref reset()
-    /// function must be used, otherwise \ref resetParams() is sufficient.
-    ///
-    /// See \ref resetParams() for examples.
-    ///
-    /// \return <tt>(*this)</tt>
-    ///
-    /// \see resetParams(), run()
-    CapacityScaling& reset() {
-      // Resize vectors
-      _node_num = countNodes(_graph);
-      _arc_num = countArcs(_graph);
-      _res_arc_num = 2 * (_arc_num + _node_num);
-      _root = _node_num;
-      ++_node_num;
-
-      _first_out.resize(_node_num + 1);
-      _forward.resize(_res_arc_num);
-      _source.resize(_res_arc_num);
-      _target.resize(_res_arc_num);
-      _reverse.resize(_res_arc_num);
-
-      _lower.resize(_res_arc_num);
-      _upper.resize(_res_arc_num);
-      _cost.resize(_res_arc_num);
-      _supply.resize(_node_num);
-      
-      _res_cap.resize(_res_arc_num);
-      _pi.resize(_node_num);
-      _excess.resize(_node_num);
-      _pred.resize(_node_num);
-
-      // Copy the graph
-      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
-      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-        _node_id[n] = i;
-      }
-      i = 0;
-      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-        _first_out[i] = j;
-        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
-          _arc_idf[a] = j;
-          _forward[j] = true;
-          _source[j] = i;
-          _target[j] = _node_id[_graph.runningNode(a)];
-        }
-        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
-          _arc_idb[a] = j;
-          _forward[j] = false;
-          _source[j] = i;
-          _target[j] = _node_id[_graph.runningNode(a)];
-        }
-        _forward[j] = false;
-        _source[j] = i;
-        _target[j] = _root;
-        _reverse[j] = k;
-        _forward[k] = true;
-        _source[k] = _root;
-        _target[k] = i;
-        _reverse[k] = j;
-        ++j; ++k;
-      }
-      _first_out[i] = j;
-      _first_out[_node_num] = k;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        int fi = _arc_idf[a];
-        int bi = _arc_idb[a];
-        _reverse[fi] = bi;
-        _reverse[bi] = fi;
-      }
-      
-      // Reset parameters
-      resetParams();
-      return *this;
-    }
-
-    /// @}
-
-    /// \name Query Functions
-    /// The results of the algorithm can be obtained using these
-    /// functions.\n
-    /// The \ref run() function must be called before using them.
-
-    /// @{
-
-    /// \brief Return the total cost of the found flow.
-    ///
-    /// This function returns the total cost of the found flow.
-    /// Its complexity is O(e).
-    ///
-    /// \note The return type of the function can be specified as a
-    /// template parameter. For example,
-    /// \code
-    ///   cs.totalCost<double>();
-    /// \endcode
-    /// It is useful if the total cost cannot be stored in the \c Cost
-    /// type of the algorithm, which is the default return type of the
-    /// function.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename Number>
-    Number totalCost() const {
-      Number c = 0;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        int i = _arc_idb[a];
-        c += static_cast<Number>(_res_cap[i]) *
-             (-static_cast<Number>(_cost[i]));
-      }
-      return c;
-    }
-
-#ifndef DOXYGEN
-    Cost totalCost() const {
-      return totalCost<Cost>();
-    }
-#endif
-
-    /// \brief Return the flow on the given arc.
-    ///
-    /// This function returns the flow on the given arc.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    Value flow(const Arc& a) const {
-      return _res_cap[_arc_idb[a]];
-    }
-
-    /// \brief Return the flow map (the primal solution).
-    ///
-    /// This function copies the flow value on each arc into the given
-    /// map. The \c Value type of the algorithm must be convertible to
-    /// the \c Value type of the map.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename FlowMap>
-    void flowMap(FlowMap &map) const {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        map.set(a, _res_cap[_arc_idb[a]]);
-      }
-    }
-
-    /// \brief Return the potential (dual value) of the given node.
-    ///
-    /// This function returns the potential (dual value) of the
-    /// given node.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    Cost potential(const Node& n) const {
-      return _pi[_node_id[n]];
-    }
-
-    /// \brief Return the potential map (the dual solution).
-    ///
-    /// This function copies the potential (dual value) of each node
-    /// into the given map.
-    /// The \c Cost type of the algorithm must be convertible to the
-    /// \c Value type of the map.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename PotentialMap>
-    void potentialMap(PotentialMap &map) const {
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        map.set(n, _pi[_node_id[n]]);
-      }
-    }
-
-    /// @}
-
-  private:
-
-    // Initialize the algorithm
-    ProblemType init() {
-      if (_node_num <= 1) return INFEASIBLE;
-
-      // Check the sum of supply values
-      _sum_supply = 0;
-      for (int i = 0; i != _root; ++i) {
-        _sum_supply += _supply[i];
-      }
-      if (_sum_supply > 0) return INFEASIBLE;
-      
-      // Initialize vectors
-      for (int i = 0; i != _root; ++i) {
-        _pi[i] = 0;
-        _excess[i] = _supply[i];
-      }
-
-      // Remove non-zero lower bounds
-      const Value MAX = std::numeric_limits<Value>::max();
-      int last_out;
-      if (_have_lower) {
-        for (int i = 0; i != _root; ++i) {
-          last_out = _first_out[i+1];
-          for (int j = _first_out[i]; j != last_out; ++j) {
-            if (_forward[j]) {
-              Value c = _lower[j];
-              if (c >= 0) {
-                _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
-              } else {
-                _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
-              }
-              _excess[i] -= c;
-              _excess[_target[j]] += c;
-            } else {
-              _res_cap[j] = 0;
-            }
-          }
-        }
-      } else {
-        for (int j = 0; j != _res_arc_num; ++j) {
-          _res_cap[j] = _forward[j] ? _upper[j] : 0;
-        }
-      }
-
-      // Handle negative costs
-      for (int i = 0; i != _root; ++i) {
-        last_out = _first_out[i+1] - 1;
-        for (int j = _first_out[i]; j != last_out; ++j) {
-          Value rc = _res_cap[j];
-          if (_cost[j] < 0 && rc > 0) {
-            if (rc >= MAX) return UNBOUNDED;
-            _excess[i] -= rc;
-            _excess[_target[j]] += rc;
-            _res_cap[j] = 0;
-            _res_cap[_reverse[j]] += rc;
-          }
-        }
-      }
-      
-      // Handle GEQ supply type
-      if (_sum_supply < 0) {
-        _pi[_root] = 0;
-        _excess[_root] = -_sum_supply;
-        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
-          int ra = _reverse[a];
-          _res_cap[a] = -_sum_supply + 1;
-          _res_cap[ra] = 0;
-          _cost[a] = 0;
-          _cost[ra] = 0;
-        }
-      } else {
-        _pi[_root] = 0;
-        _excess[_root] = 0;
-        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
-          int ra = _reverse[a];
-          _res_cap[a] = 1;
-          _res_cap[ra] = 0;
-          _cost[a] = 0;
-          _cost[ra] = 0;
-        }
-      }
-
-      // Initialize delta value
-      if (_factor > 1) {
-        // With scaling
-        Value max_sup = 0, max_dem = 0;
-        for (int i = 0; i != _node_num; ++i) {
-          Value ex = _excess[i];
-          if ( ex > max_sup) max_sup =  ex;
-          if (-ex > max_dem) max_dem = -ex;
-        }
-        Value max_cap = 0;
-        for (int j = 0; j != _res_arc_num; ++j) {
-          if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
-        }
-        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
-        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
-      } else {
-        // Without scaling
-        _delta = 1;
-      }
-
-      return OPTIMAL;
-    }
-
-    ProblemType start() {
-      // Execute the algorithm
-      ProblemType pt;
-      if (_delta > 1)
-        pt = startWithScaling();
-      else
-        pt = startWithoutScaling();
-
-      // Handle non-zero lower bounds
-      if (_have_lower) {
-        int limit = _first_out[_root];
-        for (int j = 0; j != limit; ++j) {
-          if (!_forward[j]) _res_cap[j] += _lower[j];
-        }
-      }
-
-      // Shift potentials if necessary
-      Cost pr = _pi[_root];
-      if (_sum_supply < 0 || pr > 0) {
-        for (int i = 0; i != _node_num; ++i) {
-          _pi[i] -= pr;
-        }        
-      }
-      
-      return pt;
-    }
-
-    // Execute the capacity scaling algorithm
-    ProblemType startWithScaling() {
-      // Perform capacity scaling phases
-      int s, t;
-      ResidualDijkstra _dijkstra(*this);
-      while (true) {
-        // Saturate all arcs not satisfying the optimality condition
-        int last_out;
-        for (int u = 0; u != _node_num; ++u) {
-          last_out = _sum_supply < 0 ?
-            _first_out[u+1] : _first_out[u+1] - 1;
-          for (int a = _first_out[u]; a != last_out; ++a) {
-            int v = _target[a];
-            Cost c = _cost[a] + _pi[u] - _pi[v];
-            Value rc = _res_cap[a];
-            if (c < 0 && rc >= _delta) {
-              _excess[u] -= rc;
-              _excess[v] += rc;
-              _res_cap[a] = 0;
-              _res_cap[_reverse[a]] += rc;
-            }
-          }
-        }
-
-        // Find excess nodes and deficit nodes
-        _excess_nodes.clear();
-        _deficit_nodes.clear();
-        for (int u = 0; u != _node_num; ++u) {
-          Value ex = _excess[u];
-          if (ex >=  _delta) _excess_nodes.push_back(u);
-          if (ex <= -_delta) _deficit_nodes.push_back(u);
-        }
-        int next_node = 0, next_def_node = 0;
-
-        // Find augmenting shortest paths
-        while (next_node < int(_excess_nodes.size())) {
-          // Check deficit nodes
-          if (_delta > 1) {
-            bool delta_deficit = false;
-            for ( ; next_def_node < int(_deficit_nodes.size());
-                    ++next_def_node ) {
-              if (_excess[_deficit_nodes[next_def_node]] <= -_delta) {
-                delta_deficit = true;
-                break;
-              }
-            }
-            if (!delta_deficit) break;
-          }
-
-          // Run Dijkstra in the residual network
-          s = _excess_nodes[next_node];
-          if ((t = _dijkstra.run(s, _delta)) == -1) {
-            if (_delta > 1) {
-              ++next_node;
-              continue;
-            }
-            return INFEASIBLE;
-          }
-
-          // Augment along a shortest path from s to t
-          Value d = std::min(_excess[s], -_excess[t]);
-          int u = t;
-          int a;
-          if (d > _delta) {
-            while ((a = _pred[u]) != -1) {
-              if (_res_cap[a] < d) d = _res_cap[a];
-              u = _source[a];
-            }
-          }
-          u = t;
-          while ((a = _pred[u]) != -1) {
-            _res_cap[a] -= d;
-            _res_cap[_reverse[a]] += d;
-            u = _source[a];
-          }
-          _excess[s] -= d;
-          _excess[t] += d;
-
-          if (_excess[s] < _delta) ++next_node;
-        }
-
-        if (_delta == 1) break;
-        _delta = _delta <= _factor ? 1 : _delta / _factor;
-      }
-
-      return OPTIMAL;
-    }
-
-    // Execute the successive shortest path algorithm
-    ProblemType startWithoutScaling() {
-      // Find excess nodes
-      _excess_nodes.clear();
-      for (int i = 0; i != _node_num; ++i) {
-        if (_excess[i] > 0) _excess_nodes.push_back(i);
-      }
-      if (_excess_nodes.size() == 0) return OPTIMAL;
-      int next_node = 0;
-
-      // Find shortest paths
-      int s, t;
-      ResidualDijkstra _dijkstra(*this);
-      while ( _excess[_excess_nodes[next_node]] > 0 ||
-              ++next_node < int(_excess_nodes.size()) )
-      {
-        // Run Dijkstra in the residual network
-        s = _excess_nodes[next_node];
-        if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE;
-
-        // Augment along a shortest path from s to t
-        Value d = std::min(_excess[s], -_excess[t]);
-        int u = t;
-        int a;
-        if (d > 1) {
-          while ((a = _pred[u]) != -1) {
-            if (_res_cap[a] < d) d = _res_cap[a];
-            u = _source[a];
-          }
-        }
-        u = t;
-        while ((a = _pred[u]) != -1) {
-          _res_cap[a] -= d;
-          _res_cap[_reverse[a]] += d;
-          u = _source[a];
-        }
-        _excess[s] -= d;
-        _excess[t] += d;
-      }
-
-      return OPTIMAL;
-    }
-
-  }; //class CapacityScaling
-
-  ///@}
-
-} //namespace lemon
-
-#endif //LEMON_CAPACITY_SCALING_H
Index: lemon/circulation.h
===================================================================
--- lemon/circulation.h	(revision 891)
+++ lemon/circulation.h	(revision 833)
@@ -174,9 +174,4 @@
      \tparam SM The type of the supply map. The default map type is
      \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
-     \tparam TR The traits class that defines various types used by the
-     algorithm. By default, it is \ref CirculationDefaultTraits
-     "CirculationDefaultTraits<GR, LM, UM, SM>".
-     In most cases, this parameter should not be set directly,
-     consider to use the named template parameters instead.
   */
 #ifdef DOXYGEN
Index: lemon/concepts/heap.h
===================================================================
--- lemon/concepts/heap.h	(revision 883)
+++ lemon/concepts/heap.h	(revision 757)
@@ -89,9 +89,5 @@
       /// handle the cross references. The assigned value must be
       /// \c PRE_HEAP (<tt>-1</tt>) for each item.
-#ifdef DOXYGEN
       explicit Heap(ItemIntMap &map) {}
-#else
-      explicit Heap(ItemIntMap&) {}
-#endif      
 
       /// \brief Constructor.
@@ -103,9 +99,5 @@
       /// \c PRE_HEAP (<tt>-1</tt>) for each item.
       /// \param comp The function object used for comparing the priorities.
-#ifdef DOXYGEN
       explicit Heap(ItemIntMap &map, const CMP &comp) {}
-#else
-      explicit Heap(ItemIntMap&, const CMP&) {}
-#endif      
 
       /// \brief The number of items stored in the heap.
@@ -135,9 +127,5 @@
       /// \param p The priority of the item.
       /// \pre \e i must not be stored in the heap.
-#ifdef DOXYGEN
       void push(const Item &i, const Prio &p) {}
-#else
-      void push(const Item&, const Prio&) {}
-#endif      
 
       /// \brief Return the item having minimum priority.
@@ -145,5 +133,5 @@
       /// This function returns the item having minimum priority.
       /// \pre The heap must be non-empty.
-      Item top() const { return Item(); }
+      Item top() const {}
 
       /// \brief The minimum priority.
@@ -151,5 +139,5 @@
       /// This function returns the minimum priority.
       /// \pre The heap must be non-empty.
-      Prio prio() const { return Prio(); }
+      Prio prio() const {}
 
       /// \brief Remove the item having minimum priority.
@@ -165,9 +153,5 @@
       /// \param i The item to delete.
       /// \pre \e i must be in the heap.
-#ifdef DOXYGEN
       void erase(const Item &i) {}
-#else
-      void erase(const Item&) {}
-#endif      
 
       /// \brief The priority of the given item.
@@ -176,9 +160,5 @@
       /// \param i The item.
       /// \pre \e i must be in the heap.
-#ifdef DOXYGEN
       Prio operator[](const Item &i) const {}
-#else
-      Prio operator[](const Item&) const { return Prio(); }
-#endif      
 
       /// \brief Set the priority of an item or insert it, if it is
@@ -191,9 +171,5 @@
       /// \param i The item.
       /// \param p The priority.
-#ifdef DOXYGEN
       void set(const Item &i, const Prio &p) {}
-#else
-      void set(const Item&, const Prio&) {}
-#endif      
 
       /// \brief Decrease the priority of an item to the given value.
@@ -203,9 +179,5 @@
       /// \param p The priority.
       /// \pre \e i must be stored in the heap with priority at least \e p.
-#ifdef DOXYGEN
       void decrease(const Item &i, const Prio &p) {}
-#else
-      void decrease(const Item&, const Prio&) {}
-#endif      
 
       /// \brief Increase the priority of an item to the given value.
@@ -215,9 +187,5 @@
       /// \param p The priority.
       /// \pre \e i must be stored in the heap with priority at most \e p.
-#ifdef DOXYGEN
       void increase(const Item &i, const Prio &p) {}
-#else
-      void increase(const Item&, const Prio&) {}
-#endif      
 
       /// \brief Return the state of an item.
@@ -229,9 +197,5 @@
       /// to the heap again.
       /// \param i The item.
-#ifdef DOXYGEN
       State state(const Item &i) const {}
-#else
-      State state(const Item&) const { return PRE_HEAP; }
-#endif      
 
       /// \brief Set the state of an item in the heap.
@@ -242,9 +206,5 @@
       /// \param i The item.
       /// \param st The state. It should not be \c IN_HEAP.
-#ifdef DOXYGEN
       void state(const Item& i, State st) {}
-#else
-      void state(const Item&, State) {}
-#endif      
 
 
Index: mon/cost_scaling.h
===================================================================
--- lemon/cost_scaling.h	(revision 899)
+++ 	(revision )
@@ -1,1196 +1,0 @@
-/* -*- C++ -*-
- *
- * This file is a part of LEMON, a generic C++ optimization library
- *
- * Copyright (C) 2003-2008
- * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
- * (Egervary Research Group on Combinatorial Optimization, EGRES).
- *
- * Permission to use, modify and distribute this software is granted
- * provided that this copyright notice appears in all copies. For
- * precise terms see the accompanying LICENSE file.
- *
- * This software is provided "AS IS" with no warranty of any kind,
- * express or implied, and with no claim as to its suitability for any
- * purpose.
- *
- */
-
-#ifndef LEMON_COST_SCALING_H
-#define LEMON_COST_SCALING_H
-
-/// \ingroup min_cost_flow_algs
-/// \file
-/// \brief Cost scaling algorithm for finding a minimum cost flow.
-
-#include <vector>
-#include <deque>
-#include <limits>
-
-#include <lemon/core.h>
-#include <lemon/maps.h>
-#include <lemon/math.h>
-#include <lemon/static_graph.h>
-#include <lemon/circulation.h>
-#include <lemon/bellman_ford.h>
-
-namespace lemon {
-
-  /// \brief Default traits class of CostScaling algorithm.
-  ///
-  /// Default traits class of CostScaling algorithm.
-  /// \tparam GR Digraph type.
-  /// \tparam V The number type used for flow amounts, capacity bounds
-  /// and supply values. By default it is \c int.
-  /// \tparam C The number type used for costs and potentials.
-  /// By default it is the same as \c V.
-#ifdef DOXYGEN
-  template <typename GR, typename V = int, typename C = V>
-#else
-  template < typename GR, typename V = int, typename C = V,
-             bool integer = std::numeric_limits<C>::is_integer >
-#endif
-  struct CostScalingDefaultTraits
-  {
-    /// The type of the digraph
-    typedef GR Digraph;
-    /// The type of the flow amounts, capacity bounds and supply values
-    typedef V Value;
-    /// The type of the arc costs
-    typedef C Cost;
-
-    /// \brief The large cost type used for internal computations
-    ///
-    /// The large cost type used for internal computations.
-    /// It is \c long \c long if the \c Cost type is integer,
-    /// otherwise it is \c double.
-    /// \c Cost must be convertible to \c LargeCost.
-    typedef double LargeCost;
-  };
-
-  // Default traits class for integer cost types
-  template <typename GR, typename V, typename C>
-  struct CostScalingDefaultTraits<GR, V, C, true>
-  {
-    typedef GR Digraph;
-    typedef V Value;
-    typedef C Cost;
-#ifdef LEMON_HAVE_LONG_LONG
-    typedef long long LargeCost;
-#else
-    typedef long LargeCost;
-#endif
-  };
-
-
-  /// \addtogroup min_cost_flow_algs
-  /// @{
-
-  /// \brief Implementation of the Cost Scaling algorithm for
-  /// finding a \ref min_cost_flow "minimum cost flow".
-  ///
-  /// \ref CostScaling implements a cost scaling algorithm that performs
-  /// push/augment and relabel operations for finding a \ref min_cost_flow
-  /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
-  /// \ref goldberg97efficient, \ref bunnagel98efficient. 
-  /// It is a highly efficient primal-dual solution method, which
-  /// can be viewed as the generalization of the \ref Preflow
-  /// "preflow push-relabel" algorithm for the maximum flow problem.
-  ///
-  /// Most of the parameters of the problem (except for the digraph)
-  /// can be given using separate functions, and the algorithm can be
-  /// executed using the \ref run() function. If some parameters are not
-  /// specified, then default values will be used.
-  ///
-  /// \tparam GR The digraph type the algorithm runs on.
-  /// \tparam V The number type used for flow amounts, capacity bounds
-  /// and supply values in the algorithm. By default, it is \c int.
-  /// \tparam C The number type used for costs and potentials in the
-  /// algorithm. By default, it is the same as \c V.
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref CostScalingDefaultTraits
-  /// "CostScalingDefaultTraits<GR, V, C>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
-  ///
-  /// \warning Both number types must be signed and all input data must
-  /// be integer.
-  /// \warning This algorithm does not support negative costs for such
-  /// arcs that have infinite upper bound.
-  ///
-  /// \note %CostScaling provides three different internal methods,
-  /// from which the most efficient one is used by default.
-  /// For more information, see \ref Method.
-#ifdef DOXYGEN
-  template <typename GR, typename V, typename C, typename TR>
-#else
-  template < typename GR, typename V = int, typename C = V,
-             typename TR = CostScalingDefaultTraits<GR, V, C> >
-#endif
-  class CostScaling
-  {
-  public:
-
-    /// The type of the digraph
-    typedef typename TR::Digraph Digraph;
-    /// The type of the flow amounts, capacity bounds and supply values
-    typedef typename TR::Value Value;
-    /// The type of the arc costs
-    typedef typename TR::Cost Cost;
-
-    /// \brief The large cost type
-    ///
-    /// The large cost type used for internal computations.
-    /// By default, it is \c long \c long if the \c Cost type is integer,
-    /// otherwise it is \c double.
-    typedef typename TR::LargeCost LargeCost;
-
-    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
-    typedef TR Traits;
-
-  public:
-
-    /// \brief Problem type constants for the \c run() function.
-    ///
-    /// Enum type containing the problem type constants that can be
-    /// returned by the \ref run() function of the algorithm.
-    enum ProblemType {
-      /// The problem has no feasible solution (flow).
-      INFEASIBLE,
-      /// The problem has optimal solution (i.e. it is feasible and
-      /// bounded), and the algorithm has found optimal flow and node
-      /// potentials (primal and dual solutions).
-      OPTIMAL,
-      /// The digraph contains an arc of negative cost and infinite
-      /// upper bound. It means that the objective function is unbounded
-      /// on that arc, however, note that it could actually be bounded
-      /// over the feasible flows, but this algroithm cannot handle
-      /// these cases.
-      UNBOUNDED
-    };
-
-    /// \brief Constants for selecting the internal method.
-    ///
-    /// Enum type containing constants for selecting the internal method
-    /// for the \ref run() function.
-    ///
-    /// \ref CostScaling provides three internal methods that differ mainly
-    /// in their base operations, which are used in conjunction with the
-    /// relabel operation.
-    /// By default, the so called \ref PARTIAL_AUGMENT
-    /// "Partial Augment-Relabel" method is used, which proved to be
-    /// the most efficient and the most robust on various test inputs.
-    /// However, the other methods can be selected using the \ref run()
-    /// function with the proper parameter.
-    enum Method {
-      /// Local push operations are used, i.e. flow is moved only on one
-      /// admissible arc at once.
-      PUSH,
-      /// Augment operations are used, i.e. flow is moved on admissible
-      /// paths from a node with excess to a node with deficit.
-      AUGMENT,
-      /// Partial augment operations are used, i.e. flow is moved on 
-      /// admissible paths started from a node with excess, but the
-      /// lengths of these paths are limited. This method can be viewed
-      /// as a combined version of the previous two operations.
-      PARTIAL_AUGMENT
-    };
-
-  private:
-
-    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
-
-    typedef std::vector<int> IntVector;
-    typedef std::vector<char> BoolVector;
-    typedef std::vector<Value> ValueVector;
-    typedef std::vector<Cost> CostVector;
-    typedef std::vector<LargeCost> LargeCostVector;
-
-  private:
-  
-    template <typename KT, typename VT>
-    class StaticVectorMap {
-    public:
-      typedef KT Key;
-      typedef VT Value;
-      
-      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
-      
-      const Value& operator[](const Key& key) const {
-        return _v[StaticDigraph::id(key)];
-      }
-
-      Value& operator[](const Key& key) {
-        return _v[StaticDigraph::id(key)];
-      }
-      
-      void set(const Key& key, const Value& val) {
-        _v[StaticDigraph::id(key)] = val;
-      }
-
-    private:
-      std::vector<Value>& _v;
-    };
-
-    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
-    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
-
-  private:
-
-    // Data related to the underlying digraph
-    const GR &_graph;
-    int _node_num;
-    int _arc_num;
-    int _res_node_num;
-    int _res_arc_num;
-    int _root;
-
-    // Parameters of the problem
-    bool _have_lower;
-    Value _sum_supply;
-
-    // Data structures for storing the digraph
-    IntNodeMap _node_id;
-    IntArcMap _arc_idf;
-    IntArcMap _arc_idb;
-    IntVector _first_out;
-    BoolVector _forward;
-    IntVector _source;
-    IntVector _target;
-    IntVector _reverse;
-
-    // Node and arc data
-    ValueVector _lower;
-    ValueVector _upper;
-    CostVector _scost;
-    ValueVector _supply;
-
-    ValueVector _res_cap;
-    LargeCostVector _cost;
-    LargeCostVector _pi;
-    ValueVector _excess;
-    IntVector _next_out;
-    std::deque<int> _active_nodes;
-
-    // Data for scaling
-    LargeCost _epsilon;
-    int _alpha;
-
-    // Data for a StaticDigraph structure
-    typedef std::pair<int, int> IntPair;
-    StaticDigraph _sgr;
-    std::vector<IntPair> _arc_vec;
-    std::vector<LargeCost> _cost_vec;
-    LargeCostArcMap _cost_map;
-    LargeCostNodeMap _pi_map;
-  
-  public:
-  
-    /// \brief Constant for infinite upper bounds (capacities).
-    ///
-    /// Constant for infinite upper bounds (capacities).
-    /// It is \c std::numeric_limits<Value>::infinity() if available,
-    /// \c std::numeric_limits<Value>::max() otherwise.
-    const Value INF;
-
-  public:
-
-    /// \name Named Template Parameters
-    /// @{
-
-    template <typename T>
-    struct SetLargeCostTraits : public Traits {
-      typedef T LargeCost;
-    };
-
-    /// \brief \ref named-templ-param "Named parameter" for setting
-    /// \c LargeCost type.
-    ///
-    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
-    /// type, which is used for internal computations in the algorithm.
-    /// \c Cost must be convertible to \c LargeCost.
-    template <typename T>
-    struct SetLargeCost
-      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
-      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
-    };
-
-    /// @}
-
-  public:
-
-    /// \brief Constructor.
-    ///
-    /// The constructor of the class.
-    ///
-    /// \param graph The digraph the algorithm runs on.
-    CostScaling(const GR& graph) :
-      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
-      _cost_map(_cost_vec), _pi_map(_pi),
-      INF(std::numeric_limits<Value>::has_infinity ?
-          std::numeric_limits<Value>::infinity() :
-          std::numeric_limits<Value>::max())
-    {
-      // Check the number types
-      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
-        "The flow type of CostScaling must be signed");
-      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
-        "The cost type of CostScaling must be signed");
-      
-      // Reset data structures
-      reset();
-    }
-
-    /// \name Parameters
-    /// The parameters of the algorithm can be specified using these
-    /// functions.
-
-    /// @{
-
-    /// \brief Set the lower bounds on the arcs.
-    ///
-    /// This function sets the lower bounds on the arcs.
-    /// If it is not used before calling \ref run(), the lower bounds
-    /// will be set to zero on all arcs.
-    ///
-    /// \param map An arc map storing the lower bounds.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template <typename LowerMap>
-    CostScaling& lowerMap(const LowerMap& map) {
-      _have_lower = true;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _lower[_arc_idf[a]] = map[a];
-        _lower[_arc_idb[a]] = map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the upper bounds (capacities) on the arcs.
-    ///
-    /// This function sets the upper bounds (capacities) on the arcs.
-    /// If it is not used before calling \ref run(), the upper bounds
-    /// will be set to \ref INF on all arcs (i.e. the flow value will be
-    /// unbounded from above).
-    ///
-    /// \param map An arc map storing the upper bounds.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename UpperMap>
-    CostScaling& upperMap(const UpperMap& map) {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _upper[_arc_idf[a]] = map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the costs of the arcs.
-    ///
-    /// This function sets the costs of the arcs.
-    /// If it is not used before calling \ref run(), the costs
-    /// will be set to \c 1 on all arcs.
-    ///
-    /// \param map An arc map storing the costs.
-    /// Its \c Value type must be convertible to the \c Cost type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename CostMap>
-    CostScaling& costMap(const CostMap& map) {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _scost[_arc_idf[a]] =  map[a];
-        _scost[_arc_idb[a]] = -map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the supply values of the nodes.
-    ///
-    /// This function sets the supply values of the nodes.
-    /// If neither this function nor \ref stSupply() is used before
-    /// calling \ref run(), the supply of each node will be set to zero.
-    ///
-    /// \param map A node map storing the supply values.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename SupplyMap>
-    CostScaling& supplyMap(const SupplyMap& map) {
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        _supply[_node_id[n]] = map[n];
-      }
-      return *this;
-    }
-
-    /// \brief Set single source and target nodes and a supply value.
-    ///
-    /// This function sets a single source node and a single target node
-    /// and the required flow value.
-    /// If neither this function nor \ref supplyMap() is used before
-    /// calling \ref run(), the supply of each node will be set to zero.
-    ///
-    /// Using this function has the same effect as using \ref supplyMap()
-    /// with such a map in which \c k is assigned to \c s, \c -k is
-    /// assigned to \c t and all other nodes have zero supply value.
-    ///
-    /// \param s The source node.
-    /// \param t The target node.
-    /// \param k The required amount of flow from node \c s to node \c t
-    /// (i.e. the supply of \c s and the demand of \c t).
-    ///
-    /// \return <tt>(*this)</tt>
-    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
-      for (int i = 0; i != _res_node_num; ++i) {
-        _supply[i] = 0;
-      }
-      _supply[_node_id[s]] =  k;
-      _supply[_node_id[t]] = -k;
-      return *this;
-    }
-    
-    /// @}
-
-    /// \name Execution control
-    /// The algorithm can be executed using \ref run().
-
-    /// @{
-
-    /// \brief Run the algorithm.
-    ///
-    /// This function runs the algorithm.
-    /// The paramters can be specified using functions \ref lowerMap(),
-    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
-    /// For example,
-    /// \code
-    ///   CostScaling<ListDigraph> cs(graph);
-    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
-    ///     .supplyMap(sup).run();
-    /// \endcode
-    ///
-    /// This function can be called more than once. All the given parameters
-    /// are kept for the next call, unless \ref resetParams() or \ref reset()
-    /// is used, thus only the modified parameters have to be set again.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class (or the last \ref reset() call), then the \ref reset()
-    /// function must be called.
-    ///
-    /// \param method The internal method that will be used in the
-    /// algorithm. For more information, see \ref Method.
-    /// \param factor The cost scaling factor. It must be larger than one.
-    ///
-    /// \return \c INFEASIBLE if no feasible flow exists,
-    /// \n \c OPTIMAL if the problem has optimal solution
-    /// (i.e. it is feasible and bounded), and the algorithm has found
-    /// optimal flow and node potentials (primal and dual solutions),
-    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
-    /// and infinite upper bound. It means that the objective function
-    /// is unbounded on that arc, however, note that it could actually be
-    /// bounded over the feasible flows, but this algroithm cannot handle
-    /// these cases.
-    ///
-    /// \see ProblemType, Method
-    /// \see resetParams(), reset()
-    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
-      _alpha = factor;
-      ProblemType pt = init();
-      if (pt != OPTIMAL) return pt;
-      start(method);
-      return OPTIMAL;
-    }
-
-    /// \brief Reset all the parameters that have been given before.
-    ///
-    /// This function resets all the paramaters that have been given
-    /// before using functions \ref lowerMap(), \ref upperMap(),
-    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
-    ///
-    /// It is useful for multiple \ref run() calls. Basically, all the given
-    /// parameters are kept for the next \ref run() call, unless
-    /// \ref resetParams() or \ref reset() is used.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class or the last \ref reset() call, then the \ref reset()
-    /// function must be used, otherwise \ref resetParams() is sufficient.
-    ///
-    /// For example,
-    /// \code
-    ///   CostScaling<ListDigraph> cs(graph);
-    ///
-    ///   // First run
-    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
-    ///     .supplyMap(sup).run();
-    ///
-    ///   // Run again with modified cost map (resetParams() is not called,
-    ///   // so only the cost map have to be set again)
-    ///   cost[e] += 100;
-    ///   cs.costMap(cost).run();
-    ///
-    ///   // Run again from scratch using resetParams()
-    ///   // (the lower bounds will be set to zero on all arcs)
-    ///   cs.resetParams();
-    ///   cs.upperMap(capacity).costMap(cost)
-    ///     .supplyMap(sup).run();
-    /// \endcode
-    ///
-    /// \return <tt>(*this)</tt>
-    ///
-    /// \see reset(), run()
-    CostScaling& resetParams() {
-      for (int i = 0; i != _res_node_num; ++i) {
-        _supply[i] = 0;
-      }
-      int limit = _first_out[_root];
-      for (int j = 0; j != limit; ++j) {
-        _lower[j] = 0;
-        _upper[j] = INF;
-        _scost[j] = _forward[j] ? 1 : -1;
-      }
-      for (int j = limit; j != _res_arc_num; ++j) {
-        _lower[j] = 0;
-        _upper[j] = INF;
-        _scost[j] = 0;
-        _scost[_reverse[j]] = 0;
-      }      
-      _have_lower = false;
-      return *this;
-    }
-
-    /// \brief Reset all the parameters that have been given before.
-    ///
-    /// This function resets all the paramaters that have been given
-    /// before using functions \ref lowerMap(), \ref upperMap(),
-    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
-    ///
-    /// It is useful for multiple run() calls. If this function is not
-    /// used, all the parameters given before are kept for the next
-    /// \ref run() call.
-    /// However, the underlying digraph must not be modified after this
-    /// class have been constructed, since it copies and extends the graph.
-    /// \return <tt>(*this)</tt>
-    CostScaling& reset() {
-      // Resize vectors
-      _node_num = countNodes(_graph);
-      _arc_num = countArcs(_graph);
-      _res_node_num = _node_num + 1;
-      _res_arc_num = 2 * (_arc_num + _node_num);
-      _root = _node_num;
-
-      _first_out.resize(_res_node_num + 1);
-      _forward.resize(_res_arc_num);
-      _source.resize(_res_arc_num);
-      _target.resize(_res_arc_num);
-      _reverse.resize(_res_arc_num);
-
-      _lower.resize(_res_arc_num);
-      _upper.resize(_res_arc_num);
-      _scost.resize(_res_arc_num);
-      _supply.resize(_res_node_num);
-      
-      _res_cap.resize(_res_arc_num);
-      _cost.resize(_res_arc_num);
-      _pi.resize(_res_node_num);
-      _excess.resize(_res_node_num);
-      _next_out.resize(_res_node_num);
-
-      _arc_vec.reserve(_res_arc_num);
-      _cost_vec.reserve(_res_arc_num);
-
-      // Copy the graph
-      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
-      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-        _node_id[n] = i;
-      }
-      i = 0;
-      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-        _first_out[i] = j;
-        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
-          _arc_idf[a] = j;
-          _forward[j] = true;
-          _source[j] = i;
-          _target[j] = _node_id[_graph.runningNode(a)];
-        }
-        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
-          _arc_idb[a] = j;
-          _forward[j] = false;
-          _source[j] = i;
-          _target[j] = _node_id[_graph.runningNode(a)];
-        }
-        _forward[j] = false;
-        _source[j] = i;
-        _target[j] = _root;
-        _reverse[j] = k;
-        _forward[k] = true;
-        _source[k] = _root;
-        _target[k] = i;
-        _reverse[k] = j;
-        ++j; ++k;
-      }
-      _first_out[i] = j;
-      _first_out[_res_node_num] = k;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        int fi = _arc_idf[a];
-        int bi = _arc_idb[a];
-        _reverse[fi] = bi;
-        _reverse[bi] = fi;
-      }
-      
-      // Reset parameters
-      resetParams();
-      return *this;
-    }
-
-    /// @}
-
-    /// \name Query Functions
-    /// The results of the algorithm can be obtained using these
-    /// functions.\n
-    /// The \ref run() function must be called before using them.
-
-    /// @{
-
-    /// \brief Return the total cost of the found flow.
-    ///
-    /// This function returns the total cost of the found flow.
-    /// Its complexity is O(e).
-    ///
-    /// \note The return type of the function can be specified as a
-    /// template parameter. For example,
-    /// \code
-    ///   cs.totalCost<double>();
-    /// \endcode
-    /// It is useful if the total cost cannot be stored in the \c Cost
-    /// type of the algorithm, which is the default return type of the
-    /// function.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename Number>
-    Number totalCost() const {
-      Number c = 0;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        int i = _arc_idb[a];
-        c += static_cast<Number>(_res_cap[i]) *
-             (-static_cast<Number>(_scost[i]));
-      }
-      return c;
-    }
-
-#ifndef DOXYGEN
-    Cost totalCost() const {
-      return totalCost<Cost>();
-    }
-#endif
-
-    /// \brief Return the flow on the given arc.
-    ///
-    /// This function returns the flow on the given arc.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    Value flow(const Arc& a) const {
-      return _res_cap[_arc_idb[a]];
-    }
-
-    /// \brief Return the flow map (the primal solution).
-    ///
-    /// This function copies the flow value on each arc into the given
-    /// map. The \c Value type of the algorithm must be convertible to
-    /// the \c Value type of the map.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename FlowMap>
-    void flowMap(FlowMap &map) const {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        map.set(a, _res_cap[_arc_idb[a]]);
-      }
-    }
-
-    /// \brief Return the potential (dual value) of the given node.
-    ///
-    /// This function returns the potential (dual value) of the
-    /// given node.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    Cost potential(const Node& n) const {
-      return static_cast<Cost>(_pi[_node_id[n]]);
-    }
-
-    /// \brief Return the potential map (the dual solution).
-    ///
-    /// This function copies the potential (dual value) of each node
-    /// into the given map.
-    /// The \c Cost type of the algorithm must be convertible to the
-    /// \c Value type of the map.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename PotentialMap>
-    void potentialMap(PotentialMap &map) const {
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
-      }
-    }
-
-    /// @}
-
-  private:
-
-    // Initialize the algorithm
-    ProblemType init() {
-      if (_res_node_num <= 1) return INFEASIBLE;
-
-      // Check the sum of supply values
-      _sum_supply = 0;
-      for (int i = 0; i != _root; ++i) {
-        _sum_supply += _supply[i];
-      }
-      if (_sum_supply > 0) return INFEASIBLE;
-      
-
-      // Initialize vectors
-      for (int i = 0; i != _res_node_num; ++i) {
-        _pi[i] = 0;
-        _excess[i] = _supply[i];
-      }
-      
-      // Remove infinite upper bounds and check negative arcs
-      const Value MAX = std::numeric_limits<Value>::max();
-      int last_out;
-      if (_have_lower) {
-        for (int i = 0; i != _root; ++i) {
-          last_out = _first_out[i+1];
-          for (int j = _first_out[i]; j != last_out; ++j) {
-            if (_forward[j]) {
-              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
-              if (c >= MAX) return UNBOUNDED;
-              _excess[i] -= c;
-              _excess[_target[j]] += c;
-            }
-          }
-        }
-      } else {
-        for (int i = 0; i != _root; ++i) {
-          last_out = _first_out[i+1];
-          for (int j = _first_out[i]; j != last_out; ++j) {
-            if (_forward[j] && _scost[j] < 0) {
-              Value c = _upper[j];
-              if (c >= MAX) return UNBOUNDED;
-              _excess[i] -= c;
-              _excess[_target[j]] += c;
-            }
-          }
-        }
-      }
-      Value ex, max_cap = 0;
-      for (int i = 0; i != _res_node_num; ++i) {
-        ex = _excess[i];
-        _excess[i] = 0;
-        if (ex < 0) max_cap -= ex;
-      }
-      for (int j = 0; j != _res_arc_num; ++j) {
-        if (_upper[j] >= MAX) _upper[j] = max_cap;
-      }
-
-      // Initialize the large cost vector and the epsilon parameter
-      _epsilon = 0;
-      LargeCost lc;
-      for (int i = 0; i != _root; ++i) {
-        last_out = _first_out[i+1];
-        for (int j = _first_out[i]; j != last_out; ++j) {
-          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
-          _cost[j] = lc;
-          if (lc > _epsilon) _epsilon = lc;
-        }
-      }
-      _epsilon /= _alpha;
-
-      // Initialize maps for Circulation and remove non-zero lower bounds
-      ConstMap<Arc, Value> low(0);
-      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
-      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
-      ValueArcMap cap(_graph), flow(_graph);
-      ValueNodeMap sup(_graph);
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        sup[n] = _supply[_node_id[n]];
-      }
-      if (_have_lower) {
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          int j = _arc_idf[a];
-          Value c = _lower[j];
-          cap[a] = _upper[j] - c;
-          sup[_graph.source(a)] -= c;
-          sup[_graph.target(a)] += c;
-        }
-      } else {
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          cap[a] = _upper[_arc_idf[a]];
-        }
-      }
-
-      // Find a feasible flow using Circulation
-      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
-        circ(_graph, low, cap, sup);
-      if (!circ.flowMap(flow).run()) return INFEASIBLE;
-
-      // Set residual capacities and handle GEQ supply type
-      if (_sum_supply < 0) {
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          Value fa = flow[a];
-          _res_cap[_arc_idf[a]] = cap[a] - fa;
-          _res_cap[_arc_idb[a]] = fa;
-          sup[_graph.source(a)] -= fa;
-          sup[_graph.target(a)] += fa;
-        }
-        for (NodeIt n(_graph); n != INVALID; ++n) {
-          _excess[_node_id[n]] = sup[n];
-        }
-        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
-          int u = _target[a];
-          int ra = _reverse[a];
-          _res_cap[a] = -_sum_supply + 1;
-          _res_cap[ra] = -_excess[u];
-          _cost[a] = 0;
-          _cost[ra] = 0;
-          _excess[u] = 0;
-        }
-      } else {
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          Value fa = flow[a];
-          _res_cap[_arc_idf[a]] = cap[a] - fa;
-          _res_cap[_arc_idb[a]] = fa;
-        }
-        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
-          int ra = _reverse[a];
-          _res_cap[a] = 1;
-          _res_cap[ra] = 0;
-          _cost[a] = 0;
-          _cost[ra] = 0;
-        }
-      }
-      
-      return OPTIMAL;
-    }
-
-    // Execute the algorithm and transform the results
-    void start(Method method) {
-      // Maximum path length for partial augment
-      const int MAX_PATH_LENGTH = 4;
-      
-      // Execute the algorithm
-      switch (method) {
-        case PUSH:
-          startPush();
-          break;
-        case AUGMENT:
-          startAugment();
-          break;
-        case PARTIAL_AUGMENT:
-          startAugment(MAX_PATH_LENGTH);
-          break;
-      }
-
-      // Compute node potentials for the original costs
-      _arc_vec.clear();
-      _cost_vec.clear();
-      for (int j = 0; j != _res_arc_num; ++j) {
-        if (_res_cap[j] > 0) {
-          _arc_vec.push_back(IntPair(_source[j], _target[j]));
-          _cost_vec.push_back(_scost[j]);
-        }
-      }
-      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
-
-      typename BellmanFord<StaticDigraph, LargeCostArcMap>
-        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);
-      bf.distMap(_pi_map);
-      bf.init(0);
-      bf.start();
-
-      // Handle non-zero lower bounds
-      if (_have_lower) {
-        int limit = _first_out[_root];
-        for (int j = 0; j != limit; ++j) {
-          if (!_forward[j]) _res_cap[j] += _lower[j];
-        }
-      }
-    }
-
-    /// Execute the algorithm performing augment and relabel operations
-    void startAugment(int max_length = std::numeric_limits<int>::max()) {
-      // Paramters for heuristics
-      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
-      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
-
-      // Perform cost scaling phases
-      IntVector pred_arc(_res_node_num);
-      std::vector<int> path_nodes;
-      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
-                                        1 : _epsilon / _alpha )
-      {
-        // "Early Termination" heuristic: use Bellman-Ford algorithm
-        // to check if the current flow is optimal
-        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
-          _arc_vec.clear();
-          _cost_vec.clear();
-          for (int j = 0; j != _res_arc_num; ++j) {
-            if (_res_cap[j] > 0) {
-              _arc_vec.push_back(IntPair(_source[j], _target[j]));
-              _cost_vec.push_back(_cost[j] + 1);
-            }
-          }
-          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
-
-          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
-          bf.init(0);
-          bool done = false;
-          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
-          for (int i = 0; i < K && !done; ++i)
-            done = bf.processNextWeakRound();
-          if (done) break;
-        }
-
-        // Saturate arcs not satisfying the optimality condition
-        for (int a = 0; a != _res_arc_num; ++a) {
-          if (_res_cap[a] > 0 &&
-              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
-            Value delta = _res_cap[a];
-            _excess[_source[a]] -= delta;
-            _excess[_target[a]] += delta;
-            _res_cap[a] = 0;
-            _res_cap[_reverse[a]] += delta;
-          }
-        }
-        
-        // Find active nodes (i.e. nodes with positive excess)
-        for (int u = 0; u != _res_node_num; ++u) {
-          if (_excess[u] > 0) _active_nodes.push_back(u);
-        }
-
-        // Initialize the next arcs
-        for (int u = 0; u != _res_node_num; ++u) {
-          _next_out[u] = _first_out[u];
-        }
-
-        // Perform partial augment and relabel operations
-        while (true) {
-          // Select an active node (FIFO selection)
-          while (_active_nodes.size() > 0 &&
-                 _excess[_active_nodes.front()] <= 0) {
-            _active_nodes.pop_front();
-          }
-          if (_active_nodes.size() == 0) break;
-          int start = _active_nodes.front();
-          path_nodes.clear();
-          path_nodes.push_back(start);
-
-          // Find an augmenting path from the start node
-          int tip = start;
-          while (_excess[tip] >= 0 &&
-                 int(path_nodes.size()) <= max_length) {
-            int u;
-            LargeCost min_red_cost, rc;
-            int last_out = _sum_supply < 0 ?
-              _first_out[tip+1] : _first_out[tip+1] - 1;
-            for (int a = _next_out[tip]; a != last_out; ++a) {
-              if (_res_cap[a] > 0 &&
-                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
-                u = _target[a];
-                pred_arc[u] = a;
-                _next_out[tip] = a;
-                tip = u;
-                path_nodes.push_back(tip);
-                goto next_step;
-              }
-            }
-
-            // Relabel tip node
-            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
-            for (int a = _first_out[tip]; a != last_out; ++a) {
-              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
-              if (_res_cap[a] > 0 && rc < min_red_cost) {
-                min_red_cost = rc;
-              }
-            }
-            _pi[tip] -= min_red_cost + _epsilon;
-
-            // Reset the next arc of tip
-            _next_out[tip] = _first_out[tip];
-
-            // Step back
-            if (tip != start) {
-              path_nodes.pop_back();
-              tip = path_nodes.back();
-            }
-
-          next_step: ;
-          }
-
-          // Augment along the found path (as much flow as possible)
-          Value delta;
-          int u, v = path_nodes.front(), pa;
-          for (int i = 1; i < int(path_nodes.size()); ++i) {
-            u = v;
-            v = path_nodes[i];
-            pa = pred_arc[v];
-            delta = std::min(_res_cap[pa], _excess[u]);
-            _res_cap[pa] -= delta;
-            _res_cap[_reverse[pa]] += delta;
-            _excess[u] -= delta;
-            _excess[v] += delta;
-            if (_excess[v] > 0 && _excess[v] <= delta)
-              _active_nodes.push_back(v);
-          }
-        }
-      }
-    }
-
-    /// Execute the algorithm performing push and relabel operations
-    void startPush() {
-      // Paramters for heuristics
-      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
-      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
-
-      // Perform cost scaling phases
-      BoolVector hyper(_res_node_num, false);
-      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
-                                        1 : _epsilon / _alpha )
-      {
-        // "Early Termination" heuristic: use Bellman-Ford algorithm
-        // to check if the current flow is optimal
-        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
-          _arc_vec.clear();
-          _cost_vec.clear();
-          for (int j = 0; j != _res_arc_num; ++j) {
-            if (_res_cap[j] > 0) {
-              _arc_vec.push_back(IntPair(_source[j], _target[j]));
-              _cost_vec.push_back(_cost[j] + 1);
-            }
-          }
-          _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
-
-          BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
-          bf.init(0);
-          bool done = false;
-          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
-          for (int i = 0; i < K && !done; ++i)
-            done = bf.processNextWeakRound();
-          if (done) break;
-        }
-
-        // Saturate arcs not satisfying the optimality condition
-        for (int a = 0; a != _res_arc_num; ++a) {
-          if (_res_cap[a] > 0 &&
-              _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
-            Value delta = _res_cap[a];
-            _excess[_source[a]] -= delta;
-            _excess[_target[a]] += delta;
-            _res_cap[a] = 0;
-            _res_cap[_reverse[a]] += delta;
-          }
-        }
-
-        // Find active nodes (i.e. nodes with positive excess)
-        for (int u = 0; u != _res_node_num; ++u) {
-          if (_excess[u] > 0) _active_nodes.push_back(u);
-        }
-
-        // Initialize the next arcs
-        for (int u = 0; u != _res_node_num; ++u) {
-          _next_out[u] = _first_out[u];
-        }
-
-        // Perform push and relabel operations
-        while (_active_nodes.size() > 0) {
-          LargeCost min_red_cost, rc;
-          Value delta;
-          int n, t, a, last_out = _res_arc_num;
-
-          // Select an active node (FIFO selection)
-        next_node:
-          n = _active_nodes.front();
-          last_out = _sum_supply < 0 ?
-            _first_out[n+1] : _first_out[n+1] - 1;
-
-          // Perform push operations if there are admissible arcs
-          if (_excess[n] > 0) {
-            for (a = _next_out[n]; a != last_out; ++a) {
-              if (_res_cap[a] > 0 &&
-                  _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
-                delta = std::min(_res_cap[a], _excess[n]);
-                t = _target[a];
-
-                // Push-look-ahead heuristic
-                Value ahead = -_excess[t];
-                int last_out_t = _sum_supply < 0 ?
-                  _first_out[t+1] : _first_out[t+1] - 1;
-                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
-                  if (_res_cap[ta] > 0 && 
-                      _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
-                    ahead += _res_cap[ta];
-                  if (ahead >= delta) break;
-                }
-                if (ahead < 0) ahead = 0;
-
-                // Push flow along the arc
-                if (ahead < delta) {
-                  _res_cap[a] -= ahead;
-                  _res_cap[_reverse[a]] += ahead;
-                  _excess[n] -= ahead;
-                  _excess[t] += ahead;
-                  _active_nodes.push_front(t);
-                  hyper[t] = true;
-                  _next_out[n] = a;
-                  goto next_node;
-                } else {
-                  _res_cap[a] -= delta;
-                  _res_cap[_reverse[a]] += delta;
-                  _excess[n] -= delta;
-                  _excess[t] += delta;
-                  if (_excess[t] > 0 && _excess[t] <= delta)
-                    _active_nodes.push_back(t);
-                }
-
-                if (_excess[n] == 0) {
-                  _next_out[n] = a;
-                  goto remove_nodes;
-                }
-              }
-            }
-            _next_out[n] = a;
-          }
-
-          // Relabel the node if it is still active (or hyper)
-          if (_excess[n] > 0 || hyper[n]) {
-            min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
-            for (int a = _first_out[n]; a != last_out; ++a) {
-              rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
-              if (_res_cap[a] > 0 && rc < min_red_cost) {
-                min_red_cost = rc;
-              }
-            }
-            _pi[n] -= min_red_cost + _epsilon;
-            hyper[n] = false;
-
-            // Reset the next arc
-            _next_out[n] = _first_out[n];
-          }
-        
-          // Remove nodes that are not active nor hyper
-        remove_nodes:
-          while ( _active_nodes.size() > 0 &&
-                  _excess[_active_nodes.front()] <= 0 &&
-                  !hyper[_active_nodes.front()] ) {
-            _active_nodes.pop_front();
-          }
-        }
-      }
-    }
-
-  }; //class CostScaling
-
-  ///@}
-
-} //namespace lemon
-
-#endif //LEMON_COST_SCALING_H
Index: mon/cycle_canceling.h
===================================================================
--- lemon/cycle_canceling.h	(revision 898)
+++ 	(revision )
@@ -1,1169 +1,0 @@
-/* -*- C++ -*-
- *
- * This file is a part of LEMON, a generic C++ optimization library
- *
- * Copyright (C) 2003-2008
- * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
- * (Egervary Research Group on Combinatorial Optimization, EGRES).
- *
- * Permission to use, modify and distribute this software is granted
- * provided that this copyright notice appears in all copies. For
- * precise terms see the accompanying LICENSE file.
- *
- * This software is provided "AS IS" with no warranty of any kind,
- * express or implied, and with no claim as to its suitability for any
- * purpose.
- *
- */
-
-#ifndef LEMON_CYCLE_CANCELING_H
-#define LEMON_CYCLE_CANCELING_H
-
-/// \ingroup min_cost_flow_algs
-/// \file
-/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
-
-#include <vector>
-#include <limits>
-
-#include <lemon/core.h>
-#include <lemon/maps.h>
-#include <lemon/path.h>
-#include <lemon/math.h>
-#include <lemon/static_graph.h>
-#include <lemon/adaptors.h>
-#include <lemon/circulation.h>
-#include <lemon/bellman_ford.h>
-#include <lemon/howard.h>
-
-namespace lemon {
-
-  /// \addtogroup min_cost_flow_algs
-  /// @{
-
-  /// \brief Implementation of cycle-canceling algorithms for
-  /// finding a \ref min_cost_flow "minimum cost flow".
-  ///
-  /// \ref CycleCanceling implements three different cycle-canceling
-  /// algorithms for finding a \ref min_cost_flow "minimum cost flow"
-  /// \ref amo93networkflows, \ref klein67primal,
-  /// \ref goldberg89cyclecanceling.
-  /// The most efficent one (both theoretically and practically)
-  /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
-  /// thus it is the default method.
-  /// It is strongly polynomial, but in practice, it is typically much
-  /// slower than the scaling algorithms and NetworkSimplex.
-  ///
-  /// Most of the parameters of the problem (except for the digraph)
-  /// can be given using separate functions, and the algorithm can be
-  /// executed using the \ref run() function. If some parameters are not
-  /// specified, then default values will be used.
-  ///
-  /// \tparam GR The digraph type the algorithm runs on.
-  /// \tparam V The number type used for flow amounts, capacity bounds
-  /// and supply values in the algorithm. By default, it is \c int.
-  /// \tparam C The number type used for costs and potentials in the
-  /// algorithm. By default, it is the same as \c V.
-  ///
-  /// \warning Both number types must be signed and all input data must
-  /// be integer.
-  /// \warning This algorithm does not support negative costs for such
-  /// arcs that have infinite upper bound.
-  ///
-  /// \note For more information about the three available methods,
-  /// see \ref Method.
-#ifdef DOXYGEN
-  template <typename GR, typename V, typename C>
-#else
-  template <typename GR, typename V = int, typename C = V>
-#endif
-  class CycleCanceling
-  {
-  public:
-
-    /// The type of the digraph
-    typedef GR Digraph;
-    /// The type of the flow amounts, capacity bounds and supply values
-    typedef V Value;
-    /// The type of the arc costs
-    typedef C Cost;
-
-  public:
-
-    /// \brief Problem type constants for the \c run() function.
-    ///
-    /// Enum type containing the problem type constants that can be
-    /// returned by the \ref run() function of the algorithm.
-    enum ProblemType {
-      /// The problem has no feasible solution (flow).
-      INFEASIBLE,
-      /// The problem has optimal solution (i.e. it is feasible and
-      /// bounded), and the algorithm has found optimal flow and node
-      /// potentials (primal and dual solutions).
-      OPTIMAL,
-      /// The digraph contains an arc of negative cost and infinite
-      /// upper bound. It means that the objective function is unbounded
-      /// on that arc, however, note that it could actually be bounded
-      /// over the feasible flows, but this algroithm cannot handle
-      /// these cases.
-      UNBOUNDED
-    };
-
-    /// \brief Constants for selecting the used method.
-    ///
-    /// Enum type containing constants for selecting the used method
-    /// for the \ref run() function.
-    ///
-    /// \ref CycleCanceling provides three different cycle-canceling
-    /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
-    /// is used, which proved to be the most efficient and the most robust
-    /// on various test inputs.
-    /// However, the other methods can be selected using the \ref run()
-    /// function with the proper parameter.
-    enum Method {
-      /// A simple cycle-canceling method, which uses the
-      /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
-      /// number for detecting negative cycles in the residual network.
-      SIMPLE_CYCLE_CANCELING,
-      /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
-      /// well-known strongly polynomial method
-      /// \ref goldberg89cyclecanceling. It improves along a
-      /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
-      /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
-      MINIMUM_MEAN_CYCLE_CANCELING,
-      /// The "Cancel And Tighten" algorithm, which can be viewed as an
-      /// improved version of the previous method
-      /// \ref goldberg89cyclecanceling.
-      /// It is faster both in theory and in practice, its running time
-      /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
-      CANCEL_AND_TIGHTEN
-    };
-
-  private:
-
-    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
-    
-    typedef std::vector<int> IntVector;
-    typedef std::vector<char> CharVector;
-    typedef std::vector<double> DoubleVector;
-    typedef std::vector<Value> ValueVector;
-    typedef std::vector<Cost> CostVector;
-
-  private:
-  
-    template <typename KT, typename VT>
-    class StaticVectorMap {
-    public:
-      typedef KT Key;
-      typedef VT Value;
-      
-      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
-      
-      const Value& operator[](const Key& key) const {
-        return _v[StaticDigraph::id(key)];
-      }
-
-      Value& operator[](const Key& key) {
-        return _v[StaticDigraph::id(key)];
-      }
-      
-      void set(const Key& key, const Value& val) {
-        _v[StaticDigraph::id(key)] = val;
-      }
-
-    private:
-      std::vector<Value>& _v;
-    };
-
-    typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap;
-    typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap;
-
-  private:
-
-
-    // Data related to the underlying digraph
-    const GR &_graph;
-    int _node_num;
-    int _arc_num;
-    int _res_node_num;
-    int _res_arc_num;
-    int _root;
-
-    // Parameters of the problem
-    bool _have_lower;
-    Value _sum_supply;
-
-    // Data structures for storing the digraph
-    IntNodeMap _node_id;
-    IntArcMap _arc_idf;
-    IntArcMap _arc_idb;
-    IntVector _first_out;
-    CharVector _forward;
-    IntVector _source;
-    IntVector _target;
-    IntVector _reverse;
-
-    // Node and arc data
-    ValueVector _lower;
-    ValueVector _upper;
-    CostVector _cost;
-    ValueVector _supply;
-
-    ValueVector _res_cap;
-    CostVector _pi;
-
-    // Data for a StaticDigraph structure
-    typedef std::pair<int, int> IntPair;
-    StaticDigraph _sgr;
-    std::vector<IntPair> _arc_vec;
-    std::vector<Cost> _cost_vec;
-    IntVector _id_vec;
-    CostArcMap _cost_map;
-    CostNodeMap _pi_map;
-  
-  public:
-  
-    /// \brief Constant for infinite upper bounds (capacities).
-    ///
-    /// Constant for infinite upper bounds (capacities).
-    /// It is \c std::numeric_limits<Value>::infinity() if available,
-    /// \c std::numeric_limits<Value>::max() otherwise.
-    const Value INF;
-
-  public:
-
-    /// \brief Constructor.
-    ///
-    /// The constructor of the class.
-    ///
-    /// \param graph The digraph the algorithm runs on.
-    CycleCanceling(const GR& graph) :
-      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
-      _cost_map(_cost_vec), _pi_map(_pi),
-      INF(std::numeric_limits<Value>::has_infinity ?
-          std::numeric_limits<Value>::infinity() :
-          std::numeric_limits<Value>::max())
-    {
-      // Check the number types
-      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
-        "The flow type of CycleCanceling must be signed");
-      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
-        "The cost type of CycleCanceling must be signed");
-
-      // Reset data structures
-      reset();
-    }
-
-    /// \name Parameters
-    /// The parameters of the algorithm can be specified using these
-    /// functions.
-
-    /// @{
-
-    /// \brief Set the lower bounds on the arcs.
-    ///
-    /// This function sets the lower bounds on the arcs.
-    /// If it is not used before calling \ref run(), the lower bounds
-    /// will be set to zero on all arcs.
-    ///
-    /// \param map An arc map storing the lower bounds.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template <typename LowerMap>
-    CycleCanceling& lowerMap(const LowerMap& map) {
-      _have_lower = true;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _lower[_arc_idf[a]] = map[a];
-        _lower[_arc_idb[a]] = map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the upper bounds (capacities) on the arcs.
-    ///
-    /// This function sets the upper bounds (capacities) on the arcs.
-    /// If it is not used before calling \ref run(), the upper bounds
-    /// will be set to \ref INF on all arcs (i.e. the flow value will be
-    /// unbounded from above).
-    ///
-    /// \param map An arc map storing the upper bounds.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename UpperMap>
-    CycleCanceling& upperMap(const UpperMap& map) {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _upper[_arc_idf[a]] = map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the costs of the arcs.
-    ///
-    /// This function sets the costs of the arcs.
-    /// If it is not used before calling \ref run(), the costs
-    /// will be set to \c 1 on all arcs.
-    ///
-    /// \param map An arc map storing the costs.
-    /// Its \c Value type must be convertible to the \c Cost type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename CostMap>
-    CycleCanceling& costMap(const CostMap& map) {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        _cost[_arc_idf[a]] =  map[a];
-        _cost[_arc_idb[a]] = -map[a];
-      }
-      return *this;
-    }
-
-    /// \brief Set the supply values of the nodes.
-    ///
-    /// This function sets the supply values of the nodes.
-    /// If neither this function nor \ref stSupply() is used before
-    /// calling \ref run(), the supply of each node will be set to zero.
-    ///
-    /// \param map A node map storing the supply values.
-    /// Its \c Value type must be convertible to the \c Value type
-    /// of the algorithm.
-    ///
-    /// \return <tt>(*this)</tt>
-    template<typename SupplyMap>
-    CycleCanceling& supplyMap(const SupplyMap& map) {
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        _supply[_node_id[n]] = map[n];
-      }
-      return *this;
-    }
-
-    /// \brief Set single source and target nodes and a supply value.
-    ///
-    /// This function sets a single source node and a single target node
-    /// and the required flow value.
-    /// If neither this function nor \ref supplyMap() is used before
-    /// calling \ref run(), the supply of each node will be set to zero.
-    ///
-    /// Using this function has the same effect as using \ref supplyMap()
-    /// with such a map in which \c k is assigned to \c s, \c -k is
-    /// assigned to \c t and all other nodes have zero supply value.
-    ///
-    /// \param s The source node.
-    /// \param t The target node.
-    /// \param k The required amount of flow from node \c s to node \c t
-    /// (i.e. the supply of \c s and the demand of \c t).
-    ///
-    /// \return <tt>(*this)</tt>
-    CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
-      for (int i = 0; i != _res_node_num; ++i) {
-        _supply[i] = 0;
-      }
-      _supply[_node_id[s]] =  k;
-      _supply[_node_id[t]] = -k;
-      return *this;
-    }
-    
-    /// @}
-
-    /// \name Execution control
-    /// The algorithm can be executed using \ref run().
-
-    /// @{
-
-    /// \brief Run the algorithm.
-    ///
-    /// This function runs the algorithm.
-    /// The paramters can be specified using functions \ref lowerMap(),
-    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
-    /// For example,
-    /// \code
-    ///   CycleCanceling<ListDigraph> cc(graph);
-    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
-    ///     .supplyMap(sup).run();
-    /// \endcode
-    ///
-    /// This function can be called more than once. All the given parameters
-    /// are kept for the next call, unless \ref resetParams() or \ref reset()
-    /// is used, thus only the modified parameters have to be set again.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class (or the last \ref reset() call), then the \ref reset()
-    /// function must be called.
-    ///
-    /// \param method The cycle-canceling method that will be used.
-    /// For more information, see \ref Method.
-    ///
-    /// \return \c INFEASIBLE if no feasible flow exists,
-    /// \n \c OPTIMAL if the problem has optimal solution
-    /// (i.e. it is feasible and bounded), and the algorithm has found
-    /// optimal flow and node potentials (primal and dual solutions),
-    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
-    /// and infinite upper bound. It means that the objective function
-    /// is unbounded on that arc, however, note that it could actually be
-    /// bounded over the feasible flows, but this algroithm cannot handle
-    /// these cases.
-    ///
-    /// \see ProblemType, Method
-    /// \see resetParams(), reset()
-    ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
-      ProblemType pt = init();
-      if (pt != OPTIMAL) return pt;
-      start(method);
-      return OPTIMAL;
-    }
-
-    /// \brief Reset all the parameters that have been given before.
-    ///
-    /// This function resets all the paramaters that have been given
-    /// before using functions \ref lowerMap(), \ref upperMap(),
-    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
-    ///
-    /// It is useful for multiple \ref run() calls. Basically, all the given
-    /// parameters are kept for the next \ref run() call, unless
-    /// \ref resetParams() or \ref reset() is used.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class or the last \ref reset() call, then the \ref reset()
-    /// function must be used, otherwise \ref resetParams() is sufficient.
-    ///
-    /// For example,
-    /// \code
-    ///   CycleCanceling<ListDigraph> cs(graph);
-    ///
-    ///   // First run
-    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
-    ///     .supplyMap(sup).run();
-    ///
-    ///   // Run again with modified cost map (resetParams() is not called,
-    ///   // so only the cost map have to be set again)
-    ///   cost[e] += 100;
-    ///   cc.costMap(cost).run();
-    ///
-    ///   // Run again from scratch using resetParams()
-    ///   // (the lower bounds will be set to zero on all arcs)
-    ///   cc.resetParams();
-    ///   cc.upperMap(capacity).costMap(cost)
-    ///     .supplyMap(sup).run();
-    /// \endcode
-    ///
-    /// \return <tt>(*this)</tt>
-    ///
-    /// \see reset(), run()
-    CycleCanceling& resetParams() {
-      for (int i = 0; i != _res_node_num; ++i) {
-        _supply[i] = 0;
-      }
-      int limit = _first_out[_root];
-      for (int j = 0; j != limit; ++j) {
-        _lower[j] = 0;
-        _upper[j] = INF;
-        _cost[j] = _forward[j] ? 1 : -1;
-      }
-      for (int j = limit; j != _res_arc_num; ++j) {
-        _lower[j] = 0;
-        _upper[j] = INF;
-        _cost[j] = 0;
-        _cost[_reverse[j]] = 0;
-      }      
-      _have_lower = false;
-      return *this;
-    }
-
-    /// \brief Reset the internal data structures and all the parameters
-    /// that have been given before.
-    ///
-    /// This function resets the internal data structures and all the
-    /// paramaters that have been given before using functions \ref lowerMap(),
-    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
-    ///
-    /// It is useful for multiple \ref run() calls. Basically, all the given
-    /// parameters are kept for the next \ref run() call, unless
-    /// \ref resetParams() or \ref reset() is used.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class or the last \ref reset() call, then the \ref reset()
-    /// function must be used, otherwise \ref resetParams() is sufficient.
-    ///
-    /// See \ref resetParams() for examples.
-    ///
-    /// \return <tt>(*this)</tt>
-    ///
-    /// \see resetParams(), run()
-    CycleCanceling& reset() {
-      // Resize vectors
-      _node_num = countNodes(_graph);
-      _arc_num = countArcs(_graph);
-      _res_node_num = _node_num + 1;
-      _res_arc_num = 2 * (_arc_num + _node_num);
-      _root = _node_num;
-
-      _first_out.resize(_res_node_num + 1);
-      _forward.resize(_res_arc_num);
-      _source.resize(_res_arc_num);
-      _target.resize(_res_arc_num);
-      _reverse.resize(_res_arc_num);
-
-      _lower.resize(_res_arc_num);
-      _upper.resize(_res_arc_num);
-      _cost.resize(_res_arc_num);
-      _supply.resize(_res_node_num);
-      
-      _res_cap.resize(_res_arc_num);
-      _pi.resize(_res_node_num);
-
-      _arc_vec.reserve(_res_arc_num);
-      _cost_vec.reserve(_res_arc_num);
-      _id_vec.reserve(_res_arc_num);
-
-      // Copy the graph
-      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
-      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-        _node_id[n] = i;
-      }
-      i = 0;
-      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-        _first_out[i] = j;
-        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
-          _arc_idf[a] = j;
-          _forward[j] = true;
-          _source[j] = i;
-          _target[j] = _node_id[_graph.runningNode(a)];
-        }
-        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
-          _arc_idb[a] = j;
-          _forward[j] = false;
-          _source[j] = i;
-          _target[j] = _node_id[_graph.runningNode(a)];
-        }
-        _forward[j] = false;
-        _source[j] = i;
-        _target[j] = _root;
-        _reverse[j] = k;
-        _forward[k] = true;
-        _source[k] = _root;
-        _target[k] = i;
-        _reverse[k] = j;
-        ++j; ++k;
-      }
-      _first_out[i] = j;
-      _first_out[_res_node_num] = k;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        int fi = _arc_idf[a];
-        int bi = _arc_idb[a];
-        _reverse[fi] = bi;
-        _reverse[bi] = fi;
-      }
-      
-      // Reset parameters
-      resetParams();
-      return *this;
-    }
-
-    /// @}
-
-    /// \name Query Functions
-    /// The results of the algorithm can be obtained using these
-    /// functions.\n
-    /// The \ref run() function must be called before using them.
-
-    /// @{
-
-    /// \brief Return the total cost of the found flow.
-    ///
-    /// This function returns the total cost of the found flow.
-    /// Its complexity is O(e).
-    ///
-    /// \note The return type of the function can be specified as a
-    /// template parameter. For example,
-    /// \code
-    ///   cc.totalCost<double>();
-    /// \endcode
-    /// It is useful if the total cost cannot be stored in the \c Cost
-    /// type of the algorithm, which is the default return type of the
-    /// function.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename Number>
-    Number totalCost() const {
-      Number c = 0;
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        int i = _arc_idb[a];
-        c += static_cast<Number>(_res_cap[i]) *
-             (-static_cast<Number>(_cost[i]));
-      }
-      return c;
-    }
-
-#ifndef DOXYGEN
-    Cost totalCost() const {
-      return totalCost<Cost>();
-    }
-#endif
-
-    /// \brief Return the flow on the given arc.
-    ///
-    /// This function returns the flow on the given arc.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    Value flow(const Arc& a) const {
-      return _res_cap[_arc_idb[a]];
-    }
-
-    /// \brief Return the flow map (the primal solution).
-    ///
-    /// This function copies the flow value on each arc into the given
-    /// map. The \c Value type of the algorithm must be convertible to
-    /// the \c Value type of the map.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename FlowMap>
-    void flowMap(FlowMap &map) const {
-      for (ArcIt a(_graph); a != INVALID; ++a) {
-        map.set(a, _res_cap[_arc_idb[a]]);
-      }
-    }
-
-    /// \brief Return the potential (dual value) of the given node.
-    ///
-    /// This function returns the potential (dual value) of the
-    /// given node.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    Cost potential(const Node& n) const {
-      return static_cast<Cost>(_pi[_node_id[n]]);
-    }
-
-    /// \brief Return the potential map (the dual solution).
-    ///
-    /// This function copies the potential (dual value) of each node
-    /// into the given map.
-    /// The \c Cost type of the algorithm must be convertible to the
-    /// \c Value type of the map.
-    ///
-    /// \pre \ref run() must be called before using this function.
-    template <typename PotentialMap>
-    void potentialMap(PotentialMap &map) const {
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
-      }
-    }
-
-    /// @}
-
-  private:
-
-    // Initialize the algorithm
-    ProblemType init() {
-      if (_res_node_num <= 1) return INFEASIBLE;
-
-      // Check the sum of supply values
-      _sum_supply = 0;
-      for (int i = 0; i != _root; ++i) {
-        _sum_supply += _supply[i];
-      }
-      if (_sum_supply > 0) return INFEASIBLE;
-      
-
-      // Initialize vectors
-      for (int i = 0; i != _res_node_num; ++i) {
-        _pi[i] = 0;
-      }
-      ValueVector excess(_supply);
-      
-      // Remove infinite upper bounds and check negative arcs
-      const Value MAX = std::numeric_limits<Value>::max();
-      int last_out;
-      if (_have_lower) {
-        for (int i = 0; i != _root; ++i) {
-          last_out = _first_out[i+1];
-          for (int j = _first_out[i]; j != last_out; ++j) {
-            if (_forward[j]) {
-              Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
-              if (c >= MAX) return UNBOUNDED;
-              excess[i] -= c;
-              excess[_target[j]] += c;
-            }
-          }
-        }
-      } else {
-        for (int i = 0; i != _root; ++i) {
-          last_out = _first_out[i+1];
-          for (int j = _first_out[i]; j != last_out; ++j) {
-            if (_forward[j] && _cost[j] < 0) {
-              Value c = _upper[j];
-              if (c >= MAX) return UNBOUNDED;
-              excess[i] -= c;
-              excess[_target[j]] += c;
-            }
-          }
-        }
-      }
-      Value ex, max_cap = 0;
-      for (int i = 0; i != _res_node_num; ++i) {
-        ex = excess[i];
-        if (ex < 0) max_cap -= ex;
-      }
-      for (int j = 0; j != _res_arc_num; ++j) {
-        if (_upper[j] >= MAX) _upper[j] = max_cap;
-      }
-
-      // Initialize maps for Circulation and remove non-zero lower bounds
-      ConstMap<Arc, Value> low(0);
-      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
-      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
-      ValueArcMap cap(_graph), flow(_graph);
-      ValueNodeMap sup(_graph);
-      for (NodeIt n(_graph); n != INVALID; ++n) {
-        sup[n] = _supply[_node_id[n]];
-      }
-      if (_have_lower) {
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          int j = _arc_idf[a];
-          Value c = _lower[j];
-          cap[a] = _upper[j] - c;
-          sup[_graph.source(a)] -= c;
-          sup[_graph.target(a)] += c;
-        }
-      } else {
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          cap[a] = _upper[_arc_idf[a]];
-        }
-      }
-
-      // Find a feasible flow using Circulation
-      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
-        circ(_graph, low, cap, sup);
-      if (!circ.flowMap(flow).run()) return INFEASIBLE;
-
-      // Set residual capacities and handle GEQ supply type
-      if (_sum_supply < 0) {
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          Value fa = flow[a];
-          _res_cap[_arc_idf[a]] = cap[a] - fa;
-          _res_cap[_arc_idb[a]] = fa;
-          sup[_graph.source(a)] -= fa;
-          sup[_graph.target(a)] += fa;
-        }
-        for (NodeIt n(_graph); n != INVALID; ++n) {
-          excess[_node_id[n]] = sup[n];
-        }
-        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
-          int u = _target[a];
-          int ra = _reverse[a];
-          _res_cap[a] = -_sum_supply + 1;
-          _res_cap[ra] = -excess[u];
-          _cost[a] = 0;
-          _cost[ra] = 0;
-        }
-      } else {
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          Value fa = flow[a];
-          _res_cap[_arc_idf[a]] = cap[a] - fa;
-          _res_cap[_arc_idb[a]] = fa;
-        }
-        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
-          int ra = _reverse[a];
-          _res_cap[a] = 1;
-          _res_cap[ra] = 0;
-          _cost[a] = 0;
-          _cost[ra] = 0;
-        }
-      }
-      
-      return OPTIMAL;
-    }
-    
-    // Build a StaticDigraph structure containing the current
-    // residual network
-    void buildResidualNetwork() {
-      _arc_vec.clear();
-      _cost_vec.clear();
-      _id_vec.clear();
-      for (int j = 0; j != _res_arc_num; ++j) {
-        if (_res_cap[j] > 0) {
-          _arc_vec.push_back(IntPair(_source[j], _target[j]));
-          _cost_vec.push_back(_cost[j]);
-          _id_vec.push_back(j);
-        }
-      }
-      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
-    }
-
-    // Execute the algorithm and transform the results
-    void start(Method method) {
-      // Execute the algorithm
-      switch (method) {
-        case SIMPLE_CYCLE_CANCELING:
-          startSimpleCycleCanceling();
-          break;
-        case MINIMUM_MEAN_CYCLE_CANCELING:
-          startMinMeanCycleCanceling();
-          break;
-        case CANCEL_AND_TIGHTEN:
-          startCancelAndTighten();
-          break;
-      }
-
-      // Compute node potentials
-      if (method != SIMPLE_CYCLE_CANCELING) {
-        buildResidualNetwork();
-        typename BellmanFord<StaticDigraph, CostArcMap>
-          ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map);
-        bf.distMap(_pi_map);
-        bf.init(0);
-        bf.start();
-      }
-
-      // Handle non-zero lower bounds
-      if (_have_lower) {
-        int limit = _first_out[_root];
-        for (int j = 0; j != limit; ++j) {
-          if (!_forward[j]) _res_cap[j] += _lower[j];
-        }
-      }
-    }
-
-    // Execute the "Simple Cycle Canceling" method
-    void startSimpleCycleCanceling() {
-      // Constants for computing the iteration limits
-      const int BF_FIRST_LIMIT  = 2;
-      const double BF_LIMIT_FACTOR = 1.5;
-      
-      typedef StaticVectorMap<StaticDigraph::Arc, Value> FilterMap;
-      typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;
-      typedef StaticVectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;
-      typedef typename BellmanFord<ResDigraph, CostArcMap>
-        ::template SetDistMap<CostNodeMap>
-        ::template SetPredMap<PredMap>::Create BF;
-      
-      // Build the residual network
-      _arc_vec.clear();
-      _cost_vec.clear();
-      for (int j = 0; j != _res_arc_num; ++j) {
-        _arc_vec.push_back(IntPair(_source[j], _target[j]));
-        _cost_vec.push_back(_cost[j]);
-      }
-      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
-
-      FilterMap filter_map(_res_cap);
-      ResDigraph rgr(_sgr, filter_map);
-      std::vector<int> cycle;
-      std::vector<StaticDigraph::Arc> pred(_res_arc_num);
-      PredMap pred_map(pred);
-      BF bf(rgr, _cost_map);
-      bf.distMap(_pi_map).predMap(pred_map);
-
-      int length_bound = BF_FIRST_LIMIT;
-      bool optimal = false;
-      while (!optimal) {
-        bf.init(0);
-        int iter_num = 0;
-        bool cycle_found = false;
-        while (!cycle_found) {
-          // Perform some iterations of the Bellman-Ford algorithm
-          int curr_iter_num = iter_num + length_bound <= _node_num ?
-            length_bound : _node_num - iter_num;
-          iter_num += curr_iter_num;
-          int real_iter_num = curr_iter_num;
-          for (int i = 0; i < curr_iter_num; ++i) {
-            if (bf.processNextWeakRound()) {
-              real_iter_num = i;
-              break;
-            }
-          }
-          if (real_iter_num < curr_iter_num) {
-            // Optimal flow is found
-            optimal = true;
-            break;
-          } else {
-            // Search for node disjoint negative cycles
-            std::vector<int> state(_res_node_num, 0);
-            int id = 0;
-            for (int u = 0; u != _res_node_num; ++u) {
-              if (state[u] != 0) continue;
-              ++id;
-              int v = u;
-              for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ?
-                   -1 : rgr.id(rgr.source(pred[v]))) {
-                state[v] = id;
-              }
-              if (v != -1 && state[v] == id) {
-                // A negative cycle is found
-                cycle_found = true;
-                cycle.clear();
-                StaticDigraph::Arc a = pred[v];
-                Value d, delta = _res_cap[rgr.id(a)];
-                cycle.push_back(rgr.id(a));
-                while (rgr.id(rgr.source(a)) != v) {
-                  a = pred_map[rgr.source(a)];
-                  d = _res_cap[rgr.id(a)];
-                  if (d < delta) delta = d;
-                  cycle.push_back(rgr.id(a));
-                }
-
-                // Augment along the cycle
-                for (int i = 0; i < int(cycle.size()); ++i) {
-                  int j = cycle[i];
-                  _res_cap[j] -= delta;
-                  _res_cap[_reverse[j]] += delta;
-                }
-              }
-            }
-          }
-
-          // Increase iteration limit if no cycle is found
-          if (!cycle_found) {
-            length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR);
-          }
-        }
-      }
-    }
-
-    // Execute the "Minimum Mean Cycle Canceling" method
-    void startMinMeanCycleCanceling() {
-      typedef SimplePath<StaticDigraph> SPath;
-      typedef typename SPath::ArcIt SPathArcIt;
-      typedef typename Howard<StaticDigraph, CostArcMap>
-        ::template SetPath<SPath>::Create MMC;
-      
-      SPath cycle;
-      MMC mmc(_sgr, _cost_map);
-      mmc.cycle(cycle);
-      buildResidualNetwork();
-      while (mmc.findMinMean() && mmc.cycleLength() < 0) {
-        // Find the cycle
-        mmc.findCycle();
-
-        // Compute delta value
-        Value delta = INF;
-        for (SPathArcIt a(cycle); a != INVALID; ++a) {
-          Value d = _res_cap[_id_vec[_sgr.id(a)]];
-          if (d < delta) delta = d;
-        }
-
-        // Augment along the cycle
-        for (SPathArcIt a(cycle); a != INVALID; ++a) {
-          int j = _id_vec[_sgr.id(a)];
-          _res_cap[j] -= delta;
-          _res_cap[_reverse[j]] += delta;
-        }
-
-        // Rebuild the residual network        
-        buildResidualNetwork();
-      }
-    }
-
-    // Execute the "Cancel And Tighten" method
-    void startCancelAndTighten() {
-      // Constants for the min mean cycle computations
-      const double LIMIT_FACTOR = 1.0;
-      const int MIN_LIMIT = 5;
-
-      // Contruct auxiliary data vectors
-      DoubleVector pi(_res_node_num, 0.0);
-      IntVector level(_res_node_num);
-      CharVector reached(_res_node_num);
-      CharVector processed(_res_node_num);
-      IntVector pred_node(_res_node_num);
-      IntVector pred_arc(_res_node_num);
-      std::vector<int> stack(_res_node_num);
-      std::vector<int> proc_vector(_res_node_num);
-
-      // Initialize epsilon
-      double epsilon = 0;
-      for (int a = 0; a != _res_arc_num; ++a) {
-        if (_res_cap[a] > 0 && -_cost[a] > epsilon)
-          epsilon = -_cost[a];
-      }
-
-      // Start phases
-      Tolerance<double> tol;
-      tol.epsilon(1e-6);
-      int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));
-      if (limit < MIN_LIMIT) limit = MIN_LIMIT;
-      int iter = limit;
-      while (epsilon * _res_node_num >= 1) {
-        // Find and cancel cycles in the admissible network using DFS
-        for (int u = 0; u != _res_node_num; ++u) {
-          reached[u] = false;
-          processed[u] = false;
-        }
-        int stack_head = -1;
-        int proc_head = -1;
-        for (int start = 0; start != _res_node_num; ++start) {
-          if (reached[start]) continue;
-
-          // New start node
-          reached[start] = true;
-          pred_arc[start] = -1;
-          pred_node[start] = -1;
-
-          // Find the first admissible outgoing arc
-          double p = pi[start];
-          int a = _first_out[start];
-          int last_out = _first_out[start+1];
-          for (; a != last_out && (_res_cap[a] == 0 ||
-               !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
-          if (a == last_out) {
-            processed[start] = true;
-            proc_vector[++proc_head] = start;
-            continue;
-          }
-          stack[++stack_head] = a;
-
-          while (stack_head >= 0) {
-            int sa = stack[stack_head];
-            int u = _source[sa];
-            int v = _target[sa];
-
-            if (!reached[v]) {
-              // A new node is reached
-              reached[v] = true;
-              pred_node[v] = u;
-              pred_arc[v] = sa;
-              p = pi[v];
-              a = _first_out[v];
-              last_out = _first_out[v+1];
-              for (; a != last_out && (_res_cap[a] == 0 ||
-                   !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
-              stack[++stack_head] = a == last_out ? -1 : a;
-            } else {
-              if (!processed[v]) {
-                // A cycle is found
-                int n, w = u;
-                Value d, delta = _res_cap[sa];
-                for (n = u; n != v; n = pred_node[n]) {
-                  d = _res_cap[pred_arc[n]];
-                  if (d <= delta) {
-                    delta = d;
-                    w = pred_node[n];
-                  }
-                }
-
-                // Augment along the cycle
-                _res_cap[sa] -= delta;
-                _res_cap[_reverse[sa]] += delta;
-                for (n = u; n != v; n = pred_node[n]) {
-                  int pa = pred_arc[n];
-                  _res_cap[pa] -= delta;
-                  _res_cap[_reverse[pa]] += delta;
-                }
-                for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
-                  --stack_head;
-                  reached[n] = false;
-                }
-                u = w;
-              }
-              v = u;
-
-              // Find the next admissible outgoing arc
-              p = pi[v];
-              a = stack[stack_head] + 1;
-              last_out = _first_out[v+1];
-              for (; a != last_out && (_res_cap[a] == 0 ||
-                   !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
-              stack[stack_head] = a == last_out ? -1 : a;
-            }
-
-            while (stack_head >= 0 && stack[stack_head] == -1) {
-              processed[v] = true;
-              proc_vector[++proc_head] = v;
-              if (--stack_head >= 0) {
-                // Find the next admissible outgoing arc
-                v = _source[stack[stack_head]];
-                p = pi[v];
-                a = stack[stack_head] + 1;
-                last_out = _first_out[v+1];
-                for (; a != last_out && (_res_cap[a] == 0 ||
-                     !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
-                stack[stack_head] = a == last_out ? -1 : a;
-              }
-            }
-          }
-        }
-
-        // Tighten potentials and epsilon
-        if (--iter > 0) {
-          for (int u = 0; u != _res_node_num; ++u) {
-            level[u] = 0;
-          }
-          for (int i = proc_head; i > 0; --i) {
-            int u = proc_vector[i];
-            double p = pi[u];
-            int l = level[u] + 1;
-            int last_out = _first_out[u+1];
-            for (int a = _first_out[u]; a != last_out; ++a) {
-              int v = _target[a];
-              if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) &&
-                  l > level[v]) level[v] = l;
-            }
-          }
-
-          // Modify potentials
-          double q = std::numeric_limits<double>::max();
-          for (int u = 0; u != _res_node_num; ++u) {
-            int lu = level[u];
-            double p, pu = pi[u];
-            int last_out = _first_out[u+1];
-            for (int a = _first_out[u]; a != last_out; ++a) {
-              if (_res_cap[a] == 0) continue;
-              int v = _target[a];
-              int ld = lu - level[v];
-              if (ld > 0) {
-                p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1);
-                if (p < q) q = p;
-              }
-            }
-          }
-          for (int u = 0; u != _res_node_num; ++u) {
-            pi[u] -= q * level[u];
-          }
-
-          // Modify epsilon
-          epsilon = 0;
-          for (int u = 0; u != _res_node_num; ++u) {
-            double curr, pu = pi[u];
-            int last_out = _first_out[u+1];
-            for (int a = _first_out[u]; a != last_out; ++a) {
-              if (_res_cap[a] == 0) continue;
-              curr = _cost[a] + pu - pi[_target[a]];
-              if (-curr > epsilon) epsilon = -curr;
-            }
-          }
-        } else {
-          typedef Howard<StaticDigraph, CostArcMap> MMC;
-          typedef typename BellmanFord<StaticDigraph, CostArcMap>
-            ::template SetDistMap<CostNodeMap>::Create BF;
-
-          // Set epsilon to the minimum cycle mean
-          buildResidualNetwork();
-          MMC mmc(_sgr, _cost_map);
-          mmc.findMinMean();
-          epsilon = -mmc.cycleMean();
-          Cost cycle_cost = mmc.cycleLength();
-          int cycle_size = mmc.cycleArcNum();
-          
-          // Compute feasible potentials for the current epsilon
-          for (int i = 0; i != int(_cost_vec.size()); ++i) {
-            _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;
-          }
-          BF bf(_sgr, _cost_map);
-          bf.distMap(_pi_map);
-          bf.init(0);
-          bf.start();
-          for (int u = 0; u != _res_node_num; ++u) {
-            pi[u] = static_cast<double>(_pi[u]) / cycle_size;
-          }
-        
-          iter = limit;
-        }
-      }
-    }
-
-  }; //class CycleCanceling
-
-  ///@}
-
-} //namespace lemon
-
-#endif //LEMON_CYCLE_CANCELING_H
Index: lemon/dfs.h
===================================================================
--- lemon/dfs.h	(revision 891)
+++ lemon/dfs.h	(revision 835)
@@ -122,9 +122,4 @@
   ///\tparam GR The type of the digraph the algorithm runs on.
   ///The default type is \ref ListDigraph.
-  ///\tparam TR The traits class that defines various types used by the
-  ///algorithm. By default, it is \ref DfsDefaultTraits
-  ///"DfsDefaultTraits<GR>".
-  ///In most cases, this parameter should not be set directly,
-  ///consider to use the named template parameters instead.
 #ifdef DOXYGEN
   template <typename GR,
@@ -893,7 +888,4 @@
   /// This class should only be used through the \ref dfs() function,
   /// which makes it easier to use the algorithm.
-  ///
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm.
   template<class TR>
   class DfsWizard : public TR
@@ -1246,9 +1238,9 @@
   /// does not observe the DFS events. If you want to observe the DFS
   /// events, you should implement your own visitor class.
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref DfsVisitDefaultTraits
-  /// "DfsVisitDefaultTraits<GR>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
+  /// \tparam TR Traits class to set various data types used by the
+  /// algorithm. The default traits class is
+  /// \ref DfsVisitDefaultTraits "DfsVisitDefaultTraits<GR>".
+  /// See \ref DfsVisitDefaultTraits for the documentation of
+  /// a DFS visit traits class.
 #ifdef DOXYGEN
   template <typename GR, typename VS, typename TR>
Index: lemon/dijkstra.h
===================================================================
--- lemon/dijkstra.h	(revision 891)
+++ lemon/dijkstra.h	(revision 835)
@@ -193,9 +193,4 @@
   ///it is necessary. The default map type is \ref
   ///concepts::Digraph::ArcMap "GR::ArcMap<int>".
-  ///\tparam TR The traits class that defines various types used by the
-  ///algorithm. By default, it is \ref DijkstraDefaultTraits
-  ///"DijkstraDefaultTraits<GR, LEN>".
-  ///In most cases, this parameter should not be set directly,
-  ///consider to use the named template parameters instead.
 #ifdef DOXYGEN
   template <typename GR, typename LEN, typename TR>
@@ -1098,7 +1093,4 @@
   /// This class should only be used through the \ref dijkstra() function,
   /// which makes it easier to use the algorithm.
-  ///
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm.
   template<class TR>
   class DijkstraWizard : public TR
Index: lemon/glpk.h
===================================================================
--- lemon/glpk.h	(revision 902)
+++ lemon/glpk.h	(revision 793)
@@ -26,26 +26,14 @@
 #include <lemon/lp_base.h>
 
+// forward declaration
+#if !defined _GLP_PROB && !defined GLP_PROB
+#define _GLP_PROB
+#define GLP_PROB
+typedef struct { double _opaque_prob; } glp_prob;
+/* LP/MIP problem object */
+#endif
+
 namespace lemon {
 
-  namespace _solver_bits {
-    class VoidPtr {
-    private:
-      void *_ptr;      
-    public:
-      VoidPtr() : _ptr(0) {}
-
-      template <typename T>
-      VoidPtr(T* ptr) : _ptr(reinterpret_cast<void*>(ptr)) {}
-
-      template <typename T>
-      VoidPtr& operator=(T* ptr) { 
-        _ptr = reinterpret_cast<void*>(ptr); 
-        return *this;
-      }
-
-      template <typename T>
-      operator T*() const { return reinterpret_cast<T*>(_ptr); }
-    };
-  }
 
   /// \brief Base interface for the GLPK LP and MIP solver
@@ -56,5 +44,6 @@
   protected:
 
-    _solver_bits::VoidPtr lp;
+    typedef glp_prob LPX;
+    glp_prob* lp;
 
     GlpkBase();
@@ -135,7 +124,7 @@
 
     ///Pointer to the underlying GLPK data structure.
-    _solver_bits::VoidPtr lpx() {return lp;}
+    LPX *lpx() {return lp;}
     ///Const pointer to the underlying GLPK data structure.
-    _solver_bits::VoidPtr lpx() const {return lp;}
+    const LPX *lpx() const {return lp;}
 
     ///Returns the constraint identifier understood by GLPK.
Index: lemon/hartmann_orlin.h
===================================================================
--- lemon/hartmann_orlin.h	(revision 891)
+++ lemon/hartmann_orlin.h	(revision 859)
@@ -107,9 +107,4 @@
   /// \tparam LEN The type of the length map. The default
   /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref HartmannOrlinDefaultTraits
-  /// "HartmannOrlinDefaultTraits<GR, LEN>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
 #ifdef DOXYGEN
   template <typename GR, typename LEN, typename TR>
@@ -133,5 +128,6 @@
     ///
     /// The large value type used for internal computations.
-    /// By default, it is \c long \c long if the \c Value type is integer,
+    /// Using the \ref HartmannOrlinDefaultTraits "default traits class",
+    /// it is \c long \c long if the \c Value type is integer,
     /// otherwise it is \c double.
     typedef typename TR::LargeValue LargeValue;
Index: lemon/howard.h
===================================================================
--- lemon/howard.h	(revision 891)
+++ lemon/howard.h	(revision 818)
@@ -107,9 +107,4 @@
   /// \tparam LEN The type of the length map. The default
   /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref HowardDefaultTraits
-  /// "HowardDefaultTraits<GR, LEN>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
 #ifdef DOXYGEN
   template <typename GR, typename LEN, typename TR>
@@ -133,5 +128,6 @@
     ///
     /// The large value type used for internal computations.
-    /// By default, it is \c long \c long if the \c Value type is integer,
+    /// Using the \ref HowardDefaultTraits "default traits class",
+    /// it is \c long \c long if the \c Value type is integer,
     /// otherwise it is \c double.
     typedef typename TR::LargeValue LargeValue;
Index: lemon/karp.h
===================================================================
--- lemon/karp.h	(revision 891)
+++ lemon/karp.h	(revision 819)
@@ -105,9 +105,4 @@
   /// \tparam LEN The type of the length map. The default
   /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref KarpDefaultTraits
-  /// "KarpDefaultTraits<GR, LEN>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
 #ifdef DOXYGEN
   template <typename GR, typename LEN, typename TR>
@@ -131,5 +126,6 @@
     ///
     /// The large value type used for internal computations.
-    /// By default, it is \c long \c long if the \c Value type is integer,
+    /// Using the \ref KarpDefaultTraits "default traits class",
+    /// it is \c long \c long if the \c Value type is integer,
     /// otherwise it is \c double.
     typedef typename TR::LargeValue LargeValue;
Index: lemon/lp_base.h
===================================================================
--- lemon/lp_base.h	(revision 833)
+++ lemon/lp_base.h	(revision 903)
@@ -1230,8 +1230,8 @@
       Row r;
       c.expr().simplify();
-      r._id = _addRowId(_addRow(c.lowerBounded()?c.lowerBound():-INF, 
+      r._id = _addRowId(_addRow(c.lowerBounded()?c.lowerBound()-*c.expr():-INF, 
                                 ExprIterator(c.expr().comps.begin(), cols),
                                 ExprIterator(c.expr().comps.end(), cols),
-                                c.upperBounded()?c.upperBound():INF));
+                                c.upperBounded()?c.upperBound()-*c.expr():INF));
       return r;
     }
Index: lemon/min_cost_arborescence.h
===================================================================
--- lemon/min_cost_arborescence.h	(revision 891)
+++ lemon/min_cost_arborescence.h	(revision 760)
@@ -113,9 +113,8 @@
   /// it is necessary. The default map type is \ref
   /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref MinCostArborescenceDefaultTraits
+  /// \param TR Traits class to set various data types used
+  /// by the algorithm. The default traits class is
+  /// \ref MinCostArborescenceDefaultTraits
   /// "MinCostArborescenceDefaultTraits<GR, CM>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
 #ifndef DOXYGEN
   template <typename GR,
@@ -124,5 +123,5 @@
               MinCostArborescenceDefaultTraits<GR, CM> >
 #else
-  template <typename GR, typename CM, typename TR>
+  template <typename GR, typename CM, typedef TR>
 #endif
   class MinCostArborescence {
Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h	(revision 898)
+++ lemon/network_simplex.h	(revision 835)
@@ -44,11 +44,11 @@
   /// \ref amo93networkflows, \ref dantzig63linearprog,
   /// \ref kellyoneill91netsimplex.
-  /// This algorithm is a highly efficient specialized version of the
-  /// linear programming simplex method directly for the minimum cost
-  /// flow problem.
+  /// This algorithm is a specialized version of the linear programming
+  /// simplex method directly for the minimum cost flow problem.
+  /// It is one of the most efficient solution methods.
   ///
-  /// In general, %NetworkSimplex is the fastest implementation available
-  /// in LEMON for this problem.
-  /// Moreover, it supports both directions of the supply/demand inequality
+  /// In general this class is the fastest implementation available
+  /// in LEMON for the minimum cost flow problem.
+  /// Moreover it supports both directions of the supply/demand inequality
   /// constraints. For more information, see \ref SupplyType.
   ///
@@ -59,10 +59,10 @@
   ///
   /// \tparam GR The digraph type the algorithm runs on.
-  /// \tparam V The number type used for flow amounts, capacity bounds
+  /// \tparam V The value type used for flow amounts, capacity bounds
   /// and supply values in the algorithm. By default, it is \c int.
-  /// \tparam C The number type used for costs and potentials in the
+  /// \tparam C The value type used for costs and potentials in the
   /// algorithm. By default, it is the same as \c V.
   ///
-  /// \warning Both number types must be signed and all input data must
+  /// \warning Both value types must be signed and all input data must
   /// be integer.
   ///
@@ -127,5 +127,5 @@
     /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
     /// proved to be the most efficient and the most robust on various
-    /// test inputs.
+    /// test inputs according to our benchmark tests.
     /// However, another pivot rule can be selected using the \ref run()
     /// function with the proper parameter.
@@ -165,5 +165,5 @@
 
     typedef std::vector<int> IntVector;
-    typedef std::vector<char> CharVector;
+    typedef std::vector<bool> BoolVector;
     typedef std::vector<Value> ValueVector;
     typedef std::vector<Cost> CostVector;
@@ -195,5 +195,4 @@
     IntVector _source;
     IntVector _target;
-    bool _arc_mixing;
 
     // Node and arc data
@@ -214,6 +213,6 @@
     IntVector _last_succ;
     IntVector _dirty_revs;
-    CharVector _forward;
-    CharVector _state;
+    BoolVector _forward;
+    IntVector _state;
     int _root;
 
@@ -223,6 +222,4 @@
     int stem, par_stem, new_stem;
     Value delta;
-    
-    const Value MAX;
 
   public:
@@ -246,5 +243,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -298,5 +295,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -337,5 +334,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -410,5 +407,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -513,5 +510,5 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const CharVector &_state;
+      const IntVector  &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -635,10 +632,9 @@
     NetworkSimplex(const GR& graph, bool arc_mixing = false) :
       _graph(graph), _node_id(graph), _arc_id(graph),
-      _arc_mixing(arc_mixing),
-      MAX(std::numeric_limits<Value>::max()),
       INF(std::numeric_limits<Value>::has_infinity ?
-          std::numeric_limits<Value>::infinity() : MAX)
+          std::numeric_limits<Value>::infinity() :
+          std::numeric_limits<Value>::max())
     {
-      // Check the number types
+      // Check the value types
       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
         "The flow type of NetworkSimplex must be signed");
@@ -646,5 +642,56 @@
         "The cost type of NetworkSimplex must be signed");
         
-      // Reset data structures
+      // Resize vectors
+      _node_num = countNodes(_graph);
+      _arc_num = countArcs(_graph);
+      int all_node_num = _node_num + 1;
+      int max_arc_num = _arc_num + 2 * _node_num;
+
+      _source.resize(max_arc_num);
+      _target.resize(max_arc_num);
+
+      _lower.resize(_arc_num);
+      _upper.resize(_arc_num);
+      _cap.resize(max_arc_num);
+      _cost.resize(max_arc_num);
+      _supply.resize(all_node_num);
+      _flow.resize(max_arc_num);
+      _pi.resize(all_node_num);
+
+      _parent.resize(all_node_num);
+      _pred.resize(all_node_num);
+      _forward.resize(all_node_num);
+      _thread.resize(all_node_num);
+      _rev_thread.resize(all_node_num);
+      _succ_num.resize(all_node_num);
+      _last_succ.resize(all_node_num);
+      _state.resize(max_arc_num);
+
+      // Copy the graph
+      int i = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+        _node_id[n] = i;
+      }
+      if (arc_mixing) {
+        // Store the arcs in a mixed order
+        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
+        int i = 0, j = 0;
+        for (ArcIt a(_graph); a != INVALID; ++a) {
+          _arc_id[a] = i;
+          _source[i] = _node_id[_graph.source(a)];
+          _target[i] = _node_id[_graph.target(a)];
+          if ((i += k) >= _arc_num) i = ++j;
+        }
+      } else {
+        // Store the arcs in the original order
+        int i = 0;
+        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
+          _arc_id[a] = i;
+          _source[i] = _node_id[_graph.source(a)];
+          _target[i] = _node_id[_graph.target(a)];
+        }
+      }
+      
+      // Reset parameters
       reset();
     }
@@ -681,5 +728,5 @@
     /// If it is not used before calling \ref run(), the upper bounds
     /// will be set to \ref INF on all arcs (i.e. the flow value will be
-    /// unbounded from above).
+    /// unbounded from above on each arc).
     ///
     /// \param map An arc map storing the upper bounds.
@@ -794,10 +841,10 @@
     /// \endcode
     ///
-    /// This function can be called more than once. All the given parameters
-    /// are kept for the next call, unless \ref resetParams() or \ref reset()
-    /// is used, thus only the modified parameters have to be set again.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class (or the last \ref reset() call), then the \ref reset()
-    /// function must be called.
+    /// This function can be called more than once. All the parameters
+    /// that have been given are kept for the next call, unless
+    /// \ref reset() is called, thus only the modified parameters
+    /// have to be set again. See \ref reset() for examples.
+    /// However, the underlying digraph must not be modified after this
+    /// class have been constructed, since it copies and extends the graph.
     ///
     /// \param pivot_rule The pivot rule that will be used during the
@@ -813,5 +860,4 @@
     ///
     /// \see ProblemType, PivotRule
-    /// \see resetParams(), reset()
     ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
       if (!init()) return INFEASIBLE;
@@ -825,10 +871,9 @@
     /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
     ///
-    /// It is useful for multiple \ref run() calls. Basically, all the given
-    /// parameters are kept for the next \ref run() call, unless
-    /// \ref resetParams() or \ref reset() is used.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class or the last \ref reset() call, then the \ref reset()
-    /// function must be used, otherwise \ref resetParams() is sufficient.
+    /// It is useful for multiple run() calls. If this function is not
+    /// used, all the parameters given before are kept for the next
+    /// \ref run() call.
+    /// However, the underlying digraph must not be modified after this
+    /// class have been constructed, since it copies and extends the graph.
     ///
     /// For example,
@@ -840,12 +885,12 @@
     ///     .supplyMap(sup).run();
     ///
-    ///   // Run again with modified cost map (resetParams() is not called,
+    ///   // Run again with modified cost map (reset() is not called,
     ///   // so only the cost map have to be set again)
     ///   cost[e] += 100;
     ///   ns.costMap(cost).run();
     ///
-    ///   // Run again from scratch using resetParams()
+    ///   // Run again from scratch using reset()
     ///   // (the lower bounds will be set to zero on all arcs)
-    ///   ns.resetParams();
+    ///   ns.reset();
     ///   ns.upperMap(capacity).costMap(cost)
     ///     .supplyMap(sup).run();
@@ -853,7 +898,5 @@
     ///
     /// \return <tt>(*this)</tt>
-    ///
-    /// \see reset(), run()
-    NetworkSimplex& resetParams() {
+    NetworkSimplex& reset() {
       for (int i = 0; i != _node_num; ++i) {
         _supply[i] = 0;
@@ -869,81 +912,4 @@
     }
 
-    /// \brief Reset the internal data structures and all the parameters
-    /// that have been given before.
-    ///
-    /// This function resets the internal data structures and all the
-    /// paramaters that have been given before using functions \ref lowerMap(),
-    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
-    /// \ref supplyType().
-    ///
-    /// It is useful for multiple \ref run() calls. Basically, all the given
-    /// parameters are kept for the next \ref run() call, unless
-    /// \ref resetParams() or \ref reset() is used.
-    /// If the underlying digraph was also modified after the construction
-    /// of the class or the last \ref reset() call, then the \ref reset()
-    /// function must be used, otherwise \ref resetParams() is sufficient.
-    ///
-    /// See \ref resetParams() for examples.
-    ///
-    /// \return <tt>(*this)</tt>
-    ///
-    /// \see resetParams(), run()
-    NetworkSimplex& reset() {
-      // Resize vectors
-      _node_num = countNodes(_graph);
-      _arc_num = countArcs(_graph);
-      int all_node_num = _node_num + 1;
-      int max_arc_num = _arc_num + 2 * _node_num;
-
-      _source.resize(max_arc_num);
-      _target.resize(max_arc_num);
-
-      _lower.resize(_arc_num);
-      _upper.resize(_arc_num);
-      _cap.resize(max_arc_num);
-      _cost.resize(max_arc_num);
-      _supply.resize(all_node_num);
-      _flow.resize(max_arc_num);
-      _pi.resize(all_node_num);
-
-      _parent.resize(all_node_num);
-      _pred.resize(all_node_num);
-      _forward.resize(all_node_num);
-      _thread.resize(all_node_num);
-      _rev_thread.resize(all_node_num);
-      _succ_num.resize(all_node_num);
-      _last_succ.resize(all_node_num);
-      _state.resize(max_arc_num);
-
-      // Copy the graph
-      int i = 0;
-      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
-        _node_id[n] = i;
-      }
-      if (_arc_mixing) {
-        // Store the arcs in a mixed order
-        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
-        int i = 0, j = 0;
-        for (ArcIt a(_graph); a != INVALID; ++a) {
-          _arc_id[a] = i;
-          _source[i] = _node_id[_graph.source(a)];
-          _target[i] = _node_id[_graph.target(a)];
-          if ((i += k) >= _arc_num) i = ++j;
-        }
-      } else {
-        // Store the arcs in the original order
-        int i = 0;
-        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
-          _arc_id[a] = i;
-          _source[i] = _node_id[_graph.source(a)];
-          _target[i] = _node_id[_graph.target(a)];
-        }
-      }
-      
-      // Reset parameters
-      resetParams();
-      return *this;
-    }
-    
     /// @}
 
@@ -1055,7 +1021,7 @@
           Value c = _lower[i];
           if (c >= 0) {
-            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
+            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
           } else {
-            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
+            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
           }
           _supply[_source[i]] -= c;
@@ -1249,5 +1215,5 @@
         e = _pred[u];
         d = _forward[u] ?
-          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
+          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
         if (d < delta) {
           delta = d;
@@ -1260,5 +1226,5 @@
         e = _pred[u];
         d = _forward[u] ? 
-          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
+          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
         if (d <= delta) {
           delta = d;
@@ -1459,5 +1425,5 @@
         findJoinNode();
         bool change = findLeavingArc();
-        if (delta >= MAX) return UNBOUNDED;
+        if (delta >= INF) return UNBOUNDED;
         changeFlow(change);
         if (change) {
Index: lemon/planarity.h
===================================================================
--- lemon/planarity.h	(revision 896)
+++ lemon/planarity.h	(revision 862)
@@ -519,7 +519,7 @@
   ///
   /// This function implements the Boyer-Myrvold algorithm for
-  /// planarity checking of an undirected simple graph. It is a simplified
+  /// planarity checking of an undirected graph. It is a simplified
   /// version of the PlanarEmbedding algorithm class because neither
-  /// the embedding nor the Kuratowski subdivisons are computed.
+  /// the embedding nor the kuratowski subdivisons are not computed.
   template <typename GR>
   bool checkPlanarity(const GR& graph) {
@@ -533,15 +533,14 @@
   ///
   /// This class implements the Boyer-Myrvold algorithm for planar
-  /// embedding of an undirected simple graph. The planar embedding is an
+  /// embedding of an undirected graph. The planar embedding is an
   /// ordering of the outgoing edges of the nodes, which is a possible
   /// configuration to draw the graph in the plane. If there is not
-  /// such ordering then the graph contains a K<sub>5</sub> (full graph
-  /// with 5 nodes) or a K<sub>3,3</sub> (complete bipartite graph on
-  /// 3 Red and 3 Blue nodes) subdivision.
+  /// such ordering then the graph contains a \f$ K_5 \f$ (full graph
+  /// with 5 nodes) or a \f$ K_{3,3} \f$ (complete bipartite graph on
+  /// 3 ANode and 3 BNode) subdivision.
   ///
   /// The current implementation calculates either an embedding or a
-  /// Kuratowski subdivision. The running time of the algorithm is O(n).
-  ///
-  /// \see PlanarDrawing, checkPlanarity()
+  /// Kuratowski subdivision. The running time of the algorithm is 
+  /// \f$ O(n) \f$.
   template <typename Graph>
   class PlanarEmbedding {
@@ -593,24 +592,20 @@
   public:
 
-    /// \brief The map type for storing the embedding
-    ///
-    /// The map type for storing the embedding.
-    /// \see embeddingMap()
+    /// \brief The map for store of embedding
     typedef typename Graph::template ArcMap<Arc> EmbeddingMap;
 
     /// \brief Constructor
     ///
-    /// Constructor.
-    /// \pre The graph must be simple, i.e. it should not
-    /// contain parallel or loop arcs.
+    /// \note The graph should be simple, i.e. parallel and loop arc
+    /// free.
     PlanarEmbedding(const Graph& graph)
       : _graph(graph), _embedding(_graph), _kuratowski(graph, false) {}
 
-    /// \brief Run the algorithm.
+    /// \brief Runs the algorithm.
     ///
-    /// This function runs the algorithm.
-    /// \param kuratowski If this parameter is set to \c false, then the
+    /// Runs the algorithm.
+    /// \param kuratowski If the parameter is false, then the
     /// algorithm does not compute a Kuratowski subdivision.
-    /// \return \c true if the graph is planar.
+    ///\return %True when the graph is planar.
     bool run(bool kuratowski = true) {
       typedef _planarity_bits::PlanarityVisitor<Graph> Visitor;
@@ -705,7 +700,7 @@
     }
 
-    /// \brief Give back the successor of an arc
+    /// \brief Gives back the successor of an arc
     ///
-    /// This function gives back the successor of an arc. It makes
+    /// Gives back the successor of an arc. This function makes
     /// possible to query the cyclic order of the outgoing arcs from
     /// a node.
@@ -714,21 +709,19 @@
     }
 
-    /// \brief Give back the calculated embedding map
+    /// \brief Gives back the calculated embedding map
     ///
-    /// This function gives back the calculated embedding map, which
-    /// contains the successor of each arc in the cyclic order of the
-    /// outgoing arcs of its source node.
+    /// The returned map contains the successor of each arc in the
+    /// graph.
     const EmbeddingMap& embeddingMap() const {
       return _embedding;
     }
 
-    /// \brief Give back \c true if the given edge is in the Kuratowski
+    /// \brief Gives back true if the undirected arc is in the
+    /// kuratowski subdivision
+    ///
+    /// Gives back true if the undirected arc is in the kuratowski
     /// subdivision
-    ///
-    /// This function gives back \c true if the given edge is in the found
-    /// Kuratowski subdivision.
-    /// \pre The \c run() function must be called with \c true parameter
-    /// before using this function.
-    bool kuratowski(const Edge& edge) const {
+    /// \note The \c run() had to be called with true value.
+    bool kuratowski(const Edge& edge) {
       return _kuratowski[edge];
     }
@@ -2067,13 +2060,11 @@
   ///
   /// The planar drawing algorithm calculates positions for the nodes
-  /// in the plane. These coordinates satisfy that if the edges are
-  /// represented with straight lines, then they will not intersect
+  /// in the plane which coordinates satisfy that if the arcs are
+  /// represented with straight lines then they will not intersect
   /// each other.
   ///
-  /// Scnyder's algorithm embeds the graph on an \c (n-2)x(n-2) size grid,
-  /// i.e. each node will be located in the \c [0..n-2]x[0..n-2] square.
+  /// Scnyder's algorithm embeds the graph on \c (n-2,n-2) size grid,
+  /// i.e. each node will be located in the \c [0,n-2]x[0,n-2] square.
   /// The time complexity of the algorithm is O(n).
-  ///
-  /// \see PlanarEmbedding
   template <typename Graph>
   class PlanarDrawing {
@@ -2082,7 +2073,7 @@
     TEMPLATE_GRAPH_TYPEDEFS(Graph);
 
-    /// \brief The point type for storing coordinates
+    /// \brief The point type for store coordinates
     typedef dim2::Point<int> Point;
-    /// \brief The map type for storing the coordinates of the nodes
+    /// \brief The map type for store coordinates
     typedef typename Graph::template NodeMap<Point> PointMap;
 
@@ -2091,6 +2082,5 @@
     ///
     /// Constructor
-    /// \pre The graph must be simple, i.e. it should not
-    /// contain parallel or loop arcs.
+    /// \pre The graph should be simple, i.e. loop and parallel arc free.
     PlanarDrawing(const Graph& graph)
       : _graph(graph), _point_map(graph) {}
@@ -2377,8 +2367,8 @@
   public:
 
-    /// \brief Calculate the node positions
+    /// \brief Calculates the node positions
     ///
-    /// This function calculates the node positions on the plane.
-    /// \return \c true if the graph is planar.
+    /// This function calculates the node positions.
+    /// \return %True if the graph is planar.
     bool run() {
       PlanarEmbedding<Graph> pe(_graph);
@@ -2389,11 +2379,10 @@
     }
 
-    /// \brief Calculate the node positions according to a
+    /// \brief Calculates the node positions according to a
     /// combinatorical embedding
     ///
-    /// This function calculates the node positions on the plane.
-    /// The given \c embedding map should contain a valid combinatorical
-    /// embedding, i.e. a valid cyclic order of the arcs.
-    /// It can be computed using PlanarEmbedding.
+    /// This function calculates the node locations. The \c embedding
+    /// parameter should contain a valid combinatorical embedding, i.e.
+    /// a valid cyclic order of the arcs.
     template <typename EmbeddingMap>
     void run(const EmbeddingMap& embedding) {
@@ -2435,13 +2424,12 @@
     /// \brief The coordinate of the given node
     ///
-    /// This function returns the coordinate of the given node.
+    /// The coordinate of the given node.
     Point operator[](const Node& node) const {
       return _point_map[node];
     }
 
-    /// \brief Return the grid embedding in a node map
+    /// \brief Returns the grid embedding in a \e NodeMap.
     ///
-    /// This function returns the grid embedding in a node map of
-    /// \c dim2::Point<int> coordinates.
+    /// Returns the grid embedding in a \e NodeMap of \c dim2::Point<int> .
     const PointMap& coords() const {
       return _point_map;
@@ -2483,20 +2471,21 @@
   ///
   /// The graph coloring problem is the coloring of the graph nodes
-  /// so that there are no adjacent nodes with the same color. The
-  /// planar graphs can always be colored with four colors, which is
-  /// proved by Appel and Haken. Their proofs provide a quadratic
+  /// that there are not adjacent nodes with the same color. The
+  /// planar graphs can be always colored with four colors, it is
+  /// proved by Appel and Haken and their proofs provide a quadratic
   /// time algorithm for four coloring, but it could not be used to
-  /// implement an efficient algorithm. The five and six coloring can be
-  /// made in linear time, but in this class, the five coloring has
+  /// implement efficient algorithm. The five and six coloring can be
+  /// made in linear time, but in this class the five coloring has
   /// quadratic worst case time complexity. The two coloring (if
   /// possible) is solvable with a graph search algorithm and it is
   /// implemented in \ref bipartitePartitions() function in LEMON. To
-  /// decide whether a planar graph is three colorable is NP-complete.
+  /// decide whether the planar graph is three colorable is
+  /// NP-complete.
   ///
   /// This class contains member functions for calculate colorings
   /// with five and six colors. The six coloring algorithm is a simple
   /// greedy coloring on the backward minimum outgoing order of nodes.
-  /// This order can be computed by selecting the node with least
-  /// outgoing arcs to unprocessed nodes in each phase. This order
+  /// This order can be computed as in each phase the node with least
+  /// outgoing arcs to unprocessed nodes is chosen. This order
   /// guarantees that when a node is chosen for coloring it has at
   /// most five already colored adjacents. The five coloring algorithm
@@ -2511,17 +2500,13 @@
     TEMPLATE_GRAPH_TYPEDEFS(Graph);
 
-    /// \brief The map type for storing color indices
+    /// \brief The map type for store color indexes
     typedef typename Graph::template NodeMap<int> IndexMap;
-    /// \brief The map type for storing colors
-    ///
-    /// The map type for storing colors.
-    /// \see Palette, Color
+    /// \brief The map type for store colors
     typedef ComposeMap<Palette, IndexMap> ColorMap;
 
     /// \brief Constructor
     ///
-    /// Constructor.
-    /// \pre The graph must be simple, i.e. it should not
-    /// contain parallel or loop arcs.
+    /// Constructor
+    /// \pre The graph should be simple, i.e. loop and parallel arc free.
     PlanarColoring(const Graph& graph)
       : _graph(graph), _color_map(graph), _palette(0) {
@@ -2534,32 +2519,32 @@
     }
 
-    /// \brief Return the node map of color indices
+    /// \brief Returns the \e NodeMap of color indexes
     ///
-    /// This function returns the node map of color indices. The values are
-    /// in the range \c [0..4] or \c [0..5] according to the coloring method.
+    /// Returns the \e NodeMap of color indexes. The values are in the
+    /// range \c [0..4] or \c [0..5] according to the coloring method.
     IndexMap colorIndexMap() const {
       return _color_map;
     }
 
-    /// \brief Return the node map of colors
+    /// \brief Returns the \e NodeMap of colors
     ///
-    /// This function returns the node map of colors. The values are among
-    /// five or six distinct \ref lemon::Color "colors".
+    /// Returns the \e NodeMap of colors. The values are five or six
+    /// distinct \ref lemon::Color "colors".
     ColorMap colorMap() const {
       return composeMap(_palette, _color_map);
     }
 
-    /// \brief Return the color index of the node
+    /// \brief Returns the color index of the node
     ///
-    /// This function returns the color index of the given node. The value is
-    /// in the range \c [0..4] or \c [0..5] according to the coloring method.
+    /// Returns the color index of the node. The values are in the
+    /// range \c [0..4] or \c [0..5] according to the coloring method.
     int colorIndex(const Node& node) const {
       return _color_map[node];
     }
 
-    /// \brief Return the color of the node
+    /// \brief Returns the color of the node
     ///
-    /// This function returns the color of the given node. The value is among
-    /// five or six distinct \ref lemon::Color "colors".
+    /// Returns the color of the node. The values are five or six
+    /// distinct \ref lemon::Color "colors".
     Color color(const Node& node) const {
       return _palette[_color_map[node]];
@@ -2567,12 +2552,12 @@
 
 
-    /// \brief Calculate a coloring with at most six colors
+    /// \brief Calculates a coloring with at most six colors
     ///
     /// This function calculates a coloring with at most six colors. The time
     /// complexity of this variant is linear in the size of the graph.
-    /// \return \c true if the algorithm could color the graph with six colors.
-    /// If the algorithm fails, then the graph is not planar.
-    /// \note This function can return \c true if the graph is not
-    /// planar, but it can be colored with at most six colors.
+    /// \return %True when the algorithm could color the graph with six color.
+    /// If the algorithm fails, then the graph could not be planar.
+    /// \note This function can return true if the graph is not
+    /// planar but it can be colored with 6 colors.
     bool runSixColoring() {
 
@@ -2676,12 +2661,9 @@
   public:
 
-    /// \brief Calculate a coloring with at most five colors
+    /// \brief Calculates a coloring with at most five colors
     ///
     /// This function calculates a coloring with at most five
     /// colors. The worst case time complexity of this variant is
     /// quadratic in the size of the graph.
-    /// \param embedding This map should contain a valid combinatorical
-    /// embedding, i.e. a valid cyclic order of the arcs.
-    /// It can be computed using PlanarEmbedding.
     template <typename EmbeddingMap>
     void runFiveColoring(const EmbeddingMap& embedding) {
@@ -2730,10 +2712,10 @@
     }
 
-    /// \brief Calculate a coloring with at most five colors
+    /// \brief Calculates a coloring with at most five colors
     ///
     /// This function calculates a coloring with at most five
     /// colors. The worst case time complexity of this variant is
     /// quadratic in the size of the graph.
-    /// \return \c true if the graph is planar.
+    /// \return %True when the graph is planar.
     bool runFiveColoring() {
       PlanarEmbedding<Graph> pe(_graph);
Index: lemon/preflow.h
===================================================================
--- lemon/preflow.h	(revision 891)
+++ lemon/preflow.h	(revision 835)
@@ -114,15 +114,7 @@
   /// second phase constructs a feasible maximum flow on each arc.
   ///
-  /// \warning This implementation cannot handle infinite or very large
-  /// capacities (e.g. the maximum value of \c CAP::Value).
-  ///
   /// \tparam GR The type of the digraph the algorithm runs on.
   /// \tparam CAP The type of the capacity map. The default map
   /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
-  /// \tparam TR The traits class that defines various types used by the
-  /// algorithm. By default, it is \ref PreflowDefaultTraits
-  /// "PreflowDefaultTraits<GR, CAP>".
-  /// In most cases, this parameter should not be set directly,
-  /// consider to use the named template parameters instead.
 #ifdef DOXYGEN
   template <typename GR, typename CAP, typename TR>
Index: lemon/unionfind.h
===================================================================
--- lemon/unionfind.h	(revision 864)
+++ lemon/unionfind.h	(revision 833)
@@ -740,5 +740,5 @@
     void clear() {
       items.clear();
-      classes.clear();
+      classes.clear;
       firstClass = firstFreeClass = firstFreeItem = -1;
     }
Index: test/min_cost_flow_test.cc
===================================================================
--- test/min_cost_flow_test.cc	(revision 898)
+++ test/min_cost_flow_test.cc	(revision 716)
@@ -25,10 +25,6 @@
 
 #include <lemon/network_simplex.h>
-#include <lemon/capacity_scaling.h>
-#include <lemon/cost_scaling.h>
-#include <lemon/cycle_canceling.h>
 
 #include <lemon/concepts/digraph.h>
-#include <lemon/concepts/heap.h>
 #include <lemon/concept_check.h>
 
@@ -37,5 +33,4 @@
 using namespace lemon;
 
-// Test networks
 char test_lgf[] =
   "@nodes\n"
@@ -82,56 +77,4 @@
   "target 12\n";
 
-char test_neg1_lgf[] =
-  "@nodes\n"
-  "label   sup\n"
-  "    1   100\n"
-  "    2     0\n"
-  "    3     0\n"
-  "    4  -100\n"
-  "    5     0\n"
-  "    6     0\n"
-  "    7     0\n"
-  "@arcs\n"
-  "      cost   low1   low2\n"
-  "1 2    100      0      0\n"
-  "1 3     30      0      0\n"
-  "2 4     20      0      0\n"
-  "3 4     80      0      0\n"
-  "3 2     50      0      0\n"
-  "5 3     10      0      0\n"
-  "5 6     80      0   1000\n"
-  "6 7     30      0  -1000\n"
-  "7 5   -120      0      0\n";
-  
-char test_neg2_lgf[] =
-  "@nodes\n"
-  "label   sup\n"
-  "    1   100\n"
-  "    2  -300\n"
-  "@arcs\n"
-  "      cost\n"
-  "1 2     -1\n";
-
-
-// Test data
-typedef ListDigraph Digraph;
-DIGRAPH_TYPEDEFS(ListDigraph);
-
-Digraph gr;
-Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
-Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
-ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
-Node v, w;
-
-Digraph neg1_gr;
-Digraph::ArcMap<int> neg1_c(neg1_gr), neg1_l1(neg1_gr), neg1_l2(neg1_gr);
-ConstMap<Arc, int> neg1_u1(std::numeric_limits<int>::max()), neg1_u2(5000);
-Digraph::NodeMap<int> neg1_s(neg1_gr);
-
-Digraph neg2_gr;
-Digraph::ArcMap<int> neg2_c(neg2_gr);
-ConstMap<Arc, int> neg2_l(0), neg2_u(1000);
-Digraph::NodeMap<int> neg2_s(neg2_gr);
-
 
 enum SupplyType {
@@ -140,5 +83,4 @@
   LEQ
 };
-
 
 // Check the interface of an MCF algorithm
@@ -158,5 +100,5 @@
       const MCF& const_mcf = mcf;
 
-      b = mcf.reset().resetParams()
+      b = mcf.reset()
              .lowerMap(me.lower)
              .upperMap(me.upper)
@@ -327,97 +269,28 @@
 }
 
-template < typename MCF, typename Param >
-void runMcfGeqTests( Param param,
-                     const std::string &test_str = "",
-                     bool full_neg_cost_support = false )
-{
-  MCF mcf1(gr), mcf2(neg1_gr), mcf3(neg2_gr);
-  
-  // Basic tests
-  mcf1.upperMap(u).costMap(c).supplyMap(s1);
-  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s1,
-           mcf1.OPTIMAL, true,     5240, test_str + "-1");
-  mcf1.stSupply(v, w, 27);
-  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s2,
-           mcf1.OPTIMAL, true,     7620, test_str + "-2");
-  mcf1.lowerMap(l2).supplyMap(s1);
-  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s1,
-           mcf1.OPTIMAL, true,     5970, test_str + "-3");
-  mcf1.stSupply(v, w, 27);
-  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s2,
-           mcf1.OPTIMAL, true,     8010, test_str + "-4");
-  mcf1.resetParams().supplyMap(s1);
-  checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s1,
-           mcf1.OPTIMAL, true,       74, test_str + "-5");
-  mcf1.lowerMap(l2).stSupply(v, w, 27);
-  checkMcf(mcf1, mcf1.run(param), gr, l2, cu, cc, s2,
-           mcf1.OPTIMAL, true,       94, test_str + "-6");
-  mcf1.reset();
-  checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s3,
-           mcf1.OPTIMAL, true,        0, test_str + "-7");
-  mcf1.lowerMap(l2).upperMap(u);
-  checkMcf(mcf1, mcf1.run(param), gr, l2, u, cc, s3,
-           mcf1.INFEASIBLE, false,    0, test_str + "-8");
-  mcf1.lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
-  checkMcf(mcf1, mcf1.run(param), gr, l3, u, c, s4,
-           mcf1.OPTIMAL, true,     6360, test_str + "-9");
-
-  // Tests for the GEQ form
-  mcf1.resetParams().upperMap(u).costMap(c).supplyMap(s5);
-  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s5,
-           mcf1.OPTIMAL, true,     3530, test_str + "-10", GEQ);
-  mcf1.lowerMap(l2);
-  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
-           mcf1.OPTIMAL, true,     4540, test_str + "-11", GEQ);
-  mcf1.supplyMap(s6);
-  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
-           mcf1.INFEASIBLE, false,    0, test_str + "-12", GEQ);
-
-  // Tests with negative costs
-  mcf2.lowerMap(neg1_l1).costMap(neg1_c).supplyMap(neg1_s);
-  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u1, neg1_c, neg1_s,
-           mcf2.UNBOUNDED, false,     0, test_str + "-13");
-  mcf2.upperMap(neg1_u2);
-  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u2, neg1_c, neg1_s,
-           mcf2.OPTIMAL, true,   -40000, test_str + "-14");
-  mcf2.resetParams().lowerMap(neg1_l2).costMap(neg1_c).supplyMap(neg1_s);
-  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l2, neg1_u1, neg1_c, neg1_s,
-           mcf2.UNBOUNDED, false,     0, test_str + "-15");
-
-  mcf3.costMap(neg2_c).supplyMap(neg2_s);
-  if (full_neg_cost_support) {
-    checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
-             mcf3.OPTIMAL, true,   -300, test_str + "-16", GEQ);
-  } else {
-    checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
-             mcf3.UNBOUNDED, false,   0, test_str + "-17", GEQ);
-  }
-  mcf3.upperMap(neg2_u);
-  checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
-           mcf3.OPTIMAL, true,     -300, test_str + "-18", GEQ);
-}
-
-template < typename MCF, typename Param >
-void runMcfLeqTests( Param param,
-                     const std::string &test_str = "" )
-{
-  // Tests for the LEQ form
-  MCF mcf1(gr);
-  mcf1.supplyType(mcf1.LEQ);
-  mcf1.upperMap(u).costMap(c).supplyMap(s6);
-  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s6,
-           mcf1.OPTIMAL, true,   5080, test_str + "-19", LEQ);
-  mcf1.lowerMap(l2);
-  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
-           mcf1.OPTIMAL, true,   5930, test_str + "-20", LEQ);
-  mcf1.supplyMap(s5);
-  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
-           mcf1.INFEASIBLE, false,  0, test_str + "-21", LEQ);
-}
-
-
 int main()
 {
-  // Read the test networks
+  // Check the interfaces
+  {
+    typedef concepts::Digraph GR;
+    checkConcept< McfClassConcept<GR, int, int>,
+                  NetworkSimplex<GR> >();
+    checkConcept< McfClassConcept<GR, double, double>,
+                  NetworkSimplex<GR, double> >();
+    checkConcept< McfClassConcept<GR, int, double>,
+                  NetworkSimplex<GR, int, double> >();
+  }
+
+  // Run various MCF tests
+  typedef ListDigraph Digraph;
+  DIGRAPH_TYPEDEFS(ListDigraph);
+
+  // Read the test digraph
+  Digraph gr;
+  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
+  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
+  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
+  Node v, w;
+
   std::istringstream input(test_lgf);
   DigraphReader<Digraph>(gr, input)
@@ -437,104 +310,139 @@
     .run();
   
-  std::istringstream neg_inp1(test_neg1_lgf);
-  DigraphReader<Digraph>(neg1_gr, neg_inp1)
-    .arcMap("cost", neg1_c)
-    .arcMap("low1", neg1_l1)
-    .arcMap("low2", neg1_l2)
-    .nodeMap("sup", neg1_s)
-    .run();
-  
-  std::istringstream neg_inp2(test_neg2_lgf);
-  DigraphReader<Digraph>(neg2_gr, neg_inp2)
-    .arcMap("cost", neg2_c)
-    .nodeMap("sup", neg2_s)
-    .run();
-  
-  // Check the interface of NetworkSimplex
+  // Build test digraphs with negative costs
+  Digraph neg_gr;
+  Node n1 = neg_gr.addNode();
+  Node n2 = neg_gr.addNode();
+  Node n3 = neg_gr.addNode();
+  Node n4 = neg_gr.addNode();
+  Node n5 = neg_gr.addNode();
+  Node n6 = neg_gr.addNode();
+  Node n7 = neg_gr.addNode();
+  
+  Arc a1 = neg_gr.addArc(n1, n2);
+  Arc a2 = neg_gr.addArc(n1, n3);
+  Arc a3 = neg_gr.addArc(n2, n4);
+  Arc a4 = neg_gr.addArc(n3, n4);
+  Arc a5 = neg_gr.addArc(n3, n2);
+  Arc a6 = neg_gr.addArc(n5, n3);
+  Arc a7 = neg_gr.addArc(n5, n6);
+  Arc a8 = neg_gr.addArc(n6, n7);
+  Arc a9 = neg_gr.addArc(n7, n5);
+  
+  Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
+  ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
+  Digraph::NodeMap<int> neg_s(neg_gr, 0);
+  
+  neg_l2[a7] =  1000;
+  neg_l2[a8] = -1000;
+  
+  neg_s[n1] =  100;
+  neg_s[n4] = -100;
+  
+  neg_c[a1] =  100;
+  neg_c[a2] =   30;
+  neg_c[a3] =   20;
+  neg_c[a4] =   80;
+  neg_c[a5] =   50;
+  neg_c[a6] =   10;
+  neg_c[a7] =   80;
+  neg_c[a8] =   30;
+  neg_c[a9] = -120;
+
+  Digraph negs_gr;
+  Digraph::NodeMap<int> negs_s(negs_gr);
+  Digraph::ArcMap<int> negs_c(negs_gr);
+  ConstMap<Arc, int> negs_l(0), negs_u(1000);
+  n1 = negs_gr.addNode();
+  n2 = negs_gr.addNode();
+  negs_s[n1] = 100;
+  negs_s[n2] = -300;
+  negs_c[negs_gr.addArc(n1, n2)] = -1;
+
+
+  // A. Test NetworkSimplex with the default pivot rule
   {
-    typedef concepts::Digraph GR;
-    checkConcept< McfClassConcept<GR, int, int>,
-                  NetworkSimplex<GR> >();
-    checkConcept< McfClassConcept<GR, double, double>,
-                  NetworkSimplex<GR, double> >();
-    checkConcept< McfClassConcept<GR, int, double>,
-                  NetworkSimplex<GR, int, double> >();
-  }
-
-  // Check the interface of CapacityScaling
+    NetworkSimplex<Digraph> mcf(gr);
+
+    // Check the equality form
+    mcf.upperMap(u).costMap(c);
+    checkMcf(mcf, mcf.supplyMap(s1).run(),
+             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
+    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
+             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
+    mcf.lowerMap(l2);
+    checkMcf(mcf, mcf.supplyMap(s1).run(),
+             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
+    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
+             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
+    mcf.reset();
+    checkMcf(mcf, mcf.supplyMap(s1).run(),
+             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
+    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
+             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
+    mcf.reset();
+    checkMcf(mcf, mcf.run(),
+             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
+    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
+             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
+    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
+    checkMcf(mcf, mcf.run(),
+             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
+
+    // Check the GEQ form
+    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
+    checkMcf(mcf, mcf.run(),
+             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
+    mcf.supplyType(mcf.GEQ);
+    checkMcf(mcf, mcf.lowerMap(l2).run(),
+             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
+    mcf.supplyMap(s6);
+    checkMcf(mcf, mcf.run(),
+             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
+
+    // Check the LEQ form
+    mcf.reset().supplyType(mcf.LEQ);
+    mcf.upperMap(u).costMap(c).supplyMap(s6);
+    checkMcf(mcf, mcf.run(),
+             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
+    checkMcf(mcf, mcf.lowerMap(l2).run(),
+             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
+    mcf.supplyMap(s5);
+    checkMcf(mcf, mcf.run(),
+             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
+
+    // Check negative costs
+    NetworkSimplex<Digraph> neg_mcf(neg_gr);
+    neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
+    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
+      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A16");
+    neg_mcf.upperMap(neg_u2);
+    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
+      neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");
+    neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
+    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,
+      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A18");
+      
+    NetworkSimplex<Digraph> negs_mcf(negs_gr);
+    negs_mcf.costMap(negs_c).supplyMap(negs_s);
+    checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
+      negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);
+  }
+
+  // B. Test NetworkSimplex with each pivot rule
   {
-    typedef concepts::Digraph GR;
-    checkConcept< McfClassConcept<GR, int, int>,
-                  CapacityScaling<GR> >();
-    checkConcept< McfClassConcept<GR, double, double>,
-                  CapacityScaling<GR, double> >();
-    checkConcept< McfClassConcept<GR, int, double>,
-                  CapacityScaling<GR, int, double> >();
-    typedef CapacityScaling<GR>::
-      SetHeap<concepts::Heap<int, RangeMap<int> > >::Create CAS;
-    checkConcept< McfClassConcept<GR, int, int>, CAS >();
-  }
-
-  // Check the interface of CostScaling
-  {
-    typedef concepts::Digraph GR;
-    checkConcept< McfClassConcept<GR, int, int>,
-                  CostScaling<GR> >();
-    checkConcept< McfClassConcept<GR, double, double>,
-                  CostScaling<GR, double> >();
-    checkConcept< McfClassConcept<GR, int, double>,
-                  CostScaling<GR, int, double> >();
-    typedef CostScaling<GR>::
-      SetLargeCost<double>::Create COS;
-    checkConcept< McfClassConcept<GR, int, int>, COS >();
-  }
-
-  // Check the interface of CycleCanceling
-  {
-    typedef concepts::Digraph GR;
-    checkConcept< McfClassConcept<GR, int, int>,
-                  CycleCanceling<GR> >();
-    checkConcept< McfClassConcept<GR, double, double>,
-                  CycleCanceling<GR, double> >();
-    checkConcept< McfClassConcept<GR, int, double>,
-                  CycleCanceling<GR, int, double> >();
-  }
-
-  // Test NetworkSimplex
-  { 
-    typedef NetworkSimplex<Digraph> MCF;
-    runMcfGeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE", true);
-    runMcfLeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE");
-    runMcfGeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE", true);
-    runMcfLeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE");
-    runMcfGeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS", true);
-    runMcfLeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS");
-    runMcfGeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL", true);
-    runMcfLeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL");
-    runMcfGeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL", true);
-    runMcfLeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL");
-  }
-  
-  // Test CapacityScaling
-  {
-    typedef CapacityScaling<Digraph> MCF;
-    runMcfGeqTests<MCF>(0, "SSP");
-    runMcfGeqTests<MCF>(2, "CAS");
-  }
-
-  // Test CostScaling
-  {
-    typedef CostScaling<Digraph> MCF;
-    runMcfGeqTests<MCF>(MCF::PUSH, "COS-PR");
-    runMcfGeqTests<MCF>(MCF::AUGMENT, "COS-AR");
-    runMcfGeqTests<MCF>(MCF::PARTIAL_AUGMENT, "COS-PAR");
-  }
-
-  // Test CycleCanceling
-  {
-    typedef CycleCanceling<Digraph> MCF;
-    runMcfGeqTests<MCF>(MCF::SIMPLE_CYCLE_CANCELING, "SCC");
-    runMcfGeqTests<MCF>(MCF::MINIMUM_MEAN_CYCLE_CANCELING, "MMCC");
-    runMcfGeqTests<MCF>(MCF::CANCEL_AND_TIGHTEN, "CAT");
+    NetworkSimplex<Digraph> mcf(gr);
+    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
+
+    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
+             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
+    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
+             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
+    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
+             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
+    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
+             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
+    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
+             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
   }
 
