# HG changeset patch
# User deba
# Date 1178538138 0
# Node ID c9218405595bbea9497f21ca3d10e0e75a2469ff
# Parent  3f1c7a6c33cd294b2191f0654da2ade1f841691b
Various min cost flow solvers

Patch from Peter Kovacs

diff -r 3f1c7a6c33cd -r c9218405595b lemon/Makefile.am
--- a/lemon/Makefile.am	Mon May 07 08:49:57 2007 +0000
+++ b/lemon/Makefile.am	Mon May 07 11:42:18 2007 +0000
@@ -42,12 +42,14 @@
 	lemon/bp_matching.h \
 	lemon/bpugraph_adaptor.h \
 	lemon/bucket_heap.h \
+	lemon/capacity_scaling.h \
 	lemon/circulation.h \
 	lemon/color.h \
 	lemon/config.h \
 	lemon/concept_check.h \
 	lemon/counter.h \
 	lemon/csp.h \
+	lemon/cycle_canceling.h \
 	lemon/dag_shortest_path.h \
 	lemon/dfs.h \
 	lemon/dijkstra.h \
@@ -89,10 +91,13 @@
 	lemon/matrix_maps.h \
 	lemon/max_matching.h \
 	lemon/min_cost_arborescence.h \
+	lemon/min_cost_flow.h \
+	lemon/min_cost_max_flow.h \
 	lemon/min_mean_cycle.h \
 	lemon/mip_glpk.h \
 	lemon/mip_cplex.h \
 	lemon/nagamochi_ibaraki.h \
+	lemon/network_simplex.h \
 	lemon/path.h \
 	lemon/path_utils.h \
 	lemon/polynomial.h \
diff -r 3f1c7a6c33cd -r c9218405595b lemon/capacity_scaling.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/capacity_scaling.h	Mon May 07 11:42:18 2007 +0000
@@ -0,0 +1,680 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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
+///
+/// \file
+/// \brief The capacity scaling algorithm for finding a minimum cost
+/// flow.
+
+#include <vector>
+#include <lemon/dijkstra.h>
+#include <lemon/graph_adaptor.h>
+
+#define WITH_SCALING
+
+namespace lemon {
+
+  /// \addtogroup min_cost_flow
+  /// @{
+
+  /// \brief Implementation of the capacity scaling version of the
+  /// succesive shortest path algorithm for finding a minimum cost flow.
+  ///
+  /// \ref lemon::CapacityScaling "CapacityScaling" implements a
+  /// capacity scaling algorithm for finding a minimum cost flow.
+  ///
+  /// \param Graph The directed graph type the algorithm runs on.
+  /// \param LowerMap The type of the lower bound map.
+  /// \param CapacityMap The type of the capacity (upper bound) map.
+  /// \param CostMap The type of the cost (length) map.
+  /// \param SupplyMap The type of the supply map.
+  ///
+  /// \warning
+  /// - Edge capacities and costs should be nonnegative integers.
+  ///	However \c CostMap::Value should be signed type.
+  /// - Supply values should be integers.
+  /// - \c LowerMap::Value must be convertible to
+  ///	\c CapacityMap::Value and \c CapacityMap::Value must be
+  ///	convertible to \c SupplyMap::Value.
+  ///
+  /// \author Peter Kovacs
+
+template < typename Graph,
+	   typename LowerMap = typename Graph::template EdgeMap<int>,
+	   typename CapacityMap = LowerMap,
+	   typename CostMap = typename Graph::template EdgeMap<int>,
+	   typename SupplyMap = typename Graph::template NodeMap
+				<typename CapacityMap::Value> >
+  class CapacityScaling
+  {
+    typedef typename Graph::Node Node;
+    typedef typename Graph::NodeIt NodeIt;
+    typedef typename Graph::Edge Edge;
+    typedef typename Graph::EdgeIt EdgeIt;
+    typedef typename Graph::InEdgeIt InEdgeIt;
+    typedef typename Graph::OutEdgeIt OutEdgeIt;
+
+    typedef typename LowerMap::Value Lower;
+    typedef typename CapacityMap::Value Capacity;
+    typedef typename CostMap::Value Cost;
+    typedef typename SupplyMap::Value Supply;
+    typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
+    typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
+
+    typedef ResGraphAdaptor< const Graph, Capacity,
+			     CapacityRefMap, CapacityRefMap > ResGraph;
+    typedef typename ResGraph::Node ResNode;
+    typedef typename ResGraph::NodeIt ResNodeIt;
+    typedef typename ResGraph::Edge ResEdge;
+    typedef typename ResGraph::EdgeIt ResEdgeIt;
+
+  public:
+
+    /// \brief The type of the flow map.
+    typedef CapacityRefMap FlowMap;
+    /// \brief The type of the potential map.
+    typedef typename Graph::template NodeMap<Cost> PotentialMap;
+
+  protected:
+
+    /// \brief Map adaptor class for handling reduced edge costs.
+    class ReducedCostMap : public MapBase<ResEdge, Cost>
+    {
+    private:
+
+      const ResGraph &gr;
+      const CostMap &cost_map;
+      const PotentialMap &pot_map;
+
+    public:
+
+      typedef typename MapBase<ResEdge, Cost>::Value Value;
+      typedef typename MapBase<ResEdge, Cost>::Key Key;
+
+      ReducedCostMap( const ResGraph &_gr,
+		      const CostMap &_cost,
+		      const PotentialMap &_pot ) :
+	gr(_gr), cost_map(_cost), pot_map(_pot) {}
+
+      Value operator[](const Key &e) const {
+	return ResGraph::forward(e) ?
+	   cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)] :
+	  -cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
+      }
+
+    }; //class ReducedCostMap
+
+    /// \brief Map adaptor for \ref lemon::Dijkstra "Dijkstra" class to
+    /// update node potentials.
+    class PotentialUpdateMap : public MapBase<ResNode, Cost>
+    {
+    private:
+
+      PotentialMap *pot;
+      typedef std::pair<ResNode, Cost> Pair;
+      std::vector<Pair> data;
+
+    public:
+
+      typedef typename MapBase<ResNode, Cost>::Value Value;
+      typedef typename MapBase<ResNode, Cost>::Key Key;
+
+      void potentialMap(PotentialMap &_pot) {
+	pot = &_pot;
+      }
+
+      void init() {
+	data.clear();
+      }
+
+      void set(const Key &n, const Value &v) {
+	data.push_back(Pair(n, v));
+      }
+
+      void update() {
+	Cost val = data[data.size()-1].second;
+	for (int i = 0; i < data.size()-1; ++i)
+	  (*pot)[data[i].first] += val - data[i].second;
+      }
+
+    }; //class PotentialUpdateMap
+
+#ifdef WITH_SCALING
+    /// \brief Map adaptor class for identifing deficit nodes.
+    class DeficitBoolMap : public MapBase<ResNode, bool>
+    {
+    private:
+
+      const SupplyRefMap &imb;
+      const Capacity &delta;
+
+    public:
+
+      DeficitBoolMap(const SupplyRefMap &_imb, const Capacity &_delta) :
+	imb(_imb), delta(_delta) {}
+
+      bool operator[](const ResNode &n) const {
+	return imb[n] <= -delta;
+      }
+
+    }; //class DeficitBoolMap
+
+    /// \brief Map adaptor class for filtering edges with at least
+    /// \c delta residual capacity
+    class DeltaFilterMap : public MapBase<ResEdge, bool>
+    {
+    private:
+
+      const ResGraph &gr;
+      const Capacity &delta;
+
+    public:
+
+      typedef typename MapBase<ResEdge, Cost>::Value Value;
+      typedef typename MapBase<ResEdge, Cost>::Key Key;
+
+      DeltaFilterMap(const ResGraph &_gr, const Capacity &_delta) :
+	gr(_gr), delta(_delta) {}
+
+      Value operator[](const Key &e) const {
+	return gr.rescap(e) >= delta;
+      }
+
+    }; //class DeltaFilterMap
+
+    typedef EdgeSubGraphAdaptor<const ResGraph, const DeltaFilterMap>
+      DeltaResGraph;
+
+    /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
+    class ResDijkstraTraits :
+      public DijkstraDefaultTraits<DeltaResGraph, ReducedCostMap>
+    {
+    public:
+
+      typedef PotentialUpdateMap DistMap;
+
+      static DistMap *createDistMap(const DeltaResGraph&) {
+	return new DistMap();
+      }
+
+    }; //class ResDijkstraTraits
+
+#else //WITHOUT_CAPACITY_SCALING
+    /// \brief Map adaptor class for identifing deficit nodes.
+    class DeficitBoolMap : public MapBase<ResNode, bool>
+    {
+    private:
+
+      const SupplyRefMap &imb;
+
+    public:
+
+      DeficitBoolMap(const SupplyRefMap &_imb) : imb(_imb) {}
+
+      bool operator[](const ResNode &n) const {
+	return imb[n] < 0;
+      }
+
+    }; //class DeficitBoolMap
+
+    /// \brief Traits class for \ref lemon::Dijkstra "Dijkstra" class.
+    class ResDijkstraTraits :
+      public DijkstraDefaultTraits<ResGraph, ReducedCostMap>
+    {
+    public:
+
+      typedef PotentialUpdateMap DistMap;
+
+      static DistMap *createDistMap(const ResGraph&) {
+	return new DistMap();
+      }
+
+    }; //class ResDijkstraTraits
+#endif
+
+  protected:
+
+    /// \brief The directed graph the algorithm runs on.
+    const Graph &graph;
+    /// \brief The original lower bound map.
+    const LowerMap *lower;
+    /// \brief The modified capacity map.
+    CapacityRefMap capacity;
+    /// \brief The cost map.
+    const CostMap &cost;
+    /// \brief The modified supply map.
+    SupplyRefMap supply;
+    /// \brief The sum of supply values equals zero.
+    bool valid_supply;
+
+    /// \brief The edge map of the current flow.
+    FlowMap flow;
+    /// \brief The potential node map.
+    PotentialMap potential;
+    /// \brief The residual graph.
+    ResGraph res_graph;
+    /// \brief The reduced cost map.
+    ReducedCostMap red_cost;
+
+    /// \brief The imbalance map.
+    SupplyRefMap imbalance;
+    /// \brief The excess nodes.
+    std::vector<ResNode> excess_nodes;
+    /// \brief The index of the next excess node.
+    int next_node;
+
+#ifdef WITH_SCALING
+    typedef Dijkstra<DeltaResGraph, ReducedCostMap, ResDijkstraTraits>
+      ResDijkstra;
+    /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
+    /// shortest paths in the residual graph with respect to the
+    /// reduced edge costs.
+    ResDijkstra dijkstra;
+
+    /// \brief The delta parameter used for capacity scaling.
+    Capacity delta;
+    /// \brief Edge filter map.
+    DeltaFilterMap delta_filter;
+    /// \brief The delta residual graph.
+    DeltaResGraph dres_graph;
+    /// \brief Map for identifing deficit nodes.
+    DeficitBoolMap delta_deficit;
+
+#else //WITHOUT_CAPACITY_SCALING
+    typedef Dijkstra<ResGraph, ReducedCostMap, ResDijkstraTraits>
+      ResDijkstra;
+    /// \brief \ref lemon::Dijkstra "Dijkstra" class for finding
+    /// shortest paths in the residual graph with respect to the
+    /// reduced edge costs.
+    ResDijkstra dijkstra;
+    /// \brief Map for identifing deficit nodes.
+    DeficitBoolMap has_deficit;
+#endif
+
+    /// \brief Pred map for the \ref lemon::Dijkstra "Dijkstra" class.
+    typename ResDijkstra::PredMap pred;
+    /// \brief Dist map for the \ref lemon::Dijkstra "Dijkstra" class to
+    /// update node potentials.
+    PotentialUpdateMap updater;
+
+  public :
+
+    /// \brief General constructor of the class (with lower bounds).
+    ///
+    /// General constructor of the class (with lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _lower The lower bounds of the edges.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _supply The supply values of the nodes (signed).
+    CapacityScaling( const Graph &_graph,
+		     const LowerMap &_lower,
+		     const CapacityMap &_capacity,
+		     const CostMap &_cost,
+		     const SupplyMap &_supply ) :
+      graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
+      supply(_graph), flow(_graph, 0), potential(_graph, 0),
+      res_graph(_graph, capacity, flow),
+      red_cost(res_graph, cost, potential), imbalance(_graph),
+#ifdef WITH_SCALING
+      delta(0), delta_filter(res_graph, delta),
+      dres_graph(res_graph, delta_filter),
+      dijkstra(dres_graph, red_cost), pred(dres_graph),
+      delta_deficit(imbalance, delta)
+#else //WITHOUT_CAPACITY_SCALING
+      dijkstra(res_graph, red_cost), pred(res_graph),
+      has_deficit(imbalance)
+#endif
+    {
+      // Removing nonzero lower bounds
+      capacity = subMap(_capacity, _lower);
+      Supply sum = 0;
+      for (NodeIt n(graph); n != INVALID; ++n) {
+	Supply s = _supply[n];
+	for (InEdgeIt e(graph, n); e != INVALID; ++e)
+	  s += _lower[e];
+	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
+	  s -= _lower[e];
+	supply[n] = imbalance[n] = s;
+	sum += s;
+      }
+      valid_supply = sum == 0;
+    }
+
+    /// \brief General constructor of the class (without lower bounds).
+    ///
+    /// General constructor of the class (without lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _supply The supply values of the nodes (signed).
+    CapacityScaling( const Graph &_graph,
+		     const CapacityMap &_capacity,
+		     const CostMap &_cost,
+		     const SupplyMap &_supply ) :
+      graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
+      supply(_supply), flow(_graph, 0), potential(_graph, 0),
+      res_graph(_graph, capacity, flow),
+      red_cost(res_graph, cost, potential), imbalance(_graph),
+#ifdef WITH_SCALING
+      delta(0), delta_filter(res_graph, delta),
+      dres_graph(res_graph, delta_filter),
+      dijkstra(dres_graph, red_cost), pred(dres_graph),
+      delta_deficit(imbalance, delta)
+#else //WITHOUT_CAPACITY_SCALING
+      dijkstra(res_graph, red_cost), pred(res_graph),
+      has_deficit(imbalance)
+#endif
+    {
+      // Checking the sum of supply values
+      Supply sum = 0;
+      for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
+      valid_supply = sum == 0;
+    }
+
+    /// \brief Simple constructor of the class (with lower bounds).
+    ///
+    /// Simple constructor of the class (with lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _lower The lower bounds of the edges.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _s The source node.
+    /// \param _t The target node.
+    /// \param _flow_value 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).
+    CapacityScaling( const Graph &_graph,
+		     const LowerMap &_lower,
+		     const CapacityMap &_capacity,
+		     const CostMap &_cost,
+		     Node _s, Node _t,
+		     Supply _flow_value ) :
+      graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
+      supply(_graph), flow(_graph, 0), potential(_graph, 0),
+      res_graph(_graph, capacity, flow),
+      red_cost(res_graph, cost, potential), imbalance(_graph),
+#ifdef WITH_SCALING
+      delta(0), delta_filter(res_graph, delta),
+      dres_graph(res_graph, delta_filter),
+      dijkstra(dres_graph, red_cost), pred(dres_graph),
+      delta_deficit(imbalance, delta)
+#else //WITHOUT_CAPACITY_SCALING
+      dijkstra(res_graph, red_cost), pred(res_graph),
+      has_deficit(imbalance)
+#endif
+    {
+      // Removing nonzero lower bounds
+      capacity = subMap(_capacity, _lower);
+      for (NodeIt n(graph); n != INVALID; ++n) {
+	Supply s = 0;
+	if (n == _s) s =  _flow_value;
+	if (n == _t) s = -_flow_value;
+	for (InEdgeIt e(graph, n); e != INVALID; ++e)
+	  s += _lower[e];
+	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
+	  s -= _lower[e];
+	supply[n] = imbalance[n] = s;
+      }
+      valid_supply = true;
+    }
+
+    /// \brief Simple constructor of the class (without lower bounds).
+    ///
+    /// Simple constructor of the class (without lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _s The source node.
+    /// \param _t The target node.
+    /// \param _flow_value 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).
+    CapacityScaling( const Graph &_graph,
+		     const CapacityMap &_capacity,
+		     const CostMap &_cost,
+		     Node _s, Node _t,
+		     Supply _flow_value ) :
+      graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
+      supply(_graph, 0), flow(_graph, 0), potential(_graph, 0),
+      res_graph(_graph, capacity, flow),
+      red_cost(res_graph, cost, potential), imbalance(_graph),
+#ifdef WITH_SCALING
+      delta(0), delta_filter(res_graph, delta),
+      dres_graph(res_graph, delta_filter),
+      dijkstra(dres_graph, red_cost), pred(dres_graph),
+      delta_deficit(imbalance, delta)
+#else //WITHOUT_CAPACITY_SCALING
+      dijkstra(res_graph, red_cost), pred(res_graph),
+      has_deficit(imbalance)
+#endif
+    {
+      supply[_s] =  _flow_value;
+      supply[_t] = -_flow_value;
+      valid_supply = true;
+    }
+
+    /// \brief Returns a const reference to the flow map.
+    ///
+    /// Returns a const reference to the flow map.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    const FlowMap& flowMap() const {
+      return flow;
+    }
+
+    /// \brief Returns a const reference to the potential map (the dual
+    /// solution).
+    ///
+    /// Returns a const reference to the potential map (the dual
+    /// solution).
+    ///
+    /// \pre \ref run() must be called before using this function.
+    const PotentialMap& potentialMap() const {
+      return potential;
+    }
+
+    /// \brief Returns the total cost of the found flow.
+    ///
+    /// Returns the total cost of the found flow. The complexity of the
+    /// function is \f$ O(e) \f$.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost totalCost() const {
+      Cost c = 0;
+      for (EdgeIt e(graph); e != INVALID; ++e)
+	c += flow[e] * cost[e];
+      return c;
+    }
+
+    /// \brief Runs the successive shortest path algorithm.
+    ///
+    /// Runs the successive shortest path algorithm.
+    ///
+    /// \return \c true if a feasible flow can be found.
+    bool run() {
+      return init() && start();
+    }
+
+  protected:
+
+    /// \brief Initializes the algorithm.
+    bool init() {
+      if (!valid_supply) return false;
+
+      // Initalizing Dijkstra class
+      updater.potentialMap(potential);
+      dijkstra.distMap(updater).predMap(pred);
+
+#ifdef WITH_SCALING
+      // Initilaizing delta value
+      Capacity max_cap = 0;
+      for (EdgeIt e(graph); e != INVALID; ++e) {
+	if (capacity[e] > max_cap) max_cap = capacity[e];
+      }
+      for (delta = 1; 2 * delta < max_cap; delta *= 2) ;
+#endif
+      return true;
+    }
+
+#ifdef WITH_SCALING
+    /// \brief Executes the capacity scaling version of the successive
+    /// shortest path algorithm.
+    bool start() {
+      typedef typename DeltaResGraph::EdgeIt DeltaResEdgeIt;
+      typedef typename DeltaResGraph::Edge DeltaResEdge;
+
+      // Processing capacity scaling phases
+      ResNode s, t;
+      for ( ; delta >= 1; delta = delta < 4 && delta > 1 ?
+				  1 : delta / 4 )
+      {
+	// Saturating edges not satisfying the optimality condition
+	Capacity r;
+	for (DeltaResEdgeIt e(dres_graph); e != INVALID; ++e) {
+	  if (red_cost[e] < 0) {
+	    res_graph.augment(e, r = res_graph.rescap(e));
+	    imbalance[dres_graph.target(e)] += r;
+	    imbalance[dres_graph.source(e)] -= r;
+	  }
+	}
+
+	// Finding excess nodes
+	excess_nodes.clear();
+	for (ResNodeIt n(res_graph); n != INVALID; ++n) {
+	  if (imbalance[n] >= delta) excess_nodes.push_back(n);
+	}
+	next_node = 0;
+
+	// Finding successive shortest paths
+	while (next_node < excess_nodes.size()) {
+	  // Running Dijkstra
+	  s = excess_nodes[next_node];
+	  updater.init();
+	  dijkstra.init();
+	  dijkstra.addSource(s);
+	  if ((t = dijkstra.start(delta_deficit)) == INVALID) {
+	    if (delta > 1) {
+	      ++next_node;
+	      continue;
+	    }
+	    return false;
+	  }
+
+	  // Updating node potentials
+	  updater.update();
+
+	  // Augment along a shortest path from s to t
+	  Capacity d = imbalance[s] < -imbalance[t] ?
+	    imbalance[s] : -imbalance[t];
+	  ResNode u = t;
+	  ResEdge e;
+	  if (d > delta) {
+	    while ((e = pred[u]) != INVALID) {
+	      if (res_graph.rescap(e) < d) d = res_graph.rescap(e);
+	      u = dres_graph.source(e);
+	    }
+	  }
+	  u = t;
+	  while ((e = pred[u]) != INVALID) {
+	    res_graph.augment(e, d);
+	    u = dres_graph.source(e);
+	  }
+	  imbalance[s] -= d;
+	  imbalance[t] += d;
+	  if (imbalance[s] < delta) ++next_node;
+	}
+      }
+
+      // Handling nonzero lower bounds
+      if (lower) {
+	for (EdgeIt e(graph); e != INVALID; ++e)
+	  flow[e] += (*lower)[e];
+      }
+      return true;
+    }
+
+#else //WITHOUT_CAPACITY_SCALING
+    /// \brief Executes the successive shortest path algorithm without
+    /// capacity scaling.
+    bool start() {
+      // Finding excess nodes
+      for (ResNodeIt n(res_graph); n != INVALID; ++n) {
+	if (imbalance[n] > 0) excess_nodes.push_back(n);
+      }
+      if (excess_nodes.size() == 0) return true;
+      next_node = 0;
+
+      // Finding successive shortest paths
+      ResNode s, t;
+      while ( imbalance[excess_nodes[next_node]] > 0 ||
+	      ++next_node < excess_nodes.size() )
+      {
+	// Running Dijkstra
+	s = excess_nodes[next_node];
+	updater.init();
+	dijkstra.init();
+	dijkstra.addSource(s);
+	if ((t = dijkstra.start(has_deficit)) == INVALID)
+	  return false;
+
+	// Updating node potentials
+	updater.update();
+
+	// Augmenting along a shortest path from s to t
+	Capacity delta = imbalance[s] < -imbalance[t] ?
+	  imbalance[s] : -imbalance[t];
+	ResNode u = t;
+	ResEdge e;
+	while ((e = pred[u]) != INVALID) {
+	  if (res_graph.rescap(e) < delta) delta = res_graph.rescap(e);
+	  u = res_graph.source(e);
+	}
+	u = t;
+	while ((e = pred[u]) != INVALID) {
+	  res_graph.augment(e, delta);
+	  u = res_graph.source(e);
+	}
+	imbalance[s] -= delta;
+	imbalance[t] += delta;
+      }
+
+      // Handling nonzero lower bounds
+      if (lower) {
+	for (EdgeIt e(graph); e != INVALID; ++e)
+	  flow[e] += (*lower)[e];
+      }
+      return true;
+    }
+#endif
+
+  }; //class CapacityScaling
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_CAPACITY_SCALING_H
diff -r 3f1c7a6c33cd -r c9218405595b lemon/cycle_canceling.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/cycle_canceling.h	Mon May 07 11:42:18 2007 +0000
@@ -0,0 +1,509 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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
+///
+/// \file
+/// \brief A cycle canceling algorithm for finding a minimum cost flow.
+
+#include <vector>
+#include <lemon/circulation.h>
+#include <lemon/graph_adaptor.h>
+
+/// \brief The used cycle canceling method.
+#define LIMITED_CYCLE_CANCELING
+//#define MIN_MEAN_CYCLE_CANCELING
+
+#ifdef LIMITED_CYCLE_CANCELING
+  #include <lemon/bellman_ford.h>
+  /// \brief The maximum number of iterations for the first execution
+  /// of the \ref lemon::BellmanFord "Bellman-Ford" algorithm.
+  #define STARTING_LIMIT	2
+  /// \brief The iteration limit for the
+  /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
+  /// <tt>ALPHA_MUL % ALPHA_DIV</tt> in every round.
+  #define ALPHA_MUL		3
+  /// \brief The iteration limit for the
+  /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
+  /// <tt>ALPHA_MUL % ALPHA_DIV</tt> in every round.
+  #define ALPHA_DIV		2
+
+//#define _ONLY_ONE_CYCLE_
+//#define _DEBUG_ITER_
+#endif
+
+#ifdef MIN_MEAN_CYCLE_CANCELING
+  #include <lemon/min_mean_cycle.h>
+  #include <lemon/path.h>
+#endif
+
+namespace lemon {
+
+  /// \addtogroup min_cost_flow
+  /// @{
+
+  /// \brief Implementation of a cycle canceling algorithm for finding
+  /// a minimum cost flow.
+  ///
+  /// \ref lemon::CycleCanceling "CycleCanceling" implements a cycle
+  /// canceling algorithm for finding a minimum cost flow.
+  ///
+  /// \param Graph The directed graph type the algorithm runs on.
+  /// \param LowerMap The type of the lower bound map.
+  /// \param CapacityMap The type of the capacity (upper bound) map.
+  /// \param CostMap The type of the cost (length) map.
+  /// \param SupplyMap The type of the supply map.
+  ///
+  /// \warning
+  /// - Edge capacities and costs should be nonnegative integers.
+  ///	However \c CostMap::Value should be signed type.
+  /// - Supply values should be integers.
+  /// - \c LowerMap::Value must be convertible to
+  ///	\c CapacityMap::Value and \c CapacityMap::Value must be
+  ///	convertible to \c SupplyMap::Value.
+  ///
+  /// \author Peter Kovacs
+
+template < typename Graph,
+	   typename LowerMap = typename Graph::template EdgeMap<int>,
+	   typename CapacityMap = LowerMap,
+	   typename CostMap = typename Graph::template EdgeMap<int>,
+	   typename SupplyMap = typename Graph::template NodeMap
+				<typename CapacityMap::Value> >
+  class CycleCanceling
+  {
+    typedef typename Graph::Node Node;
+    typedef typename Graph::NodeIt NodeIt;
+    typedef typename Graph::Edge Edge;
+    typedef typename Graph::EdgeIt EdgeIt;
+    typedef typename Graph::InEdgeIt InEdgeIt;
+    typedef typename Graph::OutEdgeIt OutEdgeIt;
+
+    typedef typename LowerMap::Value Lower;
+    typedef typename CapacityMap::Value Capacity;
+    typedef typename CostMap::Value Cost;
+    typedef typename SupplyMap::Value Supply;
+    typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
+    typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
+
+    typedef ResGraphAdaptor< const Graph, Capacity,
+			     CapacityRefMap, CapacityRefMap > ResGraph;
+    typedef typename ResGraph::Node ResNode;
+    typedef typename ResGraph::NodeIt ResNodeIt;
+    typedef typename ResGraph::Edge ResEdge;
+    typedef typename ResGraph::EdgeIt ResEdgeIt;
+
+  public:
+
+    /// \brief The type of the flow map.
+    typedef CapacityRefMap FlowMap;
+
+  protected:
+
+    /// \brief Map adaptor class for demand map.
+    class DemandMap : public MapBase<Node, Supply>
+    {
+    private:
+
+      const SupplyMap *map;
+
+    public:
+
+      typedef typename MapBase<Node, Supply>::Value Value;
+      typedef typename MapBase<Node, Supply>::Key Key;
+
+      DemandMap(const SupplyMap &_map) : map(&_map) {}
+
+      Value operator[](const Key &e) const {
+	return -(*map)[e];
+      }
+
+    }; //class DemandMap
+
+    /// \brief Map adaptor class for handling residual edge costs.
+    class ResCostMap : public MapBase<ResEdge, Cost>
+    {
+    private:
+
+      const CostMap &cost_map;
+
+    public:
+
+      typedef typename MapBase<ResEdge, Cost>::Value Value;
+      typedef typename MapBase<ResEdge, Cost>::Key Key;
+
+      ResCostMap(const CostMap &_cost) : cost_map(_cost) {}
+
+      Value operator[](const Key &e) const {
+	return ResGraph::forward(e) ? cost_map[e] : -cost_map[e];
+      }
+
+    }; //class ResCostMap
+
+  protected:
+
+    /// \brief The directed graph the algorithm runs on.
+    const Graph &graph;
+    /// \brief The original lower bound map.
+    const LowerMap *lower;
+    /// \brief The modified capacity map.
+    CapacityRefMap capacity;
+    /// \brief The cost map.
+    const CostMap &cost;
+    /// \brief The modified supply map.
+    SupplyRefMap supply;
+    /// \brief The modified demand map.
+    DemandMap demand;
+    /// \brief The sum of supply values equals zero.
+    bool valid_supply;
+
+    /// \brief The current flow.
+    FlowMap flow;
+    /// \brief The residual graph.
+    ResGraph res_graph;
+    /// \brief The residual cost map.
+    ResCostMap res_cost;
+
+  public :
+
+    /// \brief General constructor of the class (with lower bounds).
+    ///
+    /// General constructor of the class (with lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _lower The lower bounds of the edges.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _supply The supply values of the nodes (signed).
+    CycleCanceling( const Graph &_graph,
+		    const LowerMap &_lower,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    const SupplyMap &_supply ) :
+      graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
+      supply(_graph), demand(supply), flow(_graph, 0),
+      res_graph(_graph, capacity, flow), res_cost(cost)
+    {
+      // Removing nonzero lower bounds
+      capacity = subMap(_capacity, _lower);
+      Supply sum = 0;
+      for (NodeIt n(graph); n != INVALID; ++n) {
+	Supply s = _supply[n];
+	for (InEdgeIt e(graph, n); e != INVALID; ++e)
+	  s += _lower[e];
+	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
+	  s -= _lower[e];
+	sum += (supply[n] = s);
+      }
+      valid_supply = sum == 0;
+    }
+
+    /// \brief General constructor of the class (without lower bounds).
+    ///
+    /// General constructor of the class (without lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _supply The supply values of the nodes (signed).
+    CycleCanceling( const Graph &_graph,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    const SupplyMap &_supply ) :
+      graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
+      supply(_supply), demand(supply), flow(_graph, 0),
+      res_graph(_graph, capacity, flow), res_cost(cost)
+    {
+      // Checking the sum of supply values
+      Supply sum = 0;
+      for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
+      valid_supply = sum == 0;
+    }
+
+
+    /// \brief Simple constructor of the class (with lower bounds).
+    ///
+    /// Simple constructor of the class (with lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _lower The lower bounds of the edges.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _s The source node.
+    /// \param _t The target node.
+    /// \param _flow_value 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).
+    CycleCanceling( const Graph &_graph,
+		    const LowerMap &_lower,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    Node _s, Node _t,
+		    Supply _flow_value ) :
+      graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
+      supply(_graph), demand(supply), flow(_graph, 0),
+      res_graph(_graph, capacity, flow), res_cost(cost)
+    {
+      // Removing nonzero lower bounds
+      capacity = subMap(_capacity, _lower);
+      for (NodeIt n(graph); n != INVALID; ++n) {
+	Supply s = 0;
+	if (n == _s) s =  _flow_value;
+	if (n == _t) s = -_flow_value;
+	for (InEdgeIt e(graph, n); e != INVALID; ++e)
+	  s += _lower[e];
+	for (OutEdgeIt e(graph, n); e != INVALID; ++e)
+	  s -= _lower[e];
+	supply[n] = s;
+      }
+      valid_supply = true;
+    }
+
+    /// \brief Simple constructor of the class (without lower bounds).
+    ///
+    /// Simple constructor of the class (without lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _s The source node.
+    /// \param _t The target node.
+    /// \param _flow_value 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).
+    CycleCanceling( const Graph &_graph,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    Node _s, Node _t,
+		    Supply _flow_value ) :
+      graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
+      supply(_graph, 0), demand(supply), flow(_graph, 0),
+      res_graph(_graph, capacity, flow), res_cost(cost)
+    {
+      supply[_s] =  _flow_value;
+      supply[_t] = -_flow_value;
+      valid_supply = true;
+    }
+
+    /// \brief Returns a const reference to the flow map.
+    ///
+    /// Returns a const reference to the flow map.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    const FlowMap& flowMap() const {
+      return flow;
+    }
+
+    /// \brief Returns the total cost of the found flow.
+    ///
+    /// Returns the total cost of the found flow. The complexity of the
+    /// function is \f$ O(e) \f$.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost totalCost() const {
+      Cost c = 0;
+      for (EdgeIt e(graph); e != INVALID; ++e)
+	c += flow[e] * cost[e];
+      return c;
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// Runs the algorithm.
+    ///
+    /// \return \c true if a feasible flow can be found.
+    bool run() {
+      return init() && start();
+    }
+
+  protected:
+
+    /// \brief Initializes the algorithm.
+    bool init() {
+      // Checking the sum of supply values
+      Supply sum = 0;
+      for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
+      if (sum != 0) return false;
+
+      // Finding a feasible flow
+      Circulation< Graph, Capacity, FlowMap, ConstMap<Edge, Capacity>,
+		   CapacityRefMap, DemandMap >
+	circulation( graph, constMap<Edge>((Capacity)0),
+		     capacity, demand, flow );
+      return circulation.run() == -1;
+    }
+
+#ifdef LIMITED_CYCLE_CANCELING
+    /// \brief Executes a cycle canceling algorithm using
+    /// \ref lemon::BellmanFord "Bellman-Ford" algorithm with limited
+    /// iteration count.
+    bool start() {
+      typename BellmanFord<ResGraph, ResCostMap>::PredMap pred(res_graph);
+      typename ResGraph::template NodeMap<int> visited(res_graph);
+      std::vector<ResEdge> cycle;
+      int node_num = countNodes(graph);
+
+#ifdef _DEBUG_ITER_
+      int cycle_num = 0;
+#endif
+      int length_bound = STARTING_LIMIT;
+      bool optimal = false;
+      while (!optimal) {
+	BellmanFord<ResGraph, ResCostMap> bf(res_graph, res_cost);
+	bf.predMap(pred);
+	bf.init(0);
+	int iter_num = 0;
+	bool cycle_found = false;
+	while (!cycle_found) {
+	  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 = true;
+	    break;
+	  } else {
+	    // Searching for node disjoint negative cycles
+	    for (ResNodeIt n(res_graph); n != INVALID; ++n)
+	      visited[n] = 0;
+	    int id = 0;
+	    for (ResNodeIt n(res_graph); n != INVALID; ++n) {
+	      if (visited[n] > 0) continue;
+	      visited[n] = ++id;
+	      ResNode u = pred[n] == INVALID ?
+			  INVALID : res_graph.source(pred[n]);
+	      while (u != INVALID && visited[u] == 0) {
+		visited[u] = id;
+		u = pred[u] == INVALID ?
+		    INVALID : res_graph.source(pred[u]);
+	      }
+	      if (u != INVALID && visited[u] == id) {
+		// Finding the negative cycle
+		cycle_found = true;
+		cycle.clear();
+		ResEdge e = pred[u];
+		cycle.push_back(e);
+		Capacity d = res_graph.rescap(e);
+		while (res_graph.source(e) != u) {
+		  cycle.push_back(e = pred[res_graph.source(e)]);
+		  if (res_graph.rescap(e) < d)
+		    d = res_graph.rescap(e);
+		}
+#ifdef _DEBUG_ITER_
+		++cycle_num;
+#endif
+		// Augmenting along the cycle
+		for (int i = 0; i < cycle.size(); ++i)
+		  res_graph.augment(cycle[i], d);
+#ifdef _ONLY_ONE_CYCLE_
+		break;
+#endif
+	      }
+	    }
+	  }
+
+	  if (!cycle_found)
+	    length_bound = length_bound * ALPHA_MUL / ALPHA_DIV;
+	}
+      }
+
+#ifdef _DEBUG_ITER_
+      std::cout << "Limited cycle canceling algorithm finished. "
+		<< "Found " << cycle_num << " negative cycles."
+		<< std::endl;
+#endif
+
+      // Handling nonzero lower bounds
+      if (lower) {
+	for (EdgeIt e(graph); e != INVALID; ++e)
+	  flow[e] += (*lower)[e];
+      }
+      return true;
+    }
+#endif
+
+#ifdef MIN_MEAN_CYCLE_CANCELING
+    /// \brief Executes the minimum mean cycle canceling algorithm
+    /// using \ref lemon::MinMeanCycle "MinMeanCycle" class.
+    bool start() {
+      typedef Path<ResGraph> ResPath;
+      MinMeanCycle<ResGraph, ResCostMap> mmc(res_graph, res_cost);
+      ResPath cycle;
+
+#ifdef _DEBUG_ITER_
+      int cycle_num = 0;
+#endif
+      mmc.cyclePath(cycle).init();
+      if (mmc.findMinMean()) {
+	while (mmc.cycleLength() < 0) {
+#ifdef _DEBUG_ITER_
+	  ++iter;
+#endif
+	  // Finding the cycle
+	  mmc.findCycle();
+
+	  // Finding the largest flow amount that can be augmented
+	  // along the cycle
+	  Capacity delta = 0;
+	  for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e) {
+	    if (delta == 0 || res_graph.rescap(e) < delta)
+	      delta = res_graph.rescap(e);
+	  }
+
+	  // Augmenting along the cycle
+	  for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e)
+	    res_graph.augment(e, delta);
+
+	  // Finding the minimum cycle mean for the modified residual
+	  // graph
+	  mmc.reset();
+	  if (!mmc.findMinMean()) break;
+	}
+      }
+
+#ifdef _DEBUG_ITER_
+      std::cout << "Minimum mean cycle canceling algorithm finished. "
+		<< "Found " << cycle_num << " negative cycles."
+		<< std::endl;
+#endif
+
+      // Handling nonzero lower bounds
+      if (lower) {
+	for (EdgeIt e(graph); e != INVALID; ++e)
+	  flow[e] += (*lower)[e];
+      }
+      return true;
+    }
+#endif
+
+  }; //class CycleCanceling
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_CYCLE_CANCELING_H
diff -r 3f1c7a6c33cd -r c9218405595b lemon/min_cost_flow.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/min_cost_flow.h	Mon May 07 11:42:18 2007 +0000
@@ -0,0 +1,127 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_MIN_COST_FLOW_H
+#define LEMON_MIN_COST_FLOW_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief An efficient algorithm for finding a minimum cost flow.
+
+#include <lemon/network_simplex.h>
+
+namespace lemon {
+
+  /// \addtogroup min_cost_flow
+  /// @{
+
+  /// \brief An efficient algorithm for finding a minimum cost flow.
+  ///
+  /// \ref lemon::MinCostFlow "MinCostFlow" implements an efficient
+  /// algorithm for finding a minimum cost flow.
+  ///
+  /// \note \ref lemon::MinCostFlow "MinCostFlow" is just an alias for
+  /// \ref lemon::NetworkSimplex "NetworkSimplex", wich is the most
+  /// efficient implementation for the minimum cost flow problem in the
+  /// Lemon library according to our benchmark tests.
+  ///
+  /// \note There are three implementations for the minimum cost flow
+  /// problem, that can be used in the same way.
+  /// - \ref lemon::CapacityScaling "CapacityScaling"
+  /// - \ref lemon::CycleCanceling "CycleCanceling"
+  /// - \ref lemon::NetworkSimplex "NetworkSimplex"
+  ///
+  /// \param Graph The directed graph type the algorithm runs on.
+  /// \param LowerMap The type of the lower bound map.
+  /// \param CapacityMap The type of the capacity (upper bound) map.
+  /// \param CostMap The type of the cost (length) map.
+  /// \param SupplyMap The type of the supply map.
+  ///
+  /// \warning
+  /// - Edge capacities and costs should be nonnegative integers.
+  ///	However \c CostMap::Value should be signed type.
+  /// - Supply values should be integers.
+  /// - \c LowerMap::Value must be convertible to
+  ///	\c CapacityMap::Value and \c CapacityMap::Value must be
+  ///	convertible to \c SupplyMap::Value.
+  ///
+  /// \author Peter Kovacs
+
+template < typename Graph,
+	   typename LowerMap = typename Graph::template EdgeMap<int>,
+	   typename CapacityMap = LowerMap,
+	   typename CostMap = typename Graph::template EdgeMap<int>,
+	   typename SupplyMap = typename Graph::template NodeMap
+				<typename CapacityMap::Value> >
+  class MinCostFlow :
+    public NetworkSimplex< Graph,
+			   LowerMap,
+			   CapacityMap,
+			   CostMap,
+			   SupplyMap >
+  {
+    typedef NetworkSimplex< Graph, LowerMap, CapacityMap,
+			    CostMap, SupplyMap >
+      MinCostFlowImpl;
+    typedef typename Graph::Node Node;
+    typedef typename SupplyMap::Value Supply;
+
+  public:
+
+    /// \brief General constructor of the class (with lower bounds).
+    MinCostFlow( const Graph &_graph,
+		 const LowerMap &_lower,
+		 const CapacityMap &_capacity,
+		 const CostMap &_cost,
+		 const SupplyMap &_supply ) :
+      MinCostFlowImpl(_graph, _lower, _capacity, _cost, _supply) {}
+
+    /// \brief General constructor of the class (without lower bounds).
+    MinCostFlow( const Graph &_graph,
+		 const CapacityMap &_capacity,
+		 const CostMap &_cost,
+		 const SupplyMap &_supply ) :
+      MinCostFlowImpl(_graph, _capacity, _cost, _supply) {}
+
+    /// \brief Simple constructor of the class (with lower bounds).
+    MinCostFlow( const Graph &_graph,
+		 const LowerMap &_lower,
+		 const CapacityMap &_capacity,
+		 const CostMap &_cost,
+		 Node _s, Node _t,
+		 Supply _flow_value ) :
+      MinCostFlowImpl( _graph, _lower, _capacity, _cost,
+		       _s, _t, _flow_value ) {}
+
+    /// \brief Simple constructor of the class (without lower bounds).
+    MinCostFlow( const Graph &_graph,
+		 const CapacityMap &_capacity,
+		 const CostMap &_cost,
+		 Node _s, Node _t,
+		 Supply _flow_value ) :
+      MinCostFlowImpl( _graph, _capacity, _cost,
+		       _s, _t, _flow_value ) {}
+
+  }; //class MinCostFlow
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_MIN_COST_FLOW_H
diff -r 3f1c7a6c33cd -r c9218405595b lemon/min_cost_max_flow.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/min_cost_max_flow.h	Mon May 07 11:42:18 2007 +0000
@@ -0,0 +1,156 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_MIN_COST_MAX_FLOW_H
+#define LEMON_MIN_COST_MAX_FLOW_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief An efficient algorithm for finding a minimum cost maximum
+/// flow.
+
+#include <lemon/preflow.h>
+#include <lemon/network_simplex.h>
+#include <lemon/maps.h>
+
+namespace lemon {
+
+  /// \addtogroup min_cost_flow
+  /// @{
+
+  /// \brief An efficient algorithm for finding a minimum cost maximum
+  /// flow.
+  ///
+  /// \ref lemon::MinCostFlow "MinCostMaxFlow" implements an efficient
+  /// algorithm for finding a maximum flow having minimal total cost
+  /// from a given source node to a given target node in a directed
+  /// graph.
+  ///
+  /// \note \ref lemon::MinCostMaxFlow "MinCostMaxFlow" uses
+  /// \ref lemon::Preflow "Preflow" algorithm for finding the maximum
+  /// flow value and \ref lemon::NetworkSimplex "NetworkSimplex" 
+  /// algorithm for finding a minimum cost flow of that value.
+  ///
+  /// \param Graph The directed graph type the algorithm runs on.
+  /// \param CapacityMap The type of the capacity (upper bound) map.
+  /// \param CostMap The type of the cost (length) map.
+  ///
+  /// \warning
+  /// - Edge capacities and costs should be nonnegative integers.
+  ///	However \c CostMap::Value should be signed type.
+  ///
+  /// \author Peter Kovacs
+
+template < typename Graph,
+	   typename CapacityMap = typename Graph::template EdgeMap<int>,
+	   typename CostMap = typename Graph::template EdgeMap<int> >
+  class MinCostMaxFlow
+  {
+    typedef typename Graph::Node Node;
+    typedef typename Graph::Edge Edge;
+
+    typedef typename CapacityMap::Value Capacity;
+    typedef typename Graph::template NodeMap<Capacity> SupplyMap;
+    typedef NetworkSimplex< Graph, CapacityMap, CapacityMap,
+			    CostMap, SupplyMap >
+      MinCostFlowImpl;
+
+  public:
+
+    /// \brief The type of the flow map.
+    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
+
+  private:
+
+    /// \brief The directed graph the algorithm runs on.
+    const Graph &graph;
+    /// \brief The modified capacity map.
+    const CapacityMap &capacity;
+    /// \brief The cost map.
+    const CostMap &cost;
+    /// \brief The source node.
+    Node source;
+    /// \brief The target node.
+    Node target;
+    /// \brief The edge map of the found flow.
+    FlowMap flow;
+
+    typedef Preflow<Graph, Capacity, CapacityMap, FlowMap> PreflowImpl;
+    /// \brief \ref lemon::Preflow "Preflow" class for finding the
+    /// maximum flow value.
+    PreflowImpl preflow;
+
+  public:
+
+    /// \brief The constructor of the class.
+    ///
+    /// The constructor of the class.
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _s The source node.
+    /// \param _t The target node.
+    MinCostMaxFlow( const Graph &_graph,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    Node _s, Node _t ) :
+      graph(_graph), capacity(_capacity), cost(_cost),
+      source(_s), target(_t), flow(_graph),
+      preflow(_graph, _s, _t, _capacity, flow)
+    {}
+
+    /// \brief Returns a const reference to the flow map.
+    ///
+    /// Returns a const reference to the flow map.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    const FlowMap& flowMap() const {
+      return flow;
+    }
+
+    /// \brief Returns the total cost of the found flow.
+    ///
+    /// Returns the total cost of the found flow. The complexity of the
+    /// function is \f$ O(e) \f$.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost totalCost() const {
+      Cost c = 0;
+      for (EdgeIt e(graph); e != INVALID; ++e)
+	c += flow[e] * cost[e];
+      return c;
+    }
+
+    /// \brief Runs the algorithm.
+    void run() {
+      preflow.phase1();
+      MinCostFlowImpl mcf_impl( graph, capacity, cost,
+				source, target, preflow.flowValue() );
+      mcf_impl.run();
+      flow = mcf_impl.flowMap();
+    }
+
+  }; //class MinCostMaxFlow
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_MIN_COST_MAX_FLOW_H
diff -r 3f1c7a6c33cd -r c9218405595b lemon/network_simplex.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/network_simplex.h	Mon May 07 11:42:18 2007 +0000
@@ -0,0 +1,945 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_NETWORK_SIMPLEX_H
+#define LEMON_NETWORK_SIMPLEX_H
+
+/// \ingroup min_cost_flow
+///
+/// \file
+/// \brief The network simplex algorithm for finding a minimum cost
+/// flow.
+
+#include <limits>
+#include <lemon/smart_graph.h>
+#include <lemon/graph_utils.h>
+
+/// \brief The pivot rule used in the algorithm.
+#define EDGE_BLOCK_PIVOT
+//#define FIRST_ELIGIBLE_PIVOT
+//#define BEST_ELIGIBLE_PIVOT
+//#define CANDIDATE_LIST_PIVOT
+//#define SORTED_LIST_PIVOT
+
+/// \brief State constant for edges at their lower bounds.
+#define LOWER	1
+/// \brief State constant for edges in the spanning tree.
+#define TREE	0
+/// \brief State constant for edges at their upper bounds.
+#define UPPER	-1
+
+#ifdef EDGE_BLOCK_PIVOT
+  /// \brief Number of blocks for the "Edge Block" pivot rule.
+  #define BLOCK_NUM	  100
+  /// \brief Lower bound for the number of edges to use "Edge Block"
+  // pivot rule instead of "First Eligible" pivot rule.
+  #define MIN_BOUND	  1000
+#endif
+
+#ifdef CANDIDATE_LIST_PIVOT
+  #include <list>
+  /// \brief The maximum length of the edge list for the
+  /// "Candidate List" pivot rule.
+  #define LIST_LENGTH	  100
+  /// \brief The maximum number of minor iterations between two major
+  /// itarations.
+  #define MINOR_LIMIT	  50
+#endif
+
+#ifdef SORTED_LIST_PIVOT
+  #include <deque>
+  #include <algorithm>
+  /// \brief The maximum length of the edge list for the
+  /// "Sorted List" pivot rule.
+  #define LIST_LENGTH	  500
+  #define LOWER_DIV	  4
+#endif
+
+//#define _DEBUG_ITER_
+
+namespace lemon {
+
+  /// \addtogroup min_cost_flow
+  /// @{
+
+  /// \brief Implementation of the network simplex algorithm for
+  /// finding a minimum cost flow.
+  ///
+  /// \ref lemon::NetworkSimplex "NetworkSimplex" implements the
+  /// network simplex algorithm for finding a minimum cost flow.
+  ///
+  /// \param Graph The directed graph type the algorithm runs on.
+  /// \param LowerMap The type of the lower bound map.
+  /// \param CapacityMap The type of the capacity (upper bound) map.
+  /// \param CostMap The type of the cost (length) map.
+  /// \param SupplyMap The type of the supply map.
+  ///
+  /// \warning
+  /// - Edge capacities and costs should be nonnegative integers.
+  ///	However \c CostMap::Value should be signed type.
+  /// - Supply values should be integers.
+  /// - \c LowerMap::Value must be convertible to
+  ///	\c CapacityMap::Value and \c CapacityMap::Value must be
+  ///	convertible to \c SupplyMap::Value.
+  ///
+  /// \author Peter Kovacs
+
+template < typename Graph,
+	   typename LowerMap = typename Graph::template EdgeMap<int>,
+	   typename CapacityMap = LowerMap,
+	   typename CostMap = typename Graph::template EdgeMap<int>,
+	   typename SupplyMap = typename Graph::template NodeMap
+				<typename CapacityMap::Value> >
+  class NetworkSimplex
+  {
+    typedef typename LowerMap::Value Lower;
+    typedef typename CapacityMap::Value Capacity;
+    typedef typename CostMap::Value Cost;
+    typedef typename SupplyMap::Value Supply;
+
+    typedef SmartGraph SGraph;
+    typedef typename SGraph::Node Node;
+    typedef typename SGraph::NodeIt NodeIt;
+    typedef typename SGraph::Edge Edge;
+    typedef typename SGraph::EdgeIt EdgeIt;
+    typedef typename SGraph::InEdgeIt InEdgeIt;
+    typedef typename SGraph::OutEdgeIt OutEdgeIt;
+
+    typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
+    typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
+    typedef typename SGraph::template EdgeMap<Cost> SCostMap;
+    typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
+    typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
+
+    typedef typename SGraph::template NodeMap<int> IntNodeMap;
+    typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
+    typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
+    typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
+    typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
+
+    typedef typename Graph::template NodeMap<Node> NodeRefMap;
+    typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
+
+  public:
+
+    /// \brief The type of the flow map.
+    typedef typename Graph::template EdgeMap<Capacity> FlowMap;
+    /// \brief The type of the potential map.
+    typedef typename Graph::template NodeMap<Cost> PotentialMap;
+
+  protected:
+
+    /// \brief Map adaptor class for handling reduced edge costs.
+    class ReducedCostMap : public MapBase<Edge, Cost>
+    {
+    private:
+
+      const SGraph &gr;
+      const SCostMap &cost_map;
+      const SPotentialMap &pot_map;
+
+    public:
+
+      typedef typename MapBase<Edge, Cost>::Value Value;
+      typedef typename MapBase<Edge, Cost>::Key Key;
+
+      ReducedCostMap( const SGraph &_gr,
+		      const SCostMap &_cm,
+		      const SPotentialMap &_pm ) :
+	gr(_gr), cost_map(_cm), pot_map(_pm) {}
+
+      Value operator[](const Key &e) const {
+	return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
+      }
+
+    }; //class ReducedCostMap
+
+  protected:
+
+    /// \brief The directed graph the algorithm runs on.
+    SGraph graph;
+    /// \brief The original graph.
+    const Graph &graph_ref;
+    /// \brief The original lower bound map.
+    const LowerMap *lower;
+    /// \brief The capacity map.
+    SCapacityMap capacity;
+    /// \brief The cost map.
+    SCostMap cost;
+    /// \brief The supply map.
+    SSupplyMap supply;
+    /// \brief The reduced cost map.
+    ReducedCostMap red_cost;
+    /// \brief The sum of supply values equals zero.
+    bool valid_supply;
+
+    /// \brief The edge map of the current flow.
+    SCapacityMap flow;
+    /// \brief The edge map of the found flow on the original graph.
+    FlowMap flow_result;
+    /// \brief The potential node map.
+    SPotentialMap potential;
+    /// \brief The potential node map on the original graph.
+    PotentialMap potential_result;
+
+    /// \brief Node reference for the original graph.
+    NodeRefMap node_ref;
+    /// \brief Edge reference for the original graph.
+    EdgeRefMap edge_ref;
+
+    /// \brief The depth node map of the spanning tree structure.
+    IntNodeMap depth;
+    /// \brief The parent node map of the spanning tree structure.
+    NodeNodeMap parent;
+    /// \brief The pred_edge node map of the spanning tree structure.
+    EdgeNodeMap pred_edge;
+    /// \brief The thread node map of the spanning tree structure.
+    NodeNodeMap thread;
+    /// \brief The forward node map of the spanning tree structure.
+    BoolNodeMap forward;
+    /// \brief The state edge map.
+    IntEdgeMap state;
+
+
+#ifdef EDGE_BLOCK_PIVOT
+    /// \brief The size of blocks for the "Edge Block" pivot rule.
+    int block_size;
+#endif
+#ifdef CANDIDATE_LIST_PIVOT
+    /// \brief The list of candidate edges for the "Candidate List"
+    /// pivot rule.
+    std::list<Edge> candidates;
+    /// \brief The number of minor iterations.
+    int minor_count;
+#endif
+#ifdef SORTED_LIST_PIVOT
+    /// \brief The list of candidate edges for the "Sorted List"
+    /// pivot rule.
+    std::deque<Edge> candidates;
+#endif
+
+    // Root node of the starting spanning tree.
+    Node root;
+    // The entering edge of the current pivot iteration.
+    Edge in_edge;
+    // Temporary nodes used in the current pivot iteration.
+    Node join, u_in, v_in, u_out, v_out;
+    Node right, first, second, last;
+    Node stem, par_stem, new_stem;
+    // The maximum augment amount along the cycle in the current pivot
+    // iteration.
+    Capacity delta;
+
+  public :
+
+    /// \brief General constructor of the class (with lower bounds).
+    ///
+    /// General constructor of the class (with lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _lower The lower bounds of the edges.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _supply The supply values of the nodes (signed).
+    NetworkSimplex( const Graph &_graph,
+		    const LowerMap &_lower,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    const SupplyMap &_supply ) :
+      graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
+      supply(graph), flow(graph), flow_result(_graph), potential(graph),
+      potential_result(_graph), depth(graph), parent(graph),
+      pred_edge(graph), thread(graph), forward(graph), state(graph),
+      node_ref(graph_ref), edge_ref(graph_ref),
+      red_cost(graph, cost, potential)
+    {
+      // Checking the sum of supply values
+      Supply sum = 0;
+      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
+	sum += _supply[n];
+      if (!(valid_supply = sum == 0)) return;
+
+      // Copying graph_ref to graph
+      copyGraph(graph, graph_ref)
+	.edgeMap(cost, _cost)
+	.nodeRef(node_ref)
+	.edgeRef(edge_ref)
+	.run();
+
+      // Removing nonzero lower bounds
+      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
+	capacity[edge_ref[e]] = _capacity[e] - _lower[e];
+      }
+      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
+	Supply s = _supply[n];
+	for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
+	  s += _lower[e];
+	for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
+	  s -= _lower[e];
+	supply[node_ref[n]] = s;
+      }
+    }
+
+    /// \brief General constructor of the class (without lower bounds).
+    ///
+    /// General constructor of the class (without lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _supply The supply values of the nodes (signed).
+    NetworkSimplex( const Graph &_graph,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    const SupplyMap &_supply ) :
+      graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
+      supply(graph), flow(graph), flow_result(_graph), potential(graph),
+      potential_result(_graph), depth(graph), parent(graph),
+      pred_edge(graph), thread(graph), forward(graph), state(graph),
+      node_ref(graph_ref), edge_ref(graph_ref),
+      red_cost(graph, cost, potential)
+    {
+      // Checking the sum of supply values
+      Supply sum = 0;
+      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
+	sum += _supply[n];
+      if (!(valid_supply = sum == 0)) return;
+
+      // Copying graph_ref to graph
+      copyGraph(graph, graph_ref)
+	.edgeMap(capacity, _capacity)
+	.edgeMap(cost, _cost)
+	.nodeMap(supply, _supply)
+	.nodeRef(node_ref)
+	.edgeRef(edge_ref)
+	.run();
+    }
+
+    /// \brief Simple constructor of the class (with lower bounds).
+    ///
+    /// Simple constructor of the class (with lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _lower The lower bounds of the edges.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _s The source node.
+    /// \param _t The target node.
+    /// \param _flow_value 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).
+    NetworkSimplex( const Graph &_graph,
+		    const LowerMap &_lower,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    typename Graph::Node _s,
+		    typename Graph::Node _t,
+		    typename SupplyMap::Value _flow_value ) :
+      graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
+      supply(graph), flow(graph), flow_result(_graph), potential(graph),
+      potential_result(_graph), depth(graph), parent(graph),
+      pred_edge(graph), thread(graph), forward(graph), state(graph),
+      node_ref(graph_ref), edge_ref(graph_ref),
+      red_cost(graph, cost, potential)
+    {
+      // Copying graph_ref to graph
+      copyGraph(graph, graph_ref)
+	.edgeMap(cost, _cost)
+	.nodeRef(node_ref)
+	.edgeRef(edge_ref)
+	.run();
+
+      // Removing nonzero lower bounds
+      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
+	capacity[edge_ref[e]] = _capacity[e] - _lower[e];
+      }
+      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
+	Supply s = 0;
+	if (n == _s) s =  _flow_value;
+	if (n == _t) s = -_flow_value;
+	for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
+	  s += _lower[e];
+	for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
+	  s -= _lower[e];
+	supply[node_ref[n]] = s;
+      }
+      valid_supply = true;
+    }
+
+    /// \brief Simple constructor of the class (without lower bounds).
+    ///
+    /// Simple constructor of the class (without lower bounds).
+    ///
+    /// \param _graph The directed graph the algorithm runs on.
+    /// \param _capacity The capacities (upper bounds) of the edges.
+    /// \param _cost The cost (length) values of the edges.
+    /// \param _s The source node.
+    /// \param _t The target node.
+    /// \param _flow_value 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).
+    NetworkSimplex( const Graph &_graph,
+		    const CapacityMap &_capacity,
+		    const CostMap &_cost,
+		    typename Graph::Node _s,
+		    typename Graph::Node _t,
+		    typename SupplyMap::Value _flow_value ) :
+      graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
+      supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),
+      potential_result(_graph), depth(graph), parent(graph),
+      pred_edge(graph), thread(graph), forward(graph), state(graph),
+      node_ref(graph_ref), edge_ref(graph_ref),
+      red_cost(graph, cost, potential)
+    {
+      // Copying graph_ref to graph
+      copyGraph(graph, graph_ref)
+	.edgeMap(capacity, _capacity)
+	.edgeMap(cost, _cost)
+	.nodeRef(node_ref)
+	.edgeRef(edge_ref)
+	.run();
+      supply[node_ref[_s]] =  _flow_value;
+      supply[node_ref[_t]] = -_flow_value;
+      valid_supply = true;
+    }
+
+    /// \brief Returns a const reference to the flow map.
+    ///
+    /// Returns a const reference to the flow map.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    const FlowMap& flowMap() const {
+      return flow_result;
+    }
+
+    /// \brief Returns a const reference to the potential map (the dual
+    /// solution).
+    ///
+    /// Returns a const reference to the potential map (the dual
+    /// solution).
+    ///
+    /// \pre \ref run() must be called before using this function.
+    const PotentialMap& potentialMap() const {
+      return potential_result;
+    }
+
+    /// \brief Returns the total cost of the found flow.
+    ///
+    /// Returns the total cost of the found flow. The complexity of the
+    /// function is \f$ O(e) \f$.
+    ///
+    /// \pre \ref run() must be called before using this function.
+    Cost totalCost() const {
+      Cost c = 0;
+      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
+	c += flow_result[e] * cost[edge_ref[e]];
+      return c;
+    }
+
+    /// \brief Runs the algorithm.
+    ///
+    /// Runs the algorithm.
+    ///
+    /// \return \c true if a feasible flow can be found.
+    bool run() {
+      return init() && start();
+    }
+
+  protected:
+
+    /// \brief Extends the underlaying graph and initializes all the
+    /// node and edge maps.
+    bool init() {
+      if (!valid_supply) return false;
+
+      // Initializing state and flow maps
+      for (EdgeIt e(graph); e != INVALID; ++e) {
+	flow[e] = 0;
+	state[e] = LOWER;
+      }
+
+      // Adding an artificial root node to the graph
+      root = graph.addNode();
+      parent[root] = INVALID;
+      pred_edge[root] = INVALID;
+      depth[root] = supply[root] = potential[root] = 0;
+
+      // Adding artificial edges to the graph and initializing the node
+      // maps of the spanning tree data structure
+      Supply sum = 0;
+      Node last = root;
+      Edge e;
+      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
+      for (NodeIt u(graph); u != INVALID; ++u) {
+	if (u == root) continue;
+	thread[last] = u;
+	last = u;
+	parent[u] = root;
+	depth[u] = 1;
+	sum += supply[u];
+	if (supply[u] >= 0) {
+	  e = graph.addEdge(u, root);
+	  flow[e] = supply[u];
+	  forward[u] = true;
+	  potential[u] = max_cost;
+	} else {
+	  e = graph.addEdge(root, u);
+	  flow[e] = -supply[u];
+	  forward[u] = false;
+	  potential[u] = -max_cost;
+	}
+	cost[e] = max_cost;
+	capacity[e] = std::numeric_limits<Capacity>::max();
+	state[e] = TREE;
+	pred_edge[u] = e;
+      }
+      thread[last] = root;
+
+#ifdef EDGE_BLOCK_PIVOT
+      // Initializing block_size for the edge block pivot rule
+      int edge_num = countEdges(graph);
+      block_size = edge_num > MIN_BOUND ? edge_num / BLOCK_NUM + 1 : 1;
+#endif
+#ifdef CANDIDATE_LIST_PIVOT
+      minor_count = 0;
+#endif
+
+      return sum == 0;
+    }
+
+#ifdef FIRST_ELIGIBLE_PIVOT
+    /// \brief Finds entering edge according to the "First Eligible"
+    /// pivot rule.
+    bool findEnteringEdge(EdgeIt &next_edge) {
+      for (EdgeIt e = next_edge; e != INVALID; ++e) {
+	if (state[e] * red_cost[e] < 0) {
+	  in_edge = e;
+	  next_edge = ++e;
+	  return true;
+	}
+      }
+      for (EdgeIt e(graph); e != next_edge; ++e) {
+	if (state[e] * red_cost[e] < 0) {
+	  in_edge = e;
+	  next_edge = ++e;
+	  return true;
+	}
+      }
+      return false;
+    }
+#endif
+
+#ifdef BEST_ELIGIBLE_PIVOT
+    /// \brief Finds entering edge according to the "Best Eligible"
+    /// pivot rule.
+    bool findEnteringEdge() {
+      Cost min = 0;
+      for (EdgeIt e(graph); e != INVALID; ++e) {
+	if (state[e] * red_cost[e] < min) {
+	  min = state[e] * red_cost[e];
+	  in_edge = e;
+	}
+      }
+      return min < 0;
+    }
+#endif
+
+#ifdef EDGE_BLOCK_PIVOT
+    /// \brief Finds entering edge according to the "Edge Block"
+    /// pivot rule.
+    bool findEnteringEdge(EdgeIt &next_edge) {
+      if (block_size == 1) {
+	// Performing first eligible selection
+	for (EdgeIt e = next_edge; e != INVALID; ++e) {
+	  if (state[e] * red_cost[e] < 0) {
+	    in_edge = e;
+	    next_edge = ++e;
+	    return true;
+	  }
+	}
+	for (EdgeIt e(graph); e != next_edge; ++e) {
+	  if (state[e] * red_cost[e] < 0) {
+	    in_edge = e;
+	    next_edge = ++e;
+	    return true;
+	  }
+	}
+	return false;
+      } else {
+	// Performing edge block selection
+	Cost curr, min = 0;
+	EdgeIt min_edge(graph);
+	int cnt = 0;
+	for (EdgeIt e = next_edge; e != INVALID; ++e) {
+	  if ((curr = state[e] * red_cost[e]) < min) {
+	    min = curr;
+	    min_edge = e;
+	  }
+	  if (++cnt == block_size) {
+	    if (min < 0) break;
+	    cnt = 0;
+	  }
+	}
+	if (!(min < 0)) {
+	  for (EdgeIt e(graph); e != next_edge; ++e) {
+	    if ((curr = state[e] * red_cost[e]) < min) {
+	      min = curr;
+	      min_edge = e;
+	    }
+	    if (++cnt == block_size) {
+	      if (min < 0) break;
+	      cnt = 0;
+	    }
+	  }
+	}
+	in_edge = min_edge;
+	if ((next_edge = ++min_edge) == INVALID)
+	  next_edge = EdgeIt(graph);
+	return min < 0;
+      }
+    }
+#endif
+
+#ifdef CANDIDATE_LIST_PIVOT
+    /// \brief Functor class for removing non-eligible edges from the
+    /// candidate list.
+    class RemoveFunc
+    {
+    private:
+      const IntEdgeMap &st;
+      const ReducedCostMap &rc;
+    public:
+      RemoveFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
+	st(_st), rc(_rc) {}
+      bool operator()(const Edge &e) {
+	return st[e] * rc[e] >= 0;
+      }
+    };
+
+    /// \brief Finds entering edge according to the "Candidate List"
+    /// pivot rule.
+    bool findEnteringEdge() {
+      static RemoveFunc remove_func(state, red_cost);
+      typedef typename std::list<Edge>::iterator ListIt;
+
+      candidates.remove_if(remove_func);
+      if (minor_count >= MINOR_LIMIT || candidates.size() == 0) {
+	// Major iteration
+	for (EdgeIt e(graph); e != INVALID; ++e) {
+	  if (state[e] * red_cost[e] < 0) {
+	    candidates.push_back(e);
+	    if (candidates.size() == LIST_LENGTH) break;
+	  }
+	}
+	if (candidates.size() == 0) return false;
+      }
+
+      // Minor iteration
+      ++minor_count;
+      Cost min = 0;
+      for (ListIt it = candidates.begin(); it != candidates.end(); ++it) {
+	if (state[*it] * red_cost[*it] < min) {
+	  min = state[*it] * red_cost[*it];
+	  in_edge = *it;
+	}
+      }
+      return true;
+    }
+#endif
+
+#ifdef SORTED_LIST_PIVOT
+    /// \brief Functor class to compare edges during sort of the
+    /// candidate list.
+    class SortFunc
+    {
+    private:
+      const IntEdgeMap &st;
+      const ReducedCostMap &rc;
+    public:
+      SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
+	st(_st), rc(_rc) {}
+      bool operator()(const Edge &e1, const Edge &e2) {
+	return st[e1] * rc[e1] < st[e2] * rc[e2];
+      }
+    };
+
+    /// \brief Finds entering edge according to the "Sorted List"
+    /// pivot rule.
+    bool findEnteringEdge() {
+      static SortFunc sort_func(state, red_cost);
+
+      // Minor iteration
+      while (candidates.size() > 0) {
+	in_edge = candidates.front();
+	candidates.pop_front();
+	if (state[in_edge] * red_cost[in_edge] < 0) return true;
+      }
+
+      // Major iteration
+      Cost curr, min = 0;
+      for (EdgeIt e(graph); e != INVALID; ++e) {
+	if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
+	  candidates.push_back(e);
+	  if (curr < min) min = curr;
+	  if (candidates.size() >= LIST_LENGTH) break;
+	}
+      }
+      if (candidates.size() == 0) return false;
+      sort(candidates.begin(), candidates.end(), sort_func);
+      return true;
+    }
+#endif
+
+    /// \brief Finds the join node.
+    Node findJoinNode() {
+      // Finding the join node
+      Node u = graph.source(in_edge);
+      Node v = graph.target(in_edge);
+      while (u != v) {
+	if (depth[u] == depth[v]) {
+	  u = parent[u];
+	  v = parent[v];
+	}
+	else if (depth[u] > depth[v]) u = parent[u];
+	else v = parent[v];
+      }
+      return u;
+    }
+
+    /// \brief Finds the leaving edge of the cycle. Returns \c true if
+    /// the leaving edge is not the same as the entering edge.
+    bool findLeavingEdge() {
+      // Initializing first and second nodes according to the direction
+      // of the cycle
+      if (state[in_edge] == LOWER) {
+	first = graph.source(in_edge);
+	second	= graph.target(in_edge);
+      } else {
+	first = graph.target(in_edge);
+	second	= graph.source(in_edge);
+      }
+      delta = capacity[in_edge];
+      bool result = false;
+      Capacity d;
+      Edge e;
+
+      // Searching the cycle along the path form the first node to the
+      // root node
+      for (Node u = first; u != join; u = parent[u]) {
+	e = pred_edge[u];
+	d = forward[u] ? flow[e] : capacity[e] - flow[e];
+	if (d < delta) {
+	  delta = d;
+	  u_out = u;
+	  u_in = first;
+	  v_in = second;
+	  result = true;
+	}
+      }
+      // Searching the cycle along the path form the second node to the
+      // root node
+      for (Node u = second; u != join; u = parent[u]) {
+	e = pred_edge[u];
+	d = forward[u] ? capacity[e] - flow[e] : flow[e];
+	if (d <= delta) {
+	  delta = d;
+	  u_out = u;
+	  u_in = second;
+	  v_in = first;
+	  result = true;
+	}
+      }
+      return result;
+    }
+
+    /// \brief Changes flow and state edge maps.
+    void changeFlows(bool change) {
+      // Augmenting along the cycle
+      if (delta > 0) {
+	Capacity val = state[in_edge] * delta;
+	flow[in_edge] += val;
+	for (Node u = graph.source(in_edge); u != join; u = parent[u]) {
+	  flow[pred_edge[u]] += forward[u] ? -val : val;
+	}
+	for (Node u = graph.target(in_edge); u != join; u = parent[u]) {
+	  flow[pred_edge[u]] += forward[u] ? val : -val;
+	}
+      }
+      // Updating the state of the entering and leaving edges
+      if (change) {
+	state[in_edge] = TREE;
+	state[pred_edge[u_out]] =
+	  (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
+      } else {
+	state[in_edge] = -state[in_edge];
+      }
+    }
+
+    /// \brief Updates thread and parent node maps.
+    void updateThreadParent() {
+      Node u;
+      v_out = parent[u_out];
+
+      // Handling the case when join and v_out coincide
+      bool par_first = false;
+      if (join == v_out) {
+	for (u = join; u != u_in && u != v_in; u = thread[u]) ;
+	if (u == v_in) {
+	  par_first = true;
+	  while (thread[u] != u_out) u = thread[u];
+	  first = u;
+	}
+      }
+
+      // Finding the last successor of u_in (u) and the node after it
+      // (right) according to the thread index
+      for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;
+      right = thread[u];
+      if (thread[v_in] == u_out) {
+	for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
+	if (last == u_out) last = thread[last];
+      }
+      else last = thread[v_in];
+
+      // Updating stem nodes
+      thread[v_in] = stem = u_in;
+      par_stem = v_in;
+      while (stem != u_out) {
+	thread[u] = new_stem = parent[stem];
+
+	// Finding the node just before the stem node (u) according to
+	// the original thread index
+	for (u = new_stem; thread[u] != stem; u = thread[u]) ;
+	thread[u] = right;
+
+	// Changing the parent node of stem and shifting stem and
+	// par_stem nodes
+	parent[stem] = par_stem;
+	par_stem = stem;
+	stem = new_stem;
+
+	// Finding the last successor of stem (u) and the node after it
+	// (right) according to the thread index
+	for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
+	right = thread[u];
+      }
+      parent[u_out] = par_stem;
+      thread[u] = last;
+
+      if (join == v_out && par_first) {
+	if (first != v_in) thread[first] = right;
+      } else {
+	for (u = v_out; thread[u] != u_out; u = thread[u]) ;
+	thread[u] = right;
+      }
+    }
+
+    /// \brief Updates pred_edge and forward node maps.
+    void updatePredEdge() {
+      Node u = u_out, v;
+      while (u != u_in) {
+	v = parent[u];
+	pred_edge[u] = pred_edge[v];
+	forward[u] = !forward[v];
+	u = v;
+      }
+      pred_edge[u_in] = in_edge;
+      forward[u_in] = (u_in == graph.source(in_edge));
+    }
+
+    /// \brief Updates depth and potential node maps.
+    void updateDepthPotential() {
+      depth[u_in] = depth[v_in] + 1;
+      potential[u_in] = forward[u_in] ?
+	potential[v_in] + cost[pred_edge[u_in]] :
+	potential[v_in] - cost[pred_edge[u_in]];
+
+      Node u = thread[u_in], v;
+      while (true) {
+	v = parent[u];
+	if (v == INVALID) break;
+	depth[u] = depth[v] + 1;
+	potential[u] = forward[u] ?
+	  potential[v] + cost[pred_edge[u]] :
+	  potential[v] - cost[pred_edge[u]];
+	if (depth[u] <= depth[v_in]) break;
+	u = thread[u];
+      }
+    }
+
+    /// \brief Executes the algorithm.
+    bool start() {
+      // Processing pivots
+#ifdef _DEBUG_ITER_
+      int iter_num = 0;
+#endif
+#if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT)
+      EdgeIt next_edge(graph);
+      while (findEnteringEdge(next_edge))
+#else
+      while (findEnteringEdge())
+#endif
+      {
+	join = findJoinNode();
+	bool change = findLeavingEdge();
+	changeFlows(change);
+	if (change) {
+	  updateThreadParent();
+	  updatePredEdge();
+	  updateDepthPotential();
+	}
+#ifdef _DEBUG_ITER_
+	++iter_num;
+#endif
+      }
+
+#ifdef _DEBUG_ITER_
+      t_iter.stop();
+      std::cout << "Network Simplex algorithm finished. " << iter_num
+		<< " pivot iterations performed.";
+#endif
+
+      // Checking if the flow amount equals zero on all the
+      // artificial edges
+      for (InEdgeIt e(graph, root); e != INVALID; ++e)
+	if (flow[e] > 0) return false;
+      for (OutEdgeIt e(graph, root); e != INVALID; ++e)
+	if (flow[e] > 0) return false;
+
+      // Copying flow values to flow_result
+      if (lower) {
+	for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
+	  flow_result[e] = (*lower)[e] + flow[edge_ref[e]];
+      } else {
+	for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
+	  flow_result[e] = flow[edge_ref[e]];
+      }
+      // Copying potential values to potential_result
+      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
+	potential_result[n] = potential[node_ref[n]];
+
+      return true;
+    }
+
+  }; //class NetworkSimplex
+
+  ///@}
+
+} //namespace lemon
+
+#endif //LEMON_NETWORK_SIMPLEX_H