LEMON/LEMON-official Changeset - r656:e6927fe719e6

Repositories » LEMON » LEMON-official

Summary
Changelog
Files
Switch To
- loading...
Options
- Lightweight changelog
- Search
Pull Requests

Changeset - r656:e6927fe719e6

« 655:6ac5d9

657:dacc2c »

r656:e6927fe719e6

Fri, 17 Apr 2009 18:04:36

[Not Reviewed]

0 comments (0 inline)

kpeter (Peter Kovacs)
kpeter@inf.elte.hu

Support >= and <= constraints in NetworkSimplex (#219, #234) By default the same inequality constraints are supported as by Circulation (the GEQ form), but the LEQ form can also be selected using the problemType() function. The documentation of the min. cost flow module is reworked and extended with important notes and explanations about the different variants of the problem and about the dual solution and optimality conditions.

0 3 0

default

3 files changed with 377 insertions and 93 deletions:

doc/groups.dox

lemon/network_simplex.h

206

test/min_cost_flow_test.cc

↑ Collapse diff ↑

doc/groups.dox

Show inline comments

@@ -128,441 +128,509 @@
 		\endcode
 		This means that in a situation, when a <tt>const %ListDigraph&</tt>
 		reference to a graph is given, then it have to be instantiated with
 		<tt>Digraph=const %ListDigraph</tt>.
 		\code
 		int algorithm1(const ListDigraph& g) {
 		  ReverseDigraph<const ListDigraph> rg(g);
 		  return algorithm2(rg);
+		}
 		\endcode
 		*/
 		/**
 		@defgroup semi_adaptors Semi-Adaptor Classes for Graphs
 		@ingroup graphs
 		\brief Graph types between real graphs and graph adaptors.
 		This group describes some graph types between real graphs and graph adaptors.
 		These classes wrap graphs to give new functionality as the adaptors do it.
 		On the other hand they are not light-weight structures as the adaptors.
 		*/
 		/**
 		@defgroup maps Maps
 		@ingroup datas
 		\brief Map structures implemented in LEMON.
 		This group describes the map structures implemented in LEMON.
 		LEMON provides several special purpose maps and map adaptors that e.g. combine
 		new maps from existing ones.
 		<b>See also:</b> \ref map_concepts "Map Concepts".
 		*/
 		/**
 		@defgroup graph_maps Graph Maps
 		@ingroup maps
 		\brief Special graph-related maps.
 		This group describes maps that are specifically designed to assign
 		values to the nodes and arcs/edges of graphs.
 		If you are looking for the standard graph maps (\c NodeMap, \c ArcMap,
 		\c EdgeMap), see the \ref graph_concepts "Graph Structure Concepts".
 		*/
 		/**
 		\defgroup map_adaptors Map Adaptors
 		\ingroup maps
 		\brief Tools to create new maps from existing ones
 		This group describes map adaptors that are used to create "implicit"
 		maps from other maps.
 		Most of them are \ref concepts::ReadMap "read-only maps".
 		They can make arithmetic and logical operations between one or two maps
 		(negation, shifting, addition, multiplication, logical 'and', 'or',
 		'not' etc.) or e.g. convert a map to another one of different Value type.
 		The typical usage of this classes is passing implicit maps to
 		algorithms.  If a function type algorithm is called then the function
 		type map adaptors can be used comfortable. For example let's see the
 		usage of map adaptors with the \c graphToEps() function.
 		\code
 		  Color nodeColor(int deg) {
 		    if (deg >= 2) {
 		      return Color(0.5, 0.0, 0.5);
 		    } else if (deg == 1) {
 		      return Color(1.0, 0.5, 1.0);
 		    } else {
 		      return Color(0.0, 0.0, 0.0);
+		    }
+		  }
 		  Digraph::NodeMap<int> degree_map(graph);
 		  graphToEps(graph, "graph.eps")
 		    .coords(coords).scaleToA4().undirected()
 		    .nodeColors(composeMap(functorToMap(nodeColor), degree_map))
 		    .run();
 		\endcode
 		The \c functorToMap() function makes an \c int to \c Color map from the
 		\c nodeColor() function. The \c composeMap() compose the \c degree_map
 		and the previously created map. The composed map is a proper function to
 		get the color of each node.
 		The usage with class type algorithms is little bit harder. In this
 		case the function type map adaptors can not be used, because the
 		function map adaptors give back temporary objects.
 		\code
 		  Digraph graph;
 		  typedef Digraph::ArcMap<double> DoubleArcMap;
 		  DoubleArcMap length(graph);
 		  DoubleArcMap speed(graph);
 		  typedef DivMap<DoubleArcMap, DoubleArcMap> TimeMap;
 		  TimeMap time(length, speed);
 		  Dijkstra<Digraph, TimeMap> dijkstra(graph, time);
 		  dijkstra.run(source, target);
 		\endcode
 		We have a length map and a maximum speed map on the arcs of a digraph.
 		The minimum time to pass the arc can be calculated as the division of
 		the two maps which can be done implicitly with the \c DivMap template
 		class. We use the implicit minimum time map as the length map of the
 		\c Dijkstra algorithm.
 		*/
 		/**
 		@defgroup matrices Matrices
 		@ingroup datas
 		\brief Two dimensional data storages implemented in LEMON.
 		This group describes two dimensional data storages implemented in LEMON.
 		*/
 		/**
 		@defgroup paths Path Structures
 		@ingroup datas
 		\brief %Path structures implemented in LEMON.
 		This group describes the path structures implemented in LEMON.
 		LEMON provides flexible data structures to work with paths.
 		All of them have similar interfaces and they can be copied easily with
 		assignment operators and copy constructors. This makes it easy and
 		efficient to have e.g. the Dijkstra algorithm to store its result in
 		any kind of path structure.
 		\sa lemon::concepts::Path
 		*/
 		/**
 		@defgroup auxdat Auxiliary Data Structures
 		@ingroup datas
 		\brief Auxiliary data structures implemented in LEMON.
 		This group describes some data structures implemented in LEMON in
 		order to make it easier to implement combinatorial algorithms.
 		*/
 		/**
 		@defgroup algs Algorithms
 		\brief This group describes the several algorithms
 		implemented in LEMON.
 		This group describes the several algorithms
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup search Graph Search
 		@ingroup algs
 		\brief Common graph search algorithms.
 		This group describes the common graph search algorithms, namely
 		\e breadth-first \e search (BFS) and \e depth-first \e search (DFS).
 		*/
 		/**
 		@defgroup shortest_path Shortest Path Algorithms
 		@ingroup algs
 		\brief Algorithms for finding shortest paths.
 		This group describes the algorithms for finding shortest paths in digraphs.
 		 - \ref Dijkstra algorithm for finding shortest paths from a source node
 		   when all arc lengths are non-negative.
 		 - \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths
 		   from a source node when arc lenghts can be either positive or negative,
 		   but the digraph should not contain directed cycles with negative total
 		   length.
 		 - \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms
 		   for solving the \e all-pairs \e shortest \e paths \e problem when arc
 		   lenghts can be either positive or negative, but the digraph should
 		   not contain directed cycles with negative total length.
 		 - \ref Suurballe A successive shortest path algorithm for finding
 		   arc-disjoint paths between two nodes having minimum total length.
 		*/
 		/**
 		@defgroup max_flow Maximum Flow Algorithms
 		@ingroup algs
 		\brief Algorithms for finding maximum flows.
 		This group describes the algorithms for finding maximum flows and
 		feasible circulations.
 		The \e maximum \e flow \e problem is to find a flow of maximum value between
 		a single source and a single target. Formally, there is a \f$G=(V,A)\f$
 		digraph, a \f$cap:A\rightarrow\mathbf{R}^+_0\f$ capacity function and
+		digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
 		\f$s, t \in V\f$ source and target nodes.
 		A maximum flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of the
+		A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
 		following optimization problem.
 		\f[ \max\sum_{a\in\delta_{out}(s)}f(a) - \sum_{a\in\delta_{in}(s)}f(a) \f]
 		\f[ \sum_{a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)
 		    \qquad \forall v\in V\setminus\{s,t\} \f]
 		\f[ 0 \leq f(a) \leq cap(a) \qquad \forall a\in A \f]
 		\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
 		\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
 		    \quad \forall u\in V\setminus\{s,t\} \f]
 		\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
 		LEMON contains several algorithms for solving maximum flow problems:
 		- \ref EdmondsKarp Edmonds-Karp algorithm.
 		- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm.
 		- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees.
 		- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees.
 		In most cases the \ref Preflow "Preflow" algorithm provides the
 		fastest method for computing a maximum flow. All implementations
 		provides functions to also query the minimum cut, which is the dual
 		problem of the maximum flow.
 		*/
 		/**
 		@defgroup min_cost_flow Minimum Cost Flow Algorithms
 		@ingroup algs
 		\brief Algorithms for finding minimum cost flows and circulations.
-		This group describes the algorithms for finding minimum cost flows and
+		This group contains the algorithms for finding minimum cost flows and
 		circulations.
 		The \e minimum \e cost \e flow \e problem is to find a feasible flow of
 		minimum total cost from a set of supply nodes to a set of demand nodes
-		in a network with capacity constraints and arc costs.
+		in a network with capacity constraints (lower and upper bounds)
 		and arc costs.
 		Formally, let \f$G=(V,A)\f$ be a digraph,
 		\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
 		upper bounds for the flow values on the arcs,
+		upper bounds for the flow values on the arcs, for which
 		\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
 		\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
 		on the arcs, and
 		\f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values
 		of the nodes.
 		A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of
-		the following optimization problem.
+		on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
 		signed supply values of the nodes.
 		If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
 		supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
 		\f$-sup(u)\f$ demand.
 		A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
 		of the following optimization problem.
 		\f[ \min\sum_{a\in A} f(a) cost(a) \f]
 		\f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a) =
 		    supply(v) \qquad \forall v\in V \f]
 		\f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]
 		\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
 		\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
 		    sup(u) \quad \forall u\in V \f]
 		\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
 		LEMON contains several algorithms for solving minimum cost flow problems:
 		 - \ref CycleCanceling Cycle-canceling algorithms.
-		 - \ref CapacityScaling Successive shortest path algorithm with optional
+		The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
 		zero or negative in order to have a feasible solution (since the sum
 		of the expressions on the left-hand side of the inequalities is zero).
 		It means that the total demand must be greater or equal to the total
 		supply and all the supplies have to be carried out from the supply nodes,
 		but there could be demands that are not satisfied.
 		If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
 		constraints have to be satisfied with equality, i.e. all demands
 		have to be satisfied and all supplies have to be used.
 		If you need the opposite inequalities in the supply/demand constraints
 		(i.e. the total demand is less than the total supply and all the demands
 		have to be satisfied while there could be supplies that are not used),
 		then you could easily transform the problem to the above form by reversing
 		the direction of the arcs and taking the negative of the supply values
 		(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
 		However \ref NetworkSimplex algorithm also supports this form directly
 		for the sake of convenience.
 		A feasible solution for this problem can be found using \ref Circulation.
 		Note that the above formulation is actually more general than the usual
 		definition of the minimum cost flow problem, in which strict equalities
 		are required in the supply/demand contraints, i.e.
 		\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
 		    sup(u) \quad \forall u\in V. \f]
 		However if the sum of the supply values is zero, then these two problems
 		are equivalent. So if you need the equality form, you have to ensure this
 		additional contraint for the algorithms.
 		The dual solution of the minimum cost flow problem is represented by node
 		potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
 		An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
 		is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
 		node potentials the following \e complementary \e slackness optimality
 		conditions hold.
 		 - For all \f$uv\in A\f$ arcs:
 		   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
 		   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
 		   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
 		 - For all \f$u\in V\f$:
 		   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
 		     then \f$\pi(u)=0\f$.
 		Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
 		\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
 		\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
 		All algorithms provide dual solution (node potentials) as well
 		if an optimal flow is found.
 		LEMON contains several algorithms for solving minimum cost flow problems.
 		 - \ref NetworkSimplex Primal Network Simplex algorithm with various
 		   pivot strategies.
 		 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
 		   cost scaling.
 		 - \ref CapacityScaling Successive Shortest %Path algorithm with optional
 		   capacity scaling.
 		 - \ref CostScaling Push-relabel and augment-relabel algorithms based on
 		   cost scaling.
 		 - \ref NetworkSimplex Primal network simplex algorithm with various
 		   pivot strategies.
 		 - \ref CancelAndTighten The Cancel and Tighten algorithm.
 		 - \ref CycleCanceling Cycle-Canceling algorithms.
 		Most of these implementations support the general inequality form of the
 		minimum cost flow problem, but CancelAndTighten and CycleCanceling
 		only support the equality form due to the primal method they use.
 		In general NetworkSimplex is the most efficient implementation,
 		but in special cases other algorithms could be faster.
 		For example, if the total supply and/or capacities are rather small,
 		CapacityScaling is usually the fastest algorithm (without effective scaling).
 		*/
 		/**
 		@defgroup min_cut Minimum Cut Algorithms
 		@ingroup algs
 		\brief Algorithms for finding minimum cut in graphs.
 		This group describes the algorithms for finding minimum cut in graphs.
 		The \e minimum \e cut \e problem is to find a non-empty and non-complete
 		\f$X\f$ subset of the nodes with minimum overall capacity on
 		outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a
 		\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum
 		cut is the \f$X\f$ solution of the next optimization problem:
 		\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}}
 		    \sum_{uv\in A, u\in X, v\not\in X}cap(uv) \f]
 		LEMON contains several algorithms related to minimum cut problems:
 		- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut
 		  in directed graphs.
 		- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for
 		  calculating minimum cut in undirected graphs.
 		- \ref GomoryHuTree "Gomory-Hu tree computation" for calculating
 		  all-pairs minimum cut in undirected graphs.
 		If you want to find minimum cut just between two distinict nodes,
 		see the \ref max_flow "maximum flow problem".
 		*/
 		/**
 		@defgroup graph_prop Connectivity and Other Graph Properties
 		@ingroup algs
 		\brief Algorithms for discovering the graph properties
 		This group describes the algorithms for discovering the graph properties
 		like connectivity, bipartiteness, euler property, simplicity etc.
 		\image html edge_biconnected_components.png
 		\image latex edge_biconnected_components.eps "bi-edge-connected components" width=\textwidth
 		*/
 		/**
 		@defgroup planar Planarity Embedding and Drawing
 		@ingroup algs
 		\brief Algorithms for planarity checking, embedding and drawing
 		This group describes the algorithms for planarity checking,
 		embedding and drawing.
 		\image html planar.png
 		\image latex planar.eps "Plane graph" width=\textwidth
 		*/
 		/**
 		@defgroup matching Matching Algorithms
 		@ingroup algs
 		\brief Algorithms for finding matchings in graphs and bipartite graphs.
 		This group contains algorithm objects and functions to calculate
 		matchings in graphs and bipartite graphs. The general matching problem is
 		finding a subset of the arcs which does not shares common endpoints.
 		There are several different algorithms for calculate matchings in
 		graphs.  The matching problems in bipartite graphs are generally
 		easier than in general graphs. The goal of the matching optimization
 		can be finding maximum cardinality, maximum weight or minimum cost
 		matching. The search can be constrained to find perfect or
 		maximum cardinality matching.
 		The matching algorithms implemented in LEMON:
 		- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm
 		  for calculating maximum cardinality matching in bipartite graphs.
 		- \ref PrBipartiteMatching Push-relabel algorithm
 		  for calculating maximum cardinality matching in bipartite graphs.
 		- \ref MaxWeightedBipartiteMatching
 		  Successive shortest path algorithm for calculating maximum weighted
 		  matching and maximum weighted bipartite matching in bipartite graphs.
 		- \ref MinCostMaxBipartiteMatching
 		  Successive shortest path algorithm for calculating minimum cost maximum
 		  matching in bipartite graphs.
 		- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating
 		  maximum cardinality matching in general graphs.
 		- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating
 		  maximum weighted matching in general graphs.
 		- \ref MaxWeightedPerfectMatching
 		  Edmond's blossom shrinking algorithm for calculating maximum weighted
 		  perfect matching in general graphs.
 		\image html bipartite_matching.png
 		\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth
 		*/
 		/**
 		@defgroup spantree Minimum Spanning Tree Algorithms
 		@ingroup algs
 		\brief Algorithms for finding a minimum cost spanning tree in a graph.
 		This group describes the algorithms for finding a minimum cost spanning
 		tree in a graph.
 		*/
 		/**
 		@defgroup auxalg Auxiliary Algorithms
 		@ingroup algs
 		\brief Auxiliary algorithms implemented in LEMON.
 		This group describes some algorithms implemented in LEMON
 		in order to make it easier to implement complex algorithms.
 		*/
 		/**
 		@defgroup approx Approximation Algorithms
 		@ingroup algs
 		\brief Approximation algorithms.
 		This group describes the approximation and heuristic algorithms
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup gen_opt_group General Optimization Tools
 		\brief This group describes some general optimization frameworks
 		implemented in LEMON.
 		This group describes some general optimization frameworks
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup lp_group Lp and Mip Solvers
 		@ingroup gen_opt_group
 		\brief Lp and Mip solver interfaces for LEMON.
 		This group describes Lp and Mip solver interfaces for LEMON. The
 		various LP solvers could be used in the same manner with this
 		interface.
 		*/
 		/**
 		@defgroup lp_utils Tools for Lp and Mip Solvers
 		@ingroup lp_group
 		\brief Helper tools to the Lp and Mip solvers.
 		This group adds some helper tools to general optimization framework
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup metah Metaheuristics
 		@ingroup gen_opt_group
 		\brief Metaheuristics for LEMON library.
 		This group describes some metaheuristic optimization tools.
 		*/
 		/**
 		@defgroup utils Tools and Utilities
 		\brief Tools and utilities for programming in LEMON
 		Tools and utilities for programming in LEMON.
 		*/
 		/**
 		@defgroup gutils Basic Graph Utilities
 		@ingroup utils
 		\brief Simple basic graph utilities.
 		This group describes some simple basic graph utilities.
 		*/
 		/**
 		@defgroup misc Miscellaneous Tools
 		@ingroup utils
 		\brief Tools for development, debugging and testing.
 		This group describes several useful tools for development,
 		debugging and testing.
 		*/
 		/**
 		@defgroup timecount Time Measuring and Counting
 		@ingroup misc
 		\brief Simple tools for measuring the performance of algorithms.
 		This group describes simple tools for measuring the performance
 		of algorithms.
 		*/
 		/**

lemon/network_simplex.h

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * 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 Network Simplex algorithm for finding a minimum cost flow.
 		#include <vector>
 		#include <limits>
 		#include <algorithm>
 		#include <lemon/core.h>
 		#include <lemon/math.h>
 		#include <lemon/maps.h>
 		#include <lemon/circulation.h>
 		#include <lemon/adaptors.h>
 		namespace lemon {
 		  /// \addtogroup min_cost_flow
 		  /// @{
 		  /// \brief Implementation of the primal Network Simplex algorithm
 		  /// for finding a \ref min_cost_flow "minimum cost flow".
 		  ///
 		  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
 		  /// for finding a \ref min_cost_flow "minimum cost flow".
 		  /// 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 this class is the fastest implementation available
 		  /// in LEMON for the minimum cost flow problem.
 		  /// Moreover it supports both direction of the supply/demand inequality
 		  /// constraints. For more information see \ref ProblemType.
 		  ///
 		  /// \tparam GR The digraph type the algorithm runs on.
 		  /// \tparam F The value type used for flow amounts, capacity bounds
 		  /// and supply values in the algorithm. By default it is \c int.
 		  /// \tparam C The value type used for costs and potentials in the
 		  /// algorithm. By default it is the same as \c F.
 		  ///
 		  /// \warning Both value types must be signed and all input data must
 		  /// be integer.
 		  ///
 		  /// \note %NetworkSimplex provides five different pivot rule
-		  /// implementations. For more information see \ref PivotRule.
+		  /// implementations, from which the most efficient one is used
 		  /// by default. For more information see \ref PivotRule.
 		  template <typename GR, typename F = int, typename C = F>
 		  class NetworkSimplex
+		  {
 		  public:
 		    /// The flow type of the algorithm
 		    typedef F Flow;
 		    /// The cost type of the algorithm
 		    typedef C Cost;
 		#ifdef DOXYGEN
 		    /// The type of the flow map
 		    typedef GR::ArcMap<Flow> FlowMap;
 		    /// The type of the potential map
 		    typedef GR::NodeMap<Cost> PotentialMap;
 		#else
 		    /// The type of the flow map
 		    typedef typename GR::template ArcMap<Flow> FlowMap;
 		    /// The type of the potential map
 		    typedef typename GR::template NodeMap<Cost> PotentialMap;
 		#endif
 		  public:
 		    /// \brief Enum type for selecting the pivot rule.
 		    ///
 		    /// Enum type for selecting the pivot rule for the \ref run()
 		    /// function.
 		    ///
 		    /// \ref NetworkSimplex provides five different pivot rule
 		    /// implementations that significantly affect the running time
 		    /// of the algorithm.
 		    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
 		    /// proved to be the most efficient and the most robust on various
 		    /// test inputs according to our benchmark tests.
 		    /// However another pivot rule can be selected using the \ref run()
 		    /// function with the proper parameter.
 		    enum PivotRule {
 		      /// The First Eligible pivot rule.
 		      /// The next eligible arc is selected in a wraparound fashion
 		      /// in every iteration.
 		      FIRST_ELIGIBLE,
 		      /// The Best Eligible pivot rule.
 		      /// The best eligible arc is selected in every iteration.
 		      BEST_ELIGIBLE,
 		      /// The Block Search pivot rule.
 		      /// A specified number of arcs are examined in every iteration
 		      /// in a wraparound fashion and the best eligible arc is selected
 		      /// from this block.
 		      BLOCK_SEARCH,
 		      /// The Candidate List pivot rule.
 		      /// In a major iteration a candidate list is built from eligible arcs
 		      /// in a wraparound fashion and in the following minor iterations
 		      /// the best eligible arc is selected from this list.
 		      CANDIDATE_LIST,
 		      /// The Altering Candidate List pivot rule.
 		      /// It is a modified version of the Candidate List method.
 		      /// It keeps only the several best eligible arcs from the former
 		      /// candidate list and extends this list in every iteration.
 		      ALTERING_LIST
 		    };
 		    /// \brief Enum type for selecting the problem type.
 		    ///
 		    /// Enum type for selecting the problem type, i.e. the direction of
 		    /// the inequalities in the supply/demand constraints of the
 		    /// \ref min_cost_flow "minimum cost flow problem".
 		    ///
 		    /// The default problem type is \c GEQ, since this form is supported
 		    /// by other minimum cost flow algorithms and the \ref Circulation
 		    /// algorithm as well.
 		    /// The \c LEQ problem type can be selected using the \ref problemType()
 		    /// function.
 		    ///
 		    /// Note that the equality form is a special case of both problem type.
 		    enum ProblemType {
 		      /// This option means that there are "<em>greater or equal</em>"
 		      /// constraints in the defintion, i.e. the exact formulation of the
 		      /// problem is the following.
 		      /**
 		          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
 		          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
 		              sup(u) \quad \forall u\in V \f]
 		          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
 		      */
 		      /// It means that the total demand must be greater or equal to the
 		      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
 		      /// negative) and all the supplies have to be carried out from
 		      /// the supply nodes, but there could be demands that are not
 		      /// satisfied.
 		      GEQ,
 		      /// It is just an alias for the \c GEQ option.
 		      CARRY_SUPPLIES = GEQ,
 		      /// This option means that there are "<em>less or equal</em>"
 		      /// constraints in the defintion, i.e. the exact formulation of the
 		      /// problem is the following.
 		      /**
 		          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
 		          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
 		              sup(u) \quad \forall u\in V \f]
 		          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
 		      */
 		      /// It means that the total demand must be less or equal to the
 		      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
 		      /// positive) and all the demands have to be satisfied, but there
 		      /// could be supplies that are not carried out from the supply
 		      /// nodes.
 		      LEQ,
 		      /// It is just an alias for the \c LEQ option.
 		      SATISFY_DEMANDS = LEQ
 		    };
 		  private:
 		    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		    typedef typename GR::template ArcMap<Flow> FlowArcMap;
 		    typedef typename GR::template ArcMap<Cost> CostArcMap;
 		    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
 		    typedef std::vector<Arc> ArcVector;
 		    typedef std::vector<Node> NodeVector;
 		    typedef std::vector<int> IntVector;
 		    typedef std::vector<bool> BoolVector;
 		    typedef std::vector<Flow> FlowVector;
 		    typedef std::vector<Cost> CostVector;
 		    // State constants for arcs
 		    enum ArcStateEnum {
 		      STATE_UPPER = -1,
 		      STATE_TREE  =  0,
 		      STATE_LOWER =  1
 		    };
 		  private:
 		    // Data related to the underlying digraph
 		    const GR &_graph;
 		    int _node_num;
 		    int _arc_num;
 		    // Parameters of the problem
 		    FlowArcMap *_plower;
 		    FlowArcMap *_pupper;
 		    CostArcMap *_pcost;
 		    FlowNodeMap *_psupply;
 		    bool _pstsup;
 		    Node _psource, _ptarget;
 		    Flow _pstflow;
 		    ProblemType _ptype;
 		    // Result maps
 		    FlowMap *_flow_map;
 		    PotentialMap *_potential_map;
 		    bool _local_flow;
 		    bool _local_potential;
 		    // Data structures for storing the digraph
 		    IntNodeMap _node_id;
 		    ArcVector _arc_ref;
 		    IntVector _source;
 		    IntVector _target;
 		    // Node and arc data
 		    FlowVector _cap;
 		    CostVector _cost;
 		    FlowVector _supply;
 		    FlowVector _flow;
 		    CostVector _pi;
 		    // Data for storing the spanning tree structure
 		    IntVector _parent;
 		    IntVector _pred;
 		    IntVector _thread;
 		    IntVector _rev_thread;
 		    IntVector _succ_num;
 		    IntVector _last_succ;
 		    IntVector _dirty_revs;
 		    BoolVector _forward;
 		    IntVector _state;
 		    int _root;
 		    // Temporary data used in the current pivot iteration
 		    int in_arc, join, u_in, v_in, u_out, v_out;
 		    int first, second, right, last;
 		    int stem, par_stem, new_stem;
 		    Flow delta;
 		  private:
 		    // Implementation of the First Eligible pivot rule
 		    class FirstEligiblePivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
 		      const CostVector &_cost;
 		      const IntVector  &_state;
 		      const CostVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		      // Pivot rule data
 		      int _next_arc;
 		    public:
 		      // Constructor
 		      FirstEligiblePivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
 		      {}
 		      // Find next entering arc
 		      bool findEnteringArc() {
 		        Cost c;
 		        for (int e = _next_arc; e < _arc_num; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < 0) {
 		            _in_arc = e;
 		            _next_arc = e + 1;
 		            return true;
+		          }
+		        }
 		        for (int e = 0; e < _next_arc; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < 0) {
 		            _in_arc = e;
 		            _next_arc = e + 1;
 		            return true;
+		          }
+		        }
 		        return false;
+		      }
 		    }; //class FirstEligiblePivotRule
 		    // Implementation of the Best Eligible pivot rule
 		    class BestEligiblePivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
 		      const CostVector &_cost;
 		      const IntVector  &_state;
 		      const CostVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		    public:
 		      // Constructor
 		      BestEligiblePivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num)
 		      {}
 		      // Find next entering arc
 		      bool findEnteringArc() {
 		        Cost c, min = 0;
 		        for (int e = 0; e < _arc_num; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < min) {
 		            min = c;
 		            _in_arc = e;
+		          }
+		        }
 		        return min < 0;
+		      }
 		    }; //class BestEligiblePivotRule
 		    // Implementation of the Block Search pivot rule
 		    class BlockSearchPivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
 		      const CostVector &_cost;
 		      const IntVector  &_state;
 		      const CostVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		      // Pivot rule data
 		      int _block_size;
 		      int _next_arc;
 		    public:
 		      // Constructor
 		      BlockSearchPivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+		      {
 		        // The main parameters of the pivot rule
 		        const double BLOCK_SIZE_FACTOR = 2.0;
 		        const int MIN_BLOCK_SIZE = 10;
 		        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
 		                                MIN_BLOCK_SIZE );
+		      }
 		      // Find next entering arc
 		      bool findEnteringArc() {
 		        Cost c, min = 0;
 		        int cnt = _block_size;
 		        int e, min_arc = _next_arc;
 		        for (e = _next_arc; e < _arc_num; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < min) {
 		            min = c;
 		            min_arc = e;
+		          }
 		          if (--cnt == 0) {
 		            if (min < 0) break;
 		            cnt = _block_size;
+		          }
+		        }
 		        if (min == 0 || cnt > 0) {
 		          for (e = 0; e < _next_arc; ++e) {
 		            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		            if (c < min) {
 		              min = c;
 		              min_arc = e;
+		            }
 		            if (--cnt == 0) {
 		              if (min < 0) break;
 		              cnt = _block_size;
+		            }
+		          }
@@ -397,1019 +463,1120 @@
 		                                 MIN_MINOR_LIMIT );
 		        _curr_length = _minor_count = 0;
 		        _candidates.resize(_list_length);
+		      }
 		      /// Find next entering arc
 		      bool findEnteringArc() {
 		        Cost min, c;
 		        int e, min_arc = _next_arc;
 		        if (_curr_length > 0 && _minor_count < _minor_limit) {
 		          // Minor iteration: select the best eligible arc from the
 		          // current candidate list
 		          ++_minor_count;
 		          min = 0;
 		          for (int i = 0; i < _curr_length; ++i) {
 		            e = _candidates[i];
 		            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		            if (c < min) {
 		              min = c;
 		              min_arc = e;
+		            }
 		            if (c >= 0) {
 		              _candidates[i--] = _candidates[--_curr_length];
+		            }
+		          }
 		          if (min < 0) {
 		            _in_arc = min_arc;
 		            return true;
+		          }
+		        }
 		        // Major iteration: build a new candidate list
 		        min = 0;
 		        _curr_length = 0;
 		        for (e = _next_arc; e < _arc_num; ++e) {
 		          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (c < 0) {
 		            _candidates[_curr_length++] = e;
 		            if (c < min) {
 		              min = c;
 		              min_arc = e;
+		            }
 		            if (_curr_length == _list_length) break;
+		          }
+		        }
 		        if (_curr_length < _list_length) {
 		          for (e = 0; e < _next_arc; ++e) {
 		            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		            if (c < 0) {
 		              _candidates[_curr_length++] = e;
 		              if (c < min) {
 		                min = c;
 		                min_arc = e;
+		              }
 		              if (_curr_length == _list_length) break;
+		            }
+		          }
+		        }
 		        if (_curr_length == 0) return false;
 		        _minor_count = 1;
 		        _in_arc = min_arc;
 		        _next_arc = e;
 		        return true;
+		      }
 		    }; //class CandidateListPivotRule
 		    // Implementation of the Altering Candidate List pivot rule
 		    class AlteringListPivotRule
+		    {
 		    private:
 		      // References to the NetworkSimplex class
 		      const IntVector  &_source;
 		      const IntVector  &_target;
 		      const CostVector &_cost;
 		      const IntVector  &_state;
 		      const CostVector &_pi;
 		      int &_in_arc;
 		      int _arc_num;
 		      // Pivot rule data
 		      int _block_size, _head_length, _curr_length;
 		      int _next_arc;
 		      IntVector _candidates;
 		      CostVector _cand_cost;
 		      // Functor class to compare arcs during sort of the candidate list
 		      class SortFunc
+		      {
 		      private:
 		        const CostVector &_map;
 		      public:
 		        SortFunc(const CostVector &map) : _map(map) {}
 		        bool operator()(int left, int right) {
 		          return _map[left] > _map[right];
+		        }
 		      };
 		      SortFunc _sort_func;
 		    public:
 		      // Constructor
 		      AlteringListPivotRule(NetworkSimplex &ns) :
 		        _source(ns._source), _target(ns._target),
 		        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
 		        _in_arc(ns.in_arc), _arc_num(ns._arc_num),
 		        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
+		      {
 		        // The main parameters of the pivot rule
 		        const double BLOCK_SIZE_FACTOR = 1.5;
 		        const int MIN_BLOCK_SIZE = 10;
 		        const double HEAD_LENGTH_FACTOR = 0.1;
 		        const int MIN_HEAD_LENGTH = 3;
 		        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
 		                                MIN_BLOCK_SIZE );
 		        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
 		                                 MIN_HEAD_LENGTH );
 		        _candidates.resize(_head_length + _block_size);
 		        _curr_length = 0;
+		      }
 		      // Find next entering arc
 		      bool findEnteringArc() {
 		        // Check the current candidate list
 		        int e;
 		        for (int i = 0; i < _curr_length; ++i) {
 		          e = _candidates[i];
 		          _cand_cost[e] = _state[e] *
 		            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (_cand_cost[e] >= 0) {
 		            _candidates[i--] = _candidates[--_curr_length];
+		          }
+		        }
 		        // Extend the list
 		        int cnt = _block_size;
 		        int last_arc = 0;
 		        int limit = _head_length;
 		        for (int e = _next_arc; e < _arc_num; ++e) {
 		          _cand_cost[e] = _state[e] *
 		            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		          if (_cand_cost[e] < 0) {
 		            _candidates[_curr_length++] = e;
 		            last_arc = e;
+		          }
 		          if (--cnt == 0) {
 		            if (_curr_length > limit) break;
 		            limit = 0;
 		            cnt = _block_size;
+		          }
+		        }
 		        if (_curr_length <= limit) {
 		          for (int e = 0; e < _next_arc; ++e) {
 		            _cand_cost[e] = _state[e] *
 		              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
 		            if (_cand_cost[e] < 0) {
 		              _candidates[_curr_length++] = e;
 		              last_arc = e;
+		            }
 		            if (--cnt == 0) {
 		              if (_curr_length > limit) break;
 		              limit = 0;
 		              cnt = _block_size;
+		            }
+		          }
+		        }
 		        if (_curr_length == 0) return false;
 		        _next_arc = last_arc + 1;
 		        // Make heap of the candidate list (approximating a partial sort)
 		        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
 		                   _sort_func );
 		        // Pop the first element of the heap
 		        _in_arc = _candidates[0];
 		        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
 		                  _sort_func );
 		        _curr_length = std::min(_head_length, _curr_length - 1);
 		        return true;
+		      }
 		    }; //class AlteringListPivotRule
 		  public:
 		    /// \brief Constructor.
 		    ///
-		    /// Constructor.
+		    /// The constructor of the class.
 		    ///
 		    /// \param graph The digraph the algorithm runs on.
 		    NetworkSimplex(const GR& graph) :
 		      _graph(graph),
 		      _plower(NULL), _pupper(NULL), _pcost(NULL),
 		      _psupply(NULL), _pstsup(false),
+		      _psupply(NULL), _pstsup(false), _ptype(GEQ),
 		      _flow_map(NULL), _potential_map(NULL),
 		      _local_flow(false), _local_potential(false),
 		      _node_id(graph)
+		    {
 		      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
 		                   std::numeric_limits<Flow>::is_signed,
 		        "The flow type of NetworkSimplex must be signed integer");
 		      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
 		                   std::numeric_limits<Cost>::is_signed,
 		        "The cost type of NetworkSimplex must be signed integer");
+		    }
 		    /// Destructor.
 		    ~NetworkSimplex() {
 		      if (_local_flow) delete _flow_map;
 		      if (_local_potential) delete _potential_map;
+		    }
 		    /// \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 neither this function nor \ref boundMaps() is 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 Flow type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template <typename LOWER>
 		    NetworkSimplex& lowerMap(const LOWER& map) {
 		      delete _plower;
 		      _plower = new FlowArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        (*_plower)[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 none of the functions \ref upperMap(), \ref capacityMap()
 		    /// and \ref boundMaps() is used before calling \ref run(),
 		    /// the upper bounds (capacities) will be set to
 		    /// \c std::numeric_limits<Flow>::max() on all arcs.
 		    ///
 		    /// \param map An arc map storing the upper bounds.
 		    /// Its \c Value type must be convertible to the \c Flow type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename UPPER>
 		    NetworkSimplex& upperMap(const UPPER& map) {
 		      delete _pupper;
 		      _pupper = new FlowArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        (*_pupper)[a] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the upper bounds (capacities) on the arcs.
 		    ///
 		    /// This function sets the upper bounds (capacities) on the arcs.
 		    /// It is just an alias for \ref upperMap().
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename CAP>
 		    NetworkSimplex& capacityMap(const CAP& map) {
 		      return upperMap(map);
+		    }
 		    /// \brief Set the lower and upper bounds on the arcs.
 		    ///
 		    /// This function sets the lower and upper bounds on the arcs.
 		    /// If neither this function nor \ref lowerMap() is used before
 		    /// calling \ref run(), the lower bounds will be set to zero
 		    /// on all arcs.
 		    /// If none of the functions \ref upperMap(), \ref capacityMap()
 		    /// and \ref boundMaps() is used before calling \ref run(),
 		    /// the upper bounds (capacities) will be set to
 		    /// \c std::numeric_limits<Flow>::max() on all arcs.
 		    ///
 		    /// \param lower An arc map storing the lower bounds.
 		    /// \param upper An arc map storing the upper bounds.
 		    ///
 		    /// The \c Value type of the maps must be convertible to the
 		    /// \c Flow type of the algorithm.
 		    ///
 		    /// \note This function is just a shortcut of calling \ref lowerMap()
 		    /// and \ref upperMap() separately.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template <typename LOWER, typename UPPER>
 		    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
 		      return lowerMap(lower).upperMap(upper);
+		    }
 		    /// \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 COST>
 		    NetworkSimplex& costMap(const COST& map) {
 		      delete _pcost;
 		      _pcost = new CostArcMap(_graph);
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        (*_pcost)[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.
 		    /// (It makes sense only if non-zero lower bounds are given.)
 		    ///
 		    /// \param map A node map storing the supply values.
 		    /// Its \c Value type must be convertible to the \c Flow type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename SUP>
 		    NetworkSimplex& supplyMap(const SUP& map) {
 		      delete _psupply;
 		      _pstsup = false;
 		      _psupply = new FlowNodeMap(_graph);
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        (*_psupply)[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.
 		    /// (It makes sense only if non-zero lower bounds are given.)
 		    ///
 		    /// \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>
 		    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
 		      delete _psupply;
 		      _psupply = NULL;
 		      _pstsup = true;
 		      _psource = s;
 		      _ptarget = t;
 		      _pstflow = k;
 		      return *this;
+		    }
 		    /// \brief Set the problem type.
 		    ///
 		    /// This function sets the problem type for the algorithm.
 		    /// If it is not used before calling \ref run(), the \ref GEQ problem
 		    /// type will be used.
 		    ///
 		    /// For more information see \ref ProblemType.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& problemType(ProblemType problem_type) {
 		      _ptype = problem_type;
 		      return *this;
+		    }
 		    /// \brief Set the flow map.
 		    ///
 		    /// This function sets the flow map.
 		    /// If it is not used before calling \ref run(), an instance will
 		    /// be allocated automatically. The destructor deallocates this
 		    /// automatically allocated map, of course.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& flowMap(FlowMap& map) {
 		      if (_local_flow) {
 		        delete _flow_map;
 		        _local_flow = false;
+		      }
 		      _flow_map = &map;
 		      return *this;
+		    }
 		    /// \brief Set the potential map.
 		    ///
 		    /// This function sets the potential map, which is used for storing
 		    /// the dual solution.
 		    /// If it is not used before calling \ref run(), an instance will
 		    /// be allocated automatically. The destructor deallocates this
 		    /// automatically allocated map, of course.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& potentialMap(PotentialMap& map) {
 		      if (_local_potential) {
 		        delete _potential_map;
 		        _local_potential = false;
+		      }
 		      _potential_map = &map;
 		      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 \ref lowerMap(),
+		    /// The paramters can be specified using functions \ref lowerMap(),
 		    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
 		    /// \ref costMap(), \ref supplyMap() and \ref stSupply()
 		    /// functions. For example,
 		    /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
 		    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
 		    /// For example,
 		    /// \code
 		    ///   NetworkSimplex<ListDigraph> ns(graph);
 		    ///   ns.boundMaps(lower, upper).costMap(cost)
 		    ///     .supplyMap(sup).run();
 		    /// \endcode
 		    ///
 		    /// 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.
 		    ///
 		    /// \param pivot_rule The pivot rule that will be used during the
 		    /// algorithm. For more information see \ref PivotRule.
 		    ///
 		    /// \return \c true if a feasible flow can be found.
 		    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
 		      return init() && start(pivot_rule);
+		    }
 		    /// \brief Reset all the parameters that have been given before.
 		    ///
 		    /// This function resets all the paramaters that have been given
 		    /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
 		    /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
-		    /// \ref stSupply() functions before.
+		    /// before using functions \ref lowerMap(), \ref upperMap(),
 		    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
 		    /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
 		    /// \ref flowMap() and \ref potentialMap().
 		    ///
 		    /// 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.
 		    ///
 		    /// For example,
 		    /// \code
 		    ///   NetworkSimplex<ListDigraph> ns(graph);
 		    ///
 		    ///   // First run
 		    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
 		    ///     .supplyMap(sup).run();
 		    ///
 		    ///   // 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 reset()
 		    ///   // (the lower bounds will be set to zero on all arcs)
 		    ///   ns.reset();
 		    ///   ns.capacityMap(cap).costMap(cost)
 		    ///     .supplyMap(sup).run();
 		    /// \endcode
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& reset() {
 		      delete _plower;
 		      delete _pupper;
 		      delete _pcost;
 		      delete _psupply;
 		      _plower = NULL;
 		      _pupper = NULL;
 		      _pcost = NULL;
 		      _psupply = NULL;
 		      _pstsup = false;
 		      _ptype = GEQ;
 		      if (_local_flow) delete _flow_map;
 		      if (_local_potential) delete _potential_map;
 		      _flow_map = NULL;
 		      _potential_map = NULL;
 		      _local_flow = false;
 		      _local_potential = false;
 		      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.
 		    /// The complexity of the function is O(e).
 		    ///
 		    /// \note The return type of the function can be specified as a
 		    /// template parameter. For example,
 		    /// \code
 		    ///   ns.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 Num>
 		    Num totalCost() const {
 		      Num c = 0;
 		      if (_pcost) {
 		        for (ArcIt e(_graph); e != INVALID; ++e)
 		          c += (*_flow_map)[e] * (*_pcost)[e];
 		      } else {
 		        for (ArcIt e(_graph); e != INVALID; ++e)
 		          c += (*_flow_map)[e];
+		      }
 		      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.
 		    Flow flow(const Arc& a) const {
 		      return (*_flow_map)[a];
+		    }
 		    /// \brief Return a const reference to the flow map.
 		    ///
 		    /// This function returns a const reference to an arc map storing
 		    /// the found flow.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    const FlowMap& flowMap() const {
 		      return *_flow_map;
+		    }
 		    /// \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 (*_potential_map)[n];
+		    }
 		    /// \brief Return a const reference to the potential map
 		    /// (the dual solution).
 		    ///
 		    /// This function returns a const reference to a node map storing
 		    /// the found potentials, which form the dual solution of the
 		    /// \ref min_cost_flow "minimum cost flow" problem.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    const PotentialMap& potentialMap() const {
 		      return *_potential_map;
+		    }
 		    /// @}
 		  private:
 		    // Initialize internal data structures
 		    bool init() {
 		      // Initialize result maps
 		      if (!_flow_map) {
 		        _flow_map = new FlowMap(_graph);
 		        _local_flow = true;
+		      }
 		      if (!_potential_map) {
 		        _potential_map = new PotentialMap(_graph);
 		        _local_potential = true;
+		      }
 		      // Initialize vectors
 		      _node_num = countNodes(_graph);
 		      _arc_num = countArcs(_graph);
 		      int all_node_num = _node_num + 1;
 		      int all_arc_num = _arc_num + _node_num;
 		      if (_node_num == 0) return false;
 		      _arc_ref.resize(_arc_num);
 		      _source.resize(all_arc_num);
 		      _target.resize(all_arc_num);
 		      _cap.resize(all_arc_num);
 		      _cost.resize(all_arc_num);
 		      _supply.resize(all_node_num);
 		      _flow.resize(all_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(all_arc_num);
 		      // Initialize node related data
 		      bool valid_supply = true;
 		      Flow sum_supply = 0;
 		      if (!_pstsup && !_psupply) {
 		        _pstsup = true;
 		        _psource = _ptarget = NodeIt(_graph);
 		        _pstflow = 0;
+		      }
 		      if (_psupply) {
 		        Flow sum = 0;
 		        int i = 0;
 		        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
 		          _node_id[n] = i;
 		          _supply[i] = (*_psupply)[n];
-		          sum += _supply[i];
+		          sum_supply += _supply[i];
+		        }
-		        valid_supply = (sum == 0);
+		        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
 		                       (_ptype == LEQ && sum_supply >= 0);
 		      } else {
 		        int i = 0;
 		        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
 		          _node_id[n] = i;
 		          _supply[i] = 0;
+		        }
 		        _supply[_node_id[_psource]] =  _pstflow;
-		        _supply[_node_id[_ptarget]]   = -_pstflow;
 		        _supply[_node_id[_ptarget]] = -_pstflow;
+		      }
 		      if (!valid_supply) return false;
 		      // Infinite capacity value
 		      Flow inf_cap =
 		        std::numeric_limits<Flow>::has_infinity ?
 		        std::numeric_limits<Flow>::infinity() :
 		        std::numeric_limits<Flow>::max();
 		      // Initialize artifical cost
 		      Cost art_cost;
 		      if (std::numeric_limits<Cost>::is_exact) {
 		        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
 		      } else {
 		        art_cost = std::numeric_limits<Cost>::min();
 		        for (int i = 0; i != _arc_num; ++i) {
 		          if (_cost[i] > art_cost) art_cost = _cost[i];
+		        }
 		        art_cost = (art_cost + 1) * _node_num;
+		      }
 		      // Run Circulation to check if a feasible solution exists
 		      typedef ConstMap<Arc, Flow> ConstArcMap;
 		      FlowNodeMap *csup = NULL;
 		      bool local_csup = false;
 		      if (_psupply) {
 		        csup = _psupply;
 		      } else {
 		        csup = new FlowNodeMap(_graph, 0);
 		        (*csup)[_psource] =  _pstflow;
 		        (*csup)[_ptarget] = -_pstflow;
 		        local_csup = true;
+		      }
 		      bool circ_result = false;
 		      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
 		        // GEQ problem type
 		        if (_plower) {
 		          if (_pupper) {
 		            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
 		              circ(_graph, *_plower, *_pupper, *csup);
 		            circ_result = circ.run();
 		          } else {
 		            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
 		              circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
 		            circ_result = circ.run();
+		          }
 		        } else {
 		          if (_pupper) {
 		            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
 		              circ(_graph, ConstArcMap(0), *_pupper, *csup);
 		            circ_result = circ.run();
 		          } else {
 		            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
 		              circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
 		            circ_result = circ.run();
+		          }
+		        }
 		      } else {
 		        // LEQ problem type
 		        typedef ReverseDigraph<const GR> RevGraph;
 		        typedef NegMap<FlowNodeMap> NegNodeMap;
 		        RevGraph rgraph(_graph);
 		        NegNodeMap neg_csup(*csup);
 		        if (_plower) {
 		          if (_pupper) {
 		            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
 		              circ(rgraph, *_plower, *_pupper, neg_csup);
 		            circ_result = circ.run();
 		          } else {
 		            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
 		              circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
 		            circ_result = circ.run();
+		          }
 		        } else {
 		          if (_pupper) {
 		            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
 		              circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
 		            circ_result = circ.run();
 		          } else {
 		            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
 		              circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
 		            circ_result = circ.run();
+		          }
+		        }
+		      }
 		      if (local_csup) delete csup;
 		      if (!circ_result) return false;
 		      // Set data for the artificial root node
 		      _root = _node_num;
 		      _parent[_root] = -1;
 		      _pred[_root] = -1;
 		      _thread[_root] = 0;
 		      _rev_thread[0] = _root;
 		      _succ_num[_root] = all_node_num;
 		      _last_succ[_root] = _root - 1;
 		      _supply[_root] = 0;
 		      _pi[_root] = 0;
 		      _supply[_root] = -sum_supply;
 		      if (sum_supply < 0) {
 		        _pi[_root] = -art_cost;
 		      } else {
 		        _pi[_root] = art_cost;
+		      }
 		      // Store the arcs in a mixed order
 		      int k = std::max(int(sqrt(_arc_num)), 10);
 		      int i = 0;
 		      for (ArcIt e(_graph); e != INVALID; ++e) {
 		        _arc_ref[i] = e;
 		        if ((i += k) >= _arc_num) i = (i % k) + 1;
+		      }
 		      // Initialize arc maps
 		      Flow inf_cap =
 		        std::numeric_limits<Flow>::has_infinity ?
 		        std::numeric_limits<Flow>::infinity() :
 		        std::numeric_limits<Flow>::max();
 		      if (_pupper && _pcost) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
 		          _source[i] = _node_id[_graph.source(e)];
 		          _target[i] = _node_id[_graph.target(e)];
 		          _cap[i] = (*_pupper)[e];
 		          _cost[i] = (*_pcost)[e];
 		          _flow[i] = 0;
 		          _state[i] = STATE_LOWER;
+		        }
 		      } else {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
 		          _source[i] = _node_id[_graph.source(e)];
 		          _target[i] = _node_id[_graph.target(e)];
 		          _flow[i] = 0;
 		          _state[i] = STATE_LOWER;
+		        }
 		        if (_pupper) {
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cap[i] = (*_pupper)[_arc_ref[i]];
 		        } else {
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cap[i] = inf_cap;
+		        }
 		        if (_pcost) {
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cost[i] = (*_pcost)[_arc_ref[i]];
 		        } else {
 		          for (int i = 0; i != _arc_num; ++i)
 		            _cost[i] = 1;
+		        }
+		      }
 		      // Initialize artifical cost
 		      Cost art_cost;
 		      if (std::numeric_limits<Cost>::is_exact) {
 		        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
 		      } else {
 		        art_cost = std::numeric_limits<Cost>::min();
 		        for (int i = 0; i != _arc_num; ++i) {
 		          if (_cost[i] > art_cost) art_cost = _cost[i];
+		        }
 		        art_cost = (art_cost + 1) * _node_num;
+		      }
 		      // Remove non-zero lower bounds
 		      if (_plower) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Flow c = (*_plower)[_arc_ref[i]];
 		          if (c != 0) {
 		            _cap[i] -= c;
 		            _supply[_source[i]] -= c;
 		            _supply[_target[i]] += c;
+		          }
+		        }
+		      }
 		      // Add artificial arcs and initialize the spanning tree data structure
 		      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		        _thread[u] = u + 1;
 		        _rev_thread[u + 1] = u;
 		        _succ_num[u] = 1;
 		        _last_succ[u] = u;
 		        _parent[u] = _root;
 		        _pred[u] = e;
 		        _cost[e] = art_cost;
 		        _cap[e] = inf_cap;
 		        _state[e] = STATE_TREE;
 		        if (_supply[u] >= 0) {
+		        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
 		          _flow[e] = _supply[u];
 		          _forward[u] = true;
 		          _pi[u] = -art_cost;
+		          _pi[u] = -art_cost + _pi[_root];
 		        } else {
 		          _flow[e] = -_supply[u];
 		          _forward[u] = false;
 		          _pi[u] = art_cost;
+		          _pi[u] = art_cost + _pi[_root];
+		        }
+		      }
 		      return true;
+		    }
 		    // Find the join node
 		    void findJoinNode() {
 		      int u = _source[in_arc];
 		      int v = _target[in_arc];
 		      while (u != v) {
 		        if (_succ_num[u] < _succ_num[v]) {
 		          u = _parent[u];
 		        } else {
 		          v = _parent[v];
+		        }
+		      }
 		      join = u;
+		    }
 		    // Find the leaving arc of the cycle and returns true if the
 		    // leaving arc is not the same as the entering arc
 		    bool findLeavingArc() {
 		      // Initialize first and second nodes according to the direction
 		      // of the cycle
 		      if (_state[in_arc] == STATE_LOWER) {
 		        first  = _source[in_arc];
 		        second = _target[in_arc];
 		      } else {
 		        first  = _target[in_arc];
 		        second = _source[in_arc];
+		      }
 		      delta = _cap[in_arc];
 		      int result = 0;
 		      Flow d;
 		      int e;
 		      // Search the cycle along the path form the first node to the root
 		      for (int u = first; u != join; u = _parent[u]) {
 		        e = _pred[u];
 		        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
 		        if (d < delta) {
 		          delta = d;
 		          u_out = u;
 		          result = 1;
+		        }
+		      }
 		      // Search the cycle along the path form the second node to the root
 		      for (int u = second; u != join; u = _parent[u]) {
 		        e = _pred[u];
 		        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
 		        if (d <= delta) {
 		          delta = d;
 		          u_out = u;
 		          result = 2;
+		        }
+		      }
 		      if (result == 1) {
 		        u_in = first;
 		        v_in = second;
 		      } else {
 		        u_in = second;
 		        v_in = first;
+		      }
 		      return result != 0;
+		    }
 		    // Change _flow and _state vectors
 		    void changeFlow(bool change) {
 		      // Augment along the cycle
 		      if (delta > 0) {
 		        Flow val = _state[in_arc] * delta;
 		        _flow[in_arc] += val;
 		        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
 		          _flow[_pred[u]] += _forward[u] ? -val : val;
+		        }
 		        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
 		          _flow[_pred[u]] += _forward[u] ? val : -val;
+		        }
+		      }
 		      // Update the state of the entering and leaving arcs
 		      if (change) {
 		        _state[in_arc] = STATE_TREE;
 		        _state[_pred[u_out]] =
 		          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
 		      } else {
 		        _state[in_arc] = -_state[in_arc];
+		      }
+		    }
 		    // Update the tree structure
 		    void updateTreeStructure() {
 		      int u, w;
 		      int old_rev_thread = _rev_thread[u_out];
 		      int old_succ_num = _succ_num[u_out];
 		      int old_last_succ = _last_succ[u_out];
 		      v_out = _parent[u_out];
 		      u = _last_succ[u_in];  // the last successor of u_in
 		      right = _thread[u];    // the node after it
 		      // Handle the case when old_rev_thread equals to v_in
 		      // (it also means that join and v_out coincide)
 		      if (old_rev_thread == v_in) {
 		        last = _thread[_last_succ[u_out]];
 		      } else {
 		        last = _thread[v_in];
+		      }
 		      // Update _thread and _parent along the stem nodes (i.e. the nodes
 		      // between u_in and u_out, whose parent have to be changed)
 		      _thread[v_in] = stem = u_in;
 		      _dirty_revs.clear();
 		      _dirty_revs.push_back(v_in);
 		      par_stem = v_in;
 		      while (stem != u_out) {
 		        // Insert the next stem node into the thread list
 		        new_stem = _parent[stem];
 		        _thread[u] = new_stem;
 		        _dirty_revs.push_back(u);
 		        // Remove the subtree of stem from the thread list
 		        w = _rev_thread[stem];
 		        _thread[w] = right;
 		        _rev_thread[right] = w;
 		        // Change the parent node and shift stem nodes
 		        _parent[stem] = par_stem;
 		        par_stem = stem;
 		        stem = new_stem;
 		        // Update u and right
 		        u = _last_succ[stem] == _last_succ[par_stem] ?
 		          _rev_thread[par_stem] : _last_succ[stem];
 		        right = _thread[u];
+		      }
 		      _parent[u_out] = par_stem;
 		      _thread[u] = last;
 		      _rev_thread[last] = u;
 		      _last_succ[u_out] = u;
 		      // Remove the subtree of u_out from the thread list except for
 		      // the case when old_rev_thread equals to v_in
 		      // (it also means that join and v_out coincide)
 		      if (old_rev_thread != v_in) {
 		        _thread[old_rev_thread] = right;
 		        _rev_thread[right] = old_rev_thread;
+		      }
 		      // Update _rev_thread using the new _thread values
 		      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
 		        u = _dirty_revs[i];
 		        _rev_thread[_thread[u]] = u;
+		      }
 		      // Update _pred, _forward, _last_succ and _succ_num for the
 		      // stem nodes from u_out to u_in
 		      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
 		      u = u_out;
 		      while (u != u_in) {
 		        w = _parent[u];
 		        _pred[u] = _pred[w];
 		        _forward[u] = !_forward[w];
 		        tmp_sc += _succ_num[u] - _succ_num[w];
 		        _succ_num[u] = tmp_sc;
 		        _last_succ[w] = tmp_ls;
 		        u = w;
+		      }
 		      _pred[u_in] = in_arc;
 		      _forward[u_in] = (u_in == _source[in_arc]);
 		      _succ_num[u_in] = old_succ_num;
 		      // Set limits for updating _last_succ form v_in and v_out
 		      // towards the root
 		      int up_limit_in = -1;
 		      int up_limit_out = -1;
 		      if (_last_succ[join] == v_in) {
 		        up_limit_out = join;
 		      } else {
 		        up_limit_in = join;
+		      }
 		      // Update _last_succ from v_in towards the root
 		      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
 		           u = _parent[u]) {
 		        _last_succ[u] = _last_succ[u_out];
+		      }
 		      // Update _last_succ from v_out towards the root
 		      if (join != old_rev_thread && v_in != old_rev_thread) {
 		        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
 		             u = _parent[u]) {
 		          _last_succ[u] = old_rev_thread;
+		        }
 		      } else {
 		        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
 		             u = _parent[u]) {
 		          _last_succ[u] = _last_succ[u_out];
+		        }
+		      }
 		      // Update _succ_num from v_in to join
 		      for (u = v_in; u != join; u = _parent[u]) {
 		        _succ_num[u] += old_succ_num;
+		      }
 		      // Update _succ_num from v_out to join
 		      for (u = v_out; u != join; u = _parent[u]) {
 		        _succ_num[u] -= old_succ_num;
+		      }
+		    }
 		    // Update potentials
 		    void updatePotential() {
 		      Cost sigma = _forward[u_in] ?
 		        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
 		        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
 		      // Update potentials in the subtree, which has been moved
 		      int end = _thread[_last_succ[u_in]];
 		      for (int u = u_in; u != end; u = _thread[u]) {
 		        _pi[u] += sigma;
+		      }
+		    }
 		    // Execute the algorithm
 		    bool start(PivotRule pivot_rule) {
 		      // Select the pivot rule implementation
 		      switch (pivot_rule) {
 		        case FIRST_ELIGIBLE:
 		          return start<FirstEligiblePivotRule>();
 		        case BEST_ELIGIBLE:
 		          return start<BestEligiblePivotRule>();
 		        case BLOCK_SEARCH:
 		          return start<BlockSearchPivotRule>();
 		        case CANDIDATE_LIST:
 		          return start<CandidateListPivotRule>();
 		        case ALTERING_LIST:
 		          return start<AlteringListPivotRule>();
+		      }
 		      return false;
+		    }
 		    template <typename PivotRuleImpl>
 		    bool start() {
 		      PivotRuleImpl pivot(*this);
 		      // Execute the Network Simplex algorithm
 		      while (pivot.findEnteringArc()) {
 		        findJoinNode();
 		        bool change = findLeavingArc();
 		        changeFlow(change);
 		        if (change) {
 		          updateTreeStructure();
 		          updatePotential();
+		        }
+		      }
 		      // Check if the flow amount equals zero on all the artificial arcs
 		      for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
 		        if (_flow[e] > 0) return false;
+		      }
 		      // Copy flow values to _flow_map
 		      if (_plower) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Arc e = _arc_ref[i];
 		          _flow_map->set(e, (*_plower)[e] + _flow[i]);
+		        }
 		      } else {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          _flow_map->set(_arc_ref[i], _flow[i]);
+		        }
+		      }
 		      // Copy potential values to _potential_map
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        _potential_map->set(n, _pi[_node_id[n]]);
+		      }
 		      return true;
+		    }
 		  }; //class NetworkSimplex
 		  ///@}
 		} //namespace lemon
 		#endif //LEMON_NETWORK_SIMPLEX_H

test/min_cost_flow_test.cc

Show inline comments

 		/* -*- mode: C++; indent-tabs-mode: nil; -*-
+		 *
 		 * This file is a part of LEMON, a generic C++ optimization library.
+		 *
 		 * Copyright (C) 2003-2009
 		 * 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.
+		 *
 		 */
 		#include <iostream>
 		#include <fstream>
 		#include <lemon/list_graph.h>
 		#include <lemon/lgf_reader.h>
 		#include <lemon/network_simplex.h>
 		#include <lemon/concepts/digraph.h>
 		#include <lemon/concept_check.h>
 		#include "test_tools.h"
 		using namespace lemon;
 		char test_lgf[] =
 		  "@nodes\n"
 		  "label  sup1 sup2 sup3\n"
 		  "    1    20   27    0\n"
 		  "    2    -4    0    0\n"
 		  "    3     0    0    0\n"
 		  "    4     0    0    0\n"
 		  "    5     9    0    0\n"
 		  "    6    -6    0    0\n"
 		  "    7     0    0    0\n"
 		  "    8     0    0    0\n"
 		  "    9     3    0    0\n"
 		  "   10    -2    0    0\n"
 		  "   11     0    0    0\n"
-		  "   12   -20  -27    0\n"
+		  "label  sup1 sup2 sup3 sup4 sup5\n"
 		  "    1    20   27    0   20   30\n"
 		  "    2    -4    0    0   -8   -3\n"
 		  "    3     0    0    0    0    0\n"
 		  "    4     0    0    0    0    0\n"
 		  "    5     9    0    0    6   11\n"
 		  "    6    -6    0    0   -5   -6\n"
 		  "    7     0    0    0    0    0\n"
 		  "    8     0    0    0    0    3\n"
 		  "    9     3    0    0    0    0\n"
 		  "   10    -2    0    0   -7   -2\n"
 		  "   11     0    0    0  -10    0\n"
 		  "   12   -20  -27    0  -30  -20\n"
 		  "\n"
 		  "@arcs\n"
 		  "       cost  cap low1 low2\n"
 		  " 1  2    70   11    0    8\n"
 		  " 1  3   150    3    0    1\n"
 		  " 1  4    80   15    0    2\n"
 		  " 2  8    80   12    0    0\n"
 		  " 3  5   140    5    0    3\n"
 		  " 4  6    60   10    0    1\n"
 		  " 4  7    80    2    0    0\n"
 		  " 4  8   110    3    0    0\n"
 		  " 5  7    60   14    0    0\n"
 		  " 5 11   120   12    0    0\n"
 		  " 6  3     0    3    0    0\n"
 		  " 6  9   140    4    0    0\n"
 		  " 6 10    90    8    0    0\n"
 		  " 7  1    30    5    0    0\n"
 		  " 8 12    60   16    0    4\n"
 		  " 9 12    50    6    0    0\n"
 		  "10 12    70   13    0    5\n"
 		  "10  2   100    7    0    0\n"
 		  "10  7    60   10    0    0\n"
 		  "11 10    20   14    0    6\n"
 		  "12 11    30   10    0    0\n"
 		  "\n"
 		  "@attributes\n"
 		  "source 1\n"
 		  "target 12\n";
 		enum ProblemType {
 		  EQ,
 		  GEQ,
 		  LEQ
 		};
 		// Check the interface of an MCF algorithm
 		template <typename GR, typename Flow, typename Cost>
 		class McfClassConcept
+		{
 		public:
 		  template <typename MCF>
 		  struct Constraints {
 		    void constraints() {
 		      checkConcept<concepts::Digraph, GR>();
 		      MCF mcf(g);
 		      b = mcf.reset()
 		             .lowerMap(lower)
 		             .upperMap(upper)
 		             .capacityMap(upper)
 		             .boundMaps(lower, upper)
 		             .costMap(cost)
 		             .supplyMap(sup)
 		             .stSupply(n, n, k)
 		             .flowMap(flow)
 		             .potentialMap(pot)
 		             .run();
 		      const MCF& const_mcf = mcf;
 		      const typename MCF::FlowMap &fm = mcf.flowMap();
 		      const typename MCF::PotentialMap &pm = mcf.potentialMap();
 		      const typename MCF::FlowMap &fm = const_mcf.flowMap();
 		      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
 		      v = mcf.totalCost();
 		      double x = mcf.template totalCost<double>();
 		      v = mcf.flow(a);
 		      v = mcf.potential(n);
 		      mcf.flowMap(flow);
 		      mcf.potentialMap(pot);
 		      v = const_mcf.totalCost();
 		      double x = const_mcf.template totalCost<double>();
 		      v = const_mcf.flow(a);
 		      v = const_mcf.potential(n);
 		      ignore_unused_variable_warning(fm);
 		      ignore_unused_variable_warning(pm);
 		      ignore_unused_variable_warning(x);
+		    }
 		    typedef typename GR::Node Node;
 		    typedef typename GR::Arc Arc;
 		    typedef concepts::ReadMap<Node, Flow> NM;
 		    typedef concepts::ReadMap<Arc, Flow> FAM;
 		    typedef concepts::ReadMap<Arc, Cost> CAM;
 		    const GR &g;
 		    const FAM &lower;
 		    const FAM &upper;
 		    const CAM &cost;
 		    const NM &sup;
 		    const Node &n;
 		    const Arc &a;
 		    const Flow &k;
 		    Flow v;
 		    bool b;
 		    typename MCF::FlowMap &flow;
 		    typename MCF::PotentialMap &pot;
 		  };
 		};
 		// Check the feasibility of the given flow (primal soluiton)
 		template < typename GR, typename LM, typename UM,
 		           typename SM, typename FM >
 		bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
-		                const SM& supply, const FM& flow )
+		                const SM& supply, const FM& flow,
 		                ProblemType type = EQ )
+		{
 		  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		  for (ArcIt e(gr); e != INVALID; ++e) {
 		    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
+		  }
 		  for (NodeIt n(gr); n != INVALID; ++n) {
 		    typename SM::Value sum = 0;
 		    for (OutArcIt e(gr, n); e != INVALID; ++e)
 		      sum += flow[e];
 		    for (InArcIt e(gr, n); e != INVALID; ++e)
 		      sum -= flow[e];
-		    if (sum != supply[n]) return false;
+		    bool b = (type ==  EQ && sum == supply[n]) ||
 		             (type == GEQ && sum >= supply[n]) ||
 		             (type == LEQ && sum <= supply[n]);
 		    if (!b) return false;
+		  }
 		  return true;
+		}
 		// Check the feasibility of the given potentials (dual soluiton)
 		// using the "Complementary Slackness" optimality condition
 		template < typename GR, typename LM, typename UM,
 		           typename CM, typename FM, typename PM >
+		           typename CM, typename SM, typename FM, typename PM >
 		bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
-		                     const CM& cost, const FM& flow, const PM& pi )
+		                     const CM& cost, const SM& supply, const FM& flow,
 		                     const PM& pi )
+		{
 		  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		  bool opt = true;
 		  for (ArcIt e(gr); opt && e != INVALID; ++e) {
 		    typename CM::Value red_cost =
 		      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
 		    opt = red_cost == 0 ||
 		          (red_cost > 0 && flow[e] == lower[e]) ||
 		          (red_cost < 0 && flow[e] == upper[e]);
+		  }
 		  for (NodeIt n(gr); opt && n != INVALID; ++n) {
 		    typename SM::Value sum = 0;
 		    for (OutArcIt e(gr, n); e != INVALID; ++e)
 		      sum += flow[e];
 		    for (InArcIt e(gr, n); e != INVALID; ++e)
 		      sum -= flow[e];
 		    opt = (sum == supply[n]) || (pi[n] == 0);
+		  }
 		  return opt;
+		}
 		// Run a minimum cost flow algorithm and check the results
 		template < typename MCF, typename GR,
 		           typename LM, typename UM,
 		           typename CM, typename SM >
 		void checkMcf( const MCF& mcf, bool mcf_result,
 		               const GR& gr, const LM& lower, const UM& upper,
 		               const CM& cost, const SM& supply,
 		               bool result, typename CM::Value total,
-		               const std::string &test_id = "" )
+		               const std::string &test_id = "",
 		               ProblemType type = EQ )
+		{
 		  check(mcf_result == result, "Wrong result " + test_id);
 		  if (result) {
 		    check(checkFlow(gr, lower, upper, supply, mcf.flowMap()),
+		    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
 		          "The flow is not feasible " + test_id);
 		    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
 		    check(checkPotential(gr, lower, upper, cost, mcf.flowMap(),
+		    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
 		                         mcf.potentialMap()),
 		          "Wrong potentials " + test_id);
+		  }
+		}
 		int main()
+		{
 		  // Check the interfaces
+		  {
 		    typedef int Flow;
 		    typedef int Cost;
 		    // TODO: This typedef should be enabled if the standard maps are
 		    // reference maps in the graph concepts (See #190).
 		/**/
 		    //typedef concepts::Digraph GR;
 		    typedef ListDigraph GR;
 		/**/
 		    checkConcept< McfClassConcept<GR, Flow, Cost>,
 		                  NetworkSimplex<GR, Flow, Cost> >();
+		  }
 		  // 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), u(gr);
 		  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr);
+		  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
 		  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
 		  Node v, w;
 		  std::istringstream input(test_lgf);
 		  DigraphReader<Digraph>(gr, input)
 		    .arcMap("cost", c)
 		    .arcMap("cap", u)
 		    .arcMap("low1", l1)
 		    .arcMap("low2", l2)
 		    .nodeMap("sup1", s1)
 		    .nodeMap("sup2", s2)
 		    .nodeMap("sup3", s3)
 		    .nodeMap("sup4", s4)
 		    .nodeMap("sup5", s5)
 		    .node("source", v)
 		    .node("target", w)
 		    .run();
 		  // A. Test NetworkSimplex with the default pivot rule
+		  {
 		    NetworkSimplex<Digraph> mcf(gr);
 		    // Check the equality form
 		    mcf.upperMap(u).costMap(c);
 		    checkMcf(mcf, mcf.supplyMap(s1).run(),
 		             gr, l1, u, c, s1, true,  5240, "#A1");
 		    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
 		             gr, l1, u, c, s2, true,  7620, "#A2");
 		    mcf.lowerMap(l2);
 		    checkMcf(mcf, mcf.supplyMap(s1).run(),
 		             gr, l2, u, c, s1, true,  5970, "#A3");
 		    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
 		             gr, l2, u, c, s2, true,  8010, "#A4");
 		    mcf.reset();
 		    checkMcf(mcf, mcf.supplyMap(s1).run(),
 		             gr, l1, cu, cc, s1, true,  74, "#A5");
 		    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
 		             gr, l2, cu, cc, s2, true,  94, "#A6");
 		    mcf.reset();
 		    checkMcf(mcf, mcf.run(),
 		             gr, l1, cu, cc, s3, true,   0, "#A7");
 		    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
 		             gr, l2, u, cc, s3, false,   0, "#A8");
 		    // Check the GEQ form
 		    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
 		    checkMcf(mcf, mcf.run(),
 		             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
 		    mcf.problemType(mcf.GEQ);
 		    checkMcf(mcf, mcf.lowerMap(l2).run(),
 		             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
 		    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
 		    checkMcf(mcf, mcf.run(),
 		             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
 		    // Check the LEQ form
 		    mcf.reset().problemType(mcf.LEQ);
 		    mcf.upperMap(u).costMap(c).supplyMap(s5);
 		    checkMcf(mcf, mcf.run(),
 		             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
 		    checkMcf(mcf, mcf.lowerMap(l2).run(),
 		             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
 		    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
 		    checkMcf(mcf, mcf.run(),
 		             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
+		  }
 		  // B. Test NetworkSimplex with each pivot rule
+		  {
 		    NetworkSimplex<Digraph> mcf(gr);
 		    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
 		             gr, l2, u, c, s1, true,  5970, "#B1");
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
 		             gr, l2, u, c, s1, true,  5970, "#B2");
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
 		             gr, l2, u, c, s1, true,  5970, "#B3");
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
 		             gr, l2, u, c, s1, true,  5970, "#B4");
 		    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
 		             gr, l2, u, c, s1, true,  5970, "#B5");
+		  }
 		  return 0;
+		}

0 comments (0 inline)

RhodeCode

Login to your account