LEMON/LEMON-official Changeset - r710:8b0df68370a4

Repositories » LEMON » LEMON-official

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

Changeset - r710:8b0df68370a4

« 709:4d3d1a

711:cc61d0 »

r710:8b0df68370a4

Tue, 12 May 2009 12:06:40

[Not Reviewed]

0 comments (0 inline)

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

Fix the GEQ/LEQ handling in NetworkSimplex + improve doc (#291) - Fix the optimality conditions for the GEQ/LEQ form. - Fix the initialization of the algortihm. It ensures correct solutions and it is much faster for the inequality forms. - Fix the pivot rules to search all the arcs that have to be allowed to get in the basis. - Better block size for the Block Search pivot rule. - Improve documentation of the problem and move it to a separate page.

0 3 1

default

4 files changed with 338 insertions and 188 deletions:

doc/min_cost_flow.dox

153

doc/Makefile.am

doc/groups.dox

lemon/network_simplex.h

180

105

↑ Collapse diff ↑

doc/min_cost_flow.dox

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.
+		 *
 		 */
 		namespace lemon {
 		/**
 		\page min_cost_flow Minimum Cost Flow Problem
 		\section mcf_def Definition (GEQ form)
 		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 (lower and upper bounds)
 		and arc costs.
 		Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$,
 		\f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and
 		upper bounds for the flow values on the arcs, for which
 		\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
 		\f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow
 		on the arcs and \f$sup: V\rightarrow\mathbf{R}\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{R}\f$ solution
 		of the following optimization problem.
 		\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]
 		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.
 		\section mcf_algs Algorithms
 		LEMON contains several algorithms for solving this problem, for more
 		information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms".
 		A feasible solution for this problem can be found using \ref Circulation.
 		\section mcf_dual Dual Solution
 		The dual solution of the minimum cost flow problem is represented by
 		node potentials \f$\pi: V\rightarrow\mathbf{R}\f$.
 		An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal
 		if and only if for some \f$\pi: V\rightarrow\mathbf{R}\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$ nodes:
 		   - \f$\pi(u)<=0\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 potential function \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.
 		\section mcf_eq Equality Form
 		The above \ref mcf_def "definition" is actually more general than the
 		usual formulation of the minimum cost flow problem, in which strict
 		equalities are required in the supply/demand contraints.
 		\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) =
 		    sup(u) \quad \forall u\in V \f]
 		\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
 		However if the sum of the supply values is zero, then these two problems
 		are equivalent.
 		The \ref min_cost_flow_algs "algorithms" in LEMON support the general
 		form, so if you need the equality form, you have to ensure this additional
 		contraint manually.
 		\section mcf_leq Opposite Inequalites (LEQ Form)
 		Another possible definition of the minimum cost flow problem is
 		when there are <em>"less or equal"</em> (LEQ) supply/demand constraints,
 		instead of the <em>"greater or equal"</em> (GEQ) constraints.
 		\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.
 		The equality form is also a special case of this form, of course.
 		You could easily transform this case to the \ref mcf_def "GEQ form"
 		of the problem 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.
 		Note that the optimality conditions for this supply constraint type are
 		slightly differ from the conditions that are discussed for the GEQ form,
 		namely the potentials have to be non-negative instead of non-positive.
 		An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem
 		is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$
 		node potentials the following 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$ nodes:
 		   - \f$\pi(u)>=0\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$.
 		*/
+		}

doc/Makefile.am

Show inline comments

 		EXTRA_DIST += \
 			doc/Doxyfile.in \
 			doc/DoxygenLayout.xml \
 			doc/coding_style.dox \
 			doc/dirs.dox \
 			doc/groups.dox \
 			doc/lgf.dox \
 			doc/license.dox \
 			doc/mainpage.dox \
 			doc/migration.dox \
 			doc/min_cost_flow.dox \
 			doc/named-param.dox \
 			doc/namespaces.dox \
 			doc/html \
 			doc/CMakeLists.txt
 		DOC_EPS_IMAGES18 = \
 			grid_graph.eps \
 			nodeshape_0.eps \
 			nodeshape_1.eps \
 			nodeshape_2.eps \
 			nodeshape_3.eps \
 			nodeshape_4.eps
 		DOC_EPS_IMAGES27 = \
 			bipartite_matching.eps \
 			bipartite_partitions.eps \
 			connected_components.eps \
 			edge_biconnected_components.eps \
 			node_biconnected_components.eps \
 			strongly_connected_components.eps
 		DOC_EPS_IMAGES = \
 			$(DOC_EPS_IMAGES18) \
 			$(DOC_EPS_IMAGES27)
 		DOC_PNG_IMAGES = \
 			$(DOC_EPS_IMAGES:%.eps=doc/gen-images/%.png)
 		EXTRA_DIST += $(DOC_EPS_IMAGES:%=doc/images/%)
 		doc/html:
 			$(MAKE) $(AM_MAKEFLAGS) html
 		GS_COMMAND=gs -dNOPAUSE -dBATCH -q -dEPSCrop -dTextAlphaBits=4 -dGraphicsAlphaBits=4
 		$(DOC_EPS_IMAGES18:%.eps=doc/gen-images/%.png): doc/gen-images/%.png: doc/images/%.eps
 			-mkdir doc/gen-images
 			if test ${gs_found} = yes; then \
 			  $(GS_COMMAND) -sDEVICE=pngalpha -r18 -sOutputFile=$@ $<; \
 			else \
 			  echo; \
 			  echo "Ghostscript not found."; \
 			  echo; \
 			  exit 1; \
 			fi
 		$(DOC_EPS_IMAGES27:%.eps=doc/gen-images/%.png): doc/gen-images/%.png: doc/images/%.eps
 			-mkdir doc/gen-images
 			if test ${gs_found} = yes; then \
 			  $(GS_COMMAND) -sDEVICE=pngalpha -r27 -sOutputFile=$@ $<; \
 			else \
 			  echo; \
 			  echo "Ghostscript not found."; \
 			  echo; \
 			  exit 1; \
 			fi
 		html-local: $(DOC_PNG_IMAGES)
 			if test ${doxygen_found} = yes; then \
 			  cd doc; \
 			  doxygen Doxyfile; \
 			  cd ..; \
 			else \
 			  echo; \
 			  echo "Doxygen not found."; \
 			  echo; \
 			  exit 1; \
 			fi
 		clean-local:
 			-rm -rf doc/html
 			-rm -f doc/doxygen.log
 			-rm -f $(DOC_PNG_IMAGES)
 			-rm -rf doc/gen-images
 		update-external-tags:
 			wget -O doc/libstdc++.tag.tmp http://gcc.gnu.org/onlinedocs/libstdc++/latest-doxygen/libstdc++.tag && \
 			mv doc/libstdc++.tag.tmp doc/libstdc++.tag || \
 			rm doc/libstdc++.tag.tmp
 		install-html-local: doc/html
 			@$(NORMAL_INSTALL)
 			$(mkinstalldirs) $(DESTDIR)$(htmldir)/docs
 			for p in doc/html/*.{html,css,png,map,gif,tag} ; do \
 			  f="`echo $$p | sed -e 's|^.*/||'`"; \
 			  echo " $(INSTALL_DATA) $$p $(DESTDIR)$(htmldir)/docs/$$f"; \
 			  $(INSTALL_DATA) $$p $(DESTDIR)$(htmldir)/docs/$$f; \
 			done
 		uninstall-local:
 			@$(NORMAL_UNINSTALL)
 			for p in doc/html/*.{html,css,png,map,gif,tag} ; do \
 			  f="`echo $$p | sed -e 's|^.*/||'`"; \
 			  echo " rm -f $(DESTDIR)$(htmldir)/docs/$$f"; \
 			  rm -f $(DESTDIR)$(htmldir)/docs/$$f; \
 			done
 		.PHONY: update-external-tags

doc/groups.dox

Show inline comments

@@ -82,610 +82,531 @@
 		\endcode
 		is needed to run on the reverse oriented graph.  It may be expensive
 		(in time or in memory usage) to copy \c g with the reversed
 		arcs.  In this case, an adaptor class is used, which (according
 		to LEMON \ref concepts::Digraph "digraph concepts") works as a digraph.
 		The adaptor uses the original digraph structure and digraph operations when
 		methods of the reversed oriented graph are called.  This means that the adaptor
 		have minor memory usage, and do not perform sophisticated algorithmic
 		actions.  The purpose of it is to give a tool for the cases when a
 		graph have to be used in a specific alteration.  If this alteration is
 		obtained by a usual construction like filtering the node or the arc set or
 		considering a new orientation, then an adaptor is worthwhile to use.
 		To come back to the reverse oriented graph, in this situation
 		\code
 		template<typename Digraph> class ReverseDigraph;
 		\endcode
 		template class can be used. The code looks as follows
 		\code
 		ListDigraph g;
 		ReverseDigraph<ListDigraph> rg(g);
 		int result = algorithm(rg);
 		\endcode
 		During running the algorithm, the original digraph \c g is untouched.
 		This techniques give rise to an elegant code, and based on stable
 		graph adaptors, complex algorithms can be implemented easily.
 		In flow, circulation and matching problems, the residual
 		graph is of particular importance. Combining an adaptor implementing
 		this with shortest path algorithms or minimum mean cycle algorithms,
 		a range of weighted and cardinality optimization algorithms can be
 		obtained. For other examples, the interested user is referred to the
 		detailed documentation of particular adaptors.
 		The behavior of graph adaptors can be very different. Some of them keep
 		capabilities of the original graph while in other cases this would be
 		meaningless. This means that the concepts that they meet depend
 		on the graph adaptor, and the wrapped graph.
 		For example, if an arc of a reversed digraph is deleted, this is carried
 		out by deleting the corresponding arc of the original digraph, thus the
 		adaptor modifies the original digraph.
 		However in case of a residual digraph, this operation has no sense.
 		Let us stand one more example here to simplify your work.
 		ReverseDigraph has constructor
 		\code
 		ReverseDigraph(Digraph& digraph);
 		\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 maps Maps
 		@ingroup datas
 		\brief Map structures implemented in LEMON.
 		This group contains 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 contains 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 contains 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 contains two dimensional data storages implemented in LEMON.
 		*/
 		/**
 		@defgroup paths Path Structures
 		@ingroup datas
 		\brief %Path structures implemented in LEMON.
 		This group contains 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 contains some data structures implemented in LEMON in
 		order to make it easier to implement combinatorial algorithms.
 		*/
 		/**
 		@defgroup algs Algorithms
 		\brief This group contains the several algorithms
 		implemented in LEMON.
 		This group contains the several algorithms
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup search Graph Search
 		@ingroup algs
 		\brief Common graph search algorithms.
 		This group contains 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 contains 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 contains 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
 		\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
 		following optimization problem.
 		\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
 		also provide functions to query the minimum cut, which is the dual
 		problem of maximum flow.
 		\ref Circulation is a preflow push-relabel algorithm implemented directly
 		for finding feasible circulations, which is a somewhat different problem,
 		but it is strongly related to maximum flow.
 		For more information, see \ref Circulation.
 		*/
 		/**
-		@defgroup min_cost_flow Minimum Cost Flow Algorithms
+		@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
 		@ingroup algs
 		\brief Algorithms for finding minimum cost flows and circulations.
 		This group contains the algorithms for finding minimum cost flows and
 		circulations.
+		circulations. For more information about this problem and its dual
 		solution see \ref min_cost_flow "Minimum Cost Flow Problem".
 		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 (lower and upper bounds)
 		and arc costs.
 		Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
 		\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
 		upper bounds for the flow values on the arcs, for which
 		\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
 		\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
 		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}\f$ solution
 		of the following optimization problem.
 		\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]
 		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}\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$ nodes:
 		   - 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 potential function \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.
+		LEMON contains several algorithms for this problem.
 		 - \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 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 contains 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 GomoryHu "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_properties Connectivity and Other Graph Properties
 		@ingroup algs
 		\brief Algorithms for discovering the graph properties
 		This group contains 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 contains 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 the algorithms for calculating
 		matchings in graphs and bipartite graphs. The general matching problem is
 		finding a subset of the edges for which each node has at most one incident
 		edge.
 		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 minimum cost spanning trees and arborescences.
 		This group contains the algorithms for finding minimum cost spanning
 		trees and arborescences.
 		*/
 		/**
 		@defgroup auxalg Auxiliary Algorithms
 		@ingroup algs
 		\brief Auxiliary algorithms implemented in LEMON.
 		This group contains 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 contains the approximation and heuristic algorithms
 		implemented in LEMON.
 		*/
 		/**
 		@defgroup gen_opt_group General Optimization Tools
 		\brief This group contains some general optimization frameworks
 		implemented in LEMON.
 		This group contains 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 contains 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 contains 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 contains some simple basic graph utilities.
 		*/
 		/**
 		@defgroup misc Miscellaneous Tools
 		@ingroup utils
 		\brief Tools for development, debugging and testing.
 		This group contains 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 contains simple tools for measuring the performance
 		of algorithms.
 		*/
 		/**
 		@defgroup exceptions Exceptions
 		@ingroup utils
 		\brief Exceptions defined in LEMON.
 		This group contains the exceptions defined in LEMON.
 		*/
 		/**
 		@defgroup io_group Input-Output
 		\brief Graph Input-Output methods
 		This group contains the tools for importing and exporting graphs
 		and graph related data. Now it supports the \ref lgf-format
 		"LEMON Graph Format", the \c DIMACS format and the encapsulated
 		postscript (EPS) format.
 		*/
 		/**
 		@defgroup lemon_io LEMON Graph Format
 		@ingroup io_group
 		\brief Reading and writing LEMON Graph Format.
 		This group contains methods for reading and writing
 		\ref lgf-format "LEMON Graph Format".
 		*/
 		/**
 		@defgroup eps_io Postscript Exporting
 		@ingroup io_group
 		\brief General \c EPS drawer and graph exporter
 		This group contains general \c EPS drawing methods and special
 		graph exporting tools.
 		*/
 		/**
 		@defgroup dimacs_group DIMACS format
 		@ingroup io_group
 		\brief Read and write files in DIMACS format
 		Tools to read a digraph from or write it to a file in DIMACS format data.
 		*/
 		/**
 		@defgroup nauty_group NAUTY Format
 		@ingroup io_group
 		\brief Read \e Nauty format
 		Tool to read graphs from \e Nauty format data.
 		*/
 		/**
 		@defgroup concept Concepts
 		\brief Skeleton classes and concept checking classes
 		This group contains the data/algorithm skeletons and concept checking
 		classes implemented in LEMON.
 		The purpose of the classes in this group is fourfold.

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
+		/// \ingroup min_cost_flow_algs
 		///
 		/// \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>
 		namespace lemon {
-		  /// \addtogroup min_cost_flow
+		  /// \addtogroup min_cost_flow_algs
 		  /// @{
 		  /// \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 directions of the supply/demand inequality
 		  /// constraints. For more information see \ref SupplyType.
 		  ///
 		  /// Most of the parameters of the problem (except for the digraph)
 		  /// can be given using separate functions, and the algorithm can be
 		  /// executed using the \ref run() function. If some parameters are not
 		  /// specified, then default values will be used.
 		  ///
 		  /// \tparam GR The digraph type the algorithm runs on.
 		  /// \tparam V The 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 V.
 		  ///
 		  /// \warning Both value types must be signed and all input data must
 		  /// be integer.
 		  ///
 		  /// \note %NetworkSimplex provides five different pivot rule
 		  /// implementations, from which the most efficient one is used
 		  /// by default. For more information see \ref PivotRule.
 		  template <typename GR, typename V = int, typename C = V>
 		  class NetworkSimplex
+		  {
 		  public:
 		    /// The type of the flow amounts, capacity bounds and supply values
 		    typedef V Value;
 		    /// The type of the arc costs
 		    typedef C Cost;
 		  public:
 		    /// \brief Problem type constants for the \c run() function.
 		    ///
 		    /// Enum type containing the problem type constants that can be
 		    /// returned by the \ref run() function of the algorithm.
 		    enum ProblemType {
 		      /// The problem has no feasible solution (flow).
 		      INFEASIBLE,
 		      /// The problem has optimal solution (i.e. it is feasible and
 		      /// bounded), and the algorithm has found optimal flow and node
 		      /// potentials (primal and dual solutions).
 		      OPTIMAL,
 		      /// The objective function of the problem is unbounded, i.e.
 		      /// there is a directed cycle having negative total cost and
 		      /// infinite upper bound.
 		      UNBOUNDED
 		    };
 		    /// \brief Constants for selecting the type of the supply constraints.
 		    ///
 		    /// Enum type containing constants for selecting the supply 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 supply 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 supplyType()
 		    /// function.
 		    ///
-		    /// Note that the equality form is a special case of both supply types.
+		    /// The default supply type is \c GEQ, the \c LEQ type can be
 		    /// selected using \ref supplyType().
 		    /// The equality form is a special case of both supply types.
 		    enum SupplyType {
 		      /// This option means that there are <em>"greater or equal"</em>
 		      /// supply/demand constraints in the definition, 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.
+		      /// supply/demand constraints in the definition of the problem.
 		      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>
 		      /// supply/demand constraints in the definition, 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
 		      /// supply/demand constraints in the definition of the problem.
 		      LEQ
 		    };
 		    /// \brief Constants for selecting the pivot rule.
 		    ///
 		    /// Enum type containing constants 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
 		    };
 		  private:
 		    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
 		    typedef std::vector<Arc> ArcVector;
 		    typedef std::vector<Node> NodeVector;
 		    typedef std::vector<int> IntVector;
 		    typedef std::vector<bool> BoolVector;
 		    typedef std::vector<Value> ValueVector;
 		    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;
 		    int _all_arc_num;
 		    int _search_arc_num;
 		    // Parameters of the problem
 		    bool _have_lower;
 		    SupplyType _stype;
 		    Value _sum_supply;
 		    // Data structures for storing the digraph
 		    IntNodeMap _node_id;
 		    IntArcMap _arc_id;
 		    IntVector _source;
 		    IntVector _target;
 		    // Node and arc data
 		    ValueVector _lower;
 		    ValueVector _upper;
 		    ValueVector _cap;
 		    CostVector _cost;
 		    ValueVector _supply;
 		    ValueVector _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;
 		    Value delta;
 		  public:
 		    /// \brief Constant for infinite upper bounds (capacities).
 		    ///
 		    /// Constant for infinite upper bounds (capacities).
 		    /// It is \c std::numeric_limits<Value>::infinity() if available,
 		    /// \c std::numeric_limits<Value>::max() otherwise.
 		    const Value INF;
 		  private:
 		    // 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;
+		      int _search_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)
+		        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
 		        _next_arc(0)
 		      {}
 		      // Find next entering arc
 		      bool findEnteringArc() {
 		        Cost c;
-		        for (int e = _next_arc; e < _arc_num; ++e) {
+		        for (int e = _next_arc; e < _search_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;
+		      int _search_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)
+		        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
 		      {}
 		      // Find next entering arc
 		      bool findEnteringArc() {
 		        Cost c, min = 0;
-		        for (int e = 0; e < _arc_num; ++e) {
+		        for (int e = 0; e < _search_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;
+		      int _search_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)
+		        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
 		        _next_arc(0)
+		      {
 		        // The main parameters of the pivot rule
-		        const double BLOCK_SIZE_FACTOR = 2.0;
+		        const double BLOCK_SIZE_FACTOR = 0.5;
 		        const int MIN_BLOCK_SIZE = 10;
 		        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
-		                                    std::sqrt(double(_arc_num))),
+		                                    std::sqrt(double(_search_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) {
+		        for (e = _next_arc; e < _search_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;
+		            }
+		          }
+		        }
 		        if (min >= 0) return false;
 		        _in_arc = min_arc;
 		        _next_arc = e;
 		        return true;
+		      }
 		    }; //class BlockSearchPivotRule
 		    // Implementation of the Candidate List pivot rule
 		    class CandidateListPivotRule
+		    {
 		    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;
+		      int _search_arc_num;
 		      // Pivot rule data
 		      IntVector _candidates;
 		      int _list_length, _minor_limit;
 		      int _curr_length, _minor_count;
 		      int _next_arc;
 		    public:
 		      /// Constructor
 		      CandidateListPivotRule(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)
+		        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
 		        _next_arc(0)
+		      {
 		        // The main parameters of the pivot rule
 		        const double LIST_LENGTH_FACTOR = 1.0;
 		        const int MIN_LIST_LENGTH = 10;
 		        const double MINOR_LIMIT_FACTOR = 0.1;
 		        const int MIN_MINOR_LIMIT = 3;
 		        _list_length = std::max( int(LIST_LENGTH_FACTOR *
-		                                     std::sqrt(double(_arc_num))),
+		                                     std::sqrt(double(_search_arc_num))),
 		                                 MIN_LIST_LENGTH );
 		        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
 		                                 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) {
+		        for (e = _next_arc; e < _search_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;
+		      int _search_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)
 		        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
 		        _next_arc(0), _cand_cost(ns._search_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 *
-		                                    std::sqrt(double(_arc_num))),
+		                                    std::sqrt(double(_search_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) {
+		        for (int e = _next_arc; e < _search_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.
 		    ///
 		    /// The constructor of the class.
 		    ///
 		    /// \param graph The digraph the algorithm runs on.
 		    NetworkSimplex(const GR& graph) :
 		      _graph(graph), _node_id(graph), _arc_id(graph),
 		      INF(std::numeric_limits<Value>::has_infinity ?
 		          std::numeric_limits<Value>::infinity() :
 		          std::numeric_limits<Value>::max())
+		    {
 		      // Check the value types
 		      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
 		        "The flow type of NetworkSimplex must be signed");
 		      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
 		        "The cost type of NetworkSimplex must be signed");
 		      // Resize vectors
 		      _node_num = countNodes(_graph);
 		      _arc_num = countArcs(_graph);
 		      int all_node_num = _node_num + 1;
-		      int all_arc_num = _arc_num + _node_num;
+		      int max_arc_num = _arc_num + 2 * _node_num;
 		      _source.resize(all_arc_num);
 		      _target.resize(all_arc_num);
 		      _source.resize(max_arc_num);
 		      _target.resize(max_arc_num);
 		      _lower.resize(all_arc_num);
 		      _upper.resize(all_arc_num);
 		      _cap.resize(all_arc_num);
 		      _cost.resize(all_arc_num);
 		      _lower.resize(_arc_num);
 		      _upper.resize(_arc_num);
 		      _cap.resize(max_arc_num);
 		      _cost.resize(max_arc_num);
 		      _supply.resize(all_node_num);
-		      _flow.resize(all_arc_num);
+		      _flow.resize(max_arc_num);
 		      _pi.resize(all_node_num);
 		      _parent.resize(all_node_num);
 		      _pred.resize(all_node_num);
 		      _forward.resize(all_node_num);
 		      _thread.resize(all_node_num);
 		      _rev_thread.resize(all_node_num);
 		      _succ_num.resize(all_node_num);
 		      _last_succ.resize(all_node_num);
-		      _state.resize(all_arc_num);
+		      _state.resize(max_arc_num);
 		      // Copy the graph (store the arcs in a mixed order)
 		      int i = 0;
 		      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
 		        _node_id[n] = i;
+		      }
 		      int k = std::max(int(std::sqrt(double(_arc_num))), 10);
 		      i = 0;
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        _arc_id[a] = i;
 		        _source[i] = _node_id[_graph.source(a)];
 		        _target[i] = _node_id[_graph.target(a)];
 		        if ((i += k) >= _arc_num) i = (i % k) + 1;
+		      }
 		      // Initialize maps
 		      for (int i = 0; i != _node_num; ++i) {
 		        _supply[i] = 0;
+		      }
 		      for (int i = 0; i != _arc_num; ++i) {
 		        _lower[i] = 0;
 		        _upper[i] = INF;
 		        _cost[i] = 1;
+		      }
 		      _have_lower = false;
 		      _stype = GEQ;
+		    }
 		    /// \name Parameters
 		    /// The parameters of the algorithm can be specified using these
 		    /// functions.
 		    /// @{
 		    /// \brief Set the lower bounds on the arcs.
 		    ///
 		    /// This function sets the lower bounds on the arcs.
 		    /// If it is not used before calling \ref run(), the lower bounds
 		    /// will be set to zero on all arcs.
 		    ///
 		    /// \param map An arc map storing the lower bounds.
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template <typename LowerMap>
 		    NetworkSimplex& lowerMap(const LowerMap& map) {
 		      _have_lower = true;
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        _lower[_arc_id[a]] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the upper bounds (capacities) on the arcs.
 		    ///
 		    /// This function sets the upper bounds (capacities) on the arcs.
 		    /// If it is not used before calling \ref run(), the upper bounds
 		    /// will be set to \ref INF on all arcs (i.e. the flow value will be
 		    /// unbounded from above on each arc).
 		    ///
 		    /// \param map An arc map storing the upper bounds.
 		    /// Its \c Value type must be convertible to the \c Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename UpperMap>
 		    NetworkSimplex& upperMap(const UpperMap& map) {
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        _upper[_arc_id[a]] = map[a];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set the costs of the arcs.
 		    ///
 		    /// This function sets the costs of the arcs.
 		    /// If it is not used before calling \ref run(), the costs
 		    /// will be set to \c 1 on all arcs.
 		    ///
 		    /// \param map An arc map storing the costs.
 		    /// Its \c Value type must be convertible to the \c Cost type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename CostMap>
 		    NetworkSimplex& costMap(const CostMap& map) {
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        _cost[_arc_id[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 Value type
 		    /// of the algorithm.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    template<typename SupplyMap>
 		    NetworkSimplex& supplyMap(const SupplyMap& map) {
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        _supply[_node_id[n]] = map[n];
+		      }
 		      return *this;
+		    }
 		    /// \brief Set single source and target nodes and a supply value.
 		    ///
 		    /// This function sets a single source node and a single target node
 		    /// and the required flow value.
 		    /// If neither this function nor \ref supplyMap() is used before
 		    /// calling \ref run(), the supply of each node will be set to zero.
 		    /// (It makes sense only if non-zero lower bounds are given.)
 		    ///
 		    /// Using this function has the same effect as using \ref supplyMap()
 		    /// with such a map in which \c k is assigned to \c s, \c -k is
 		    /// assigned to \c t and all other nodes have zero supply value.
 		    ///
 		    /// \param s The source node.
 		    /// \param t The target node.
 		    /// \param k The required amount of flow from node \c s to node \c t
 		    /// (i.e. the supply of \c s and the demand of \c t).
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
 		      for (int i = 0; i != _node_num; ++i) {
 		        _supply[i] = 0;
+		      }
 		      _supply[_node_id[s]] =  k;
 		      _supply[_node_id[t]] = -k;
 		      return *this;
+		    }
 		    /// \brief Set the type of the supply constraints.
 		    ///
 		    /// This function sets the type of the supply/demand constraints.
 		    /// If it is not used before calling \ref run(), the \ref GEQ supply
 		    /// type will be used.
 		    ///
 		    /// For more information see \ref SupplyType.
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& supplyType(SupplyType supply_type) {
 		      _stype = supply_type;
 		      return *this;
+		    }
 		    /// @}
 		    /// \name Execution Control
 		    /// The algorithm can be executed using \ref run().
 		    /// @{
 		    /// \brief Run the algorithm.
 		    ///
 		    /// This function runs the algorithm.
 		    /// The paramters can be specified using functions \ref lowerMap(),
 		    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
 		    /// \ref supplyType().
 		    /// For example,
 		    /// \code
 		    ///   NetworkSimplex<ListDigraph> ns(graph);
 		    ///   ns.lowerMap(lower).upperMap(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.
 		    /// However the underlying digraph must not be modified after this
 		    /// class have been constructed, since it copies and extends the graph.
 		    ///
 		    /// \param pivot_rule The pivot rule that will be used during the
 		    /// algorithm. For more information see \ref PivotRule.
 		    ///
 		    /// \return \c INFEASIBLE if no feasible flow exists,
 		    /// \n \c OPTIMAL if the problem has optimal solution
 		    /// (i.e. it is feasible and bounded), and the algorithm has found
 		    /// optimal flow and node potentials (primal and dual solutions),
 		    /// \n \c UNBOUNDED if the objective function of the problem is
 		    /// unbounded, i.e. there is a directed cycle having negative total
 		    /// cost and infinite upper bound.
 		    ///
 		    /// \see ProblemType, PivotRule
 		    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
 		      if (!init()) return INFEASIBLE;
 		      return start(pivot_rule);
+		    }
 		    /// \brief Reset all the parameters that have been given before.
 		    ///
 		    /// This function resets all the paramaters that have been given
 		    /// before using functions \ref lowerMap(), \ref upperMap(),
 		    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
 		    ///
 		    /// It is useful for multiple run() calls. If this function is not
 		    /// used, all the parameters given before are kept for the next
 		    /// \ref run() call.
 		    /// However the underlying digraph must not be modified after this
 		    /// class have been constructed, since it copies and extends the graph.
 		    ///
 		    /// For example,
 		    /// \code
 		    ///   NetworkSimplex<ListDigraph> ns(graph);
 		    ///
 		    ///   // First run
 		    ///   ns.lowerMap(lower).upperMap(upper).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.upperMap(capacity).costMap(cost)
 		    ///     .supplyMap(sup).run();
 		    /// \endcode
 		    ///
 		    /// \return <tt>(*this)</tt>
 		    NetworkSimplex& reset() {
 		      for (int i = 0; i != _node_num; ++i) {
 		        _supply[i] = 0;
+		      }
 		      for (int i = 0; i != _arc_num; ++i) {
 		        _lower[i] = 0;
 		        _upper[i] = INF;
 		        _cost[i] = 1;
+		      }
 		      _have_lower = false;
 		      _stype = GEQ;
 		      return *this;
+		    }
 		    /// @}
 		    /// \name Query Functions
 		    /// The results of the algorithm can be obtained using these
 		    /// functions.\n
 		    /// The \ref run() function must be called before using them.
 		    /// @{
 		    /// \brief Return the total cost of the found flow.
 		    ///
 		    /// This function returns the total cost of the found flow.
 		    /// Its complexity is O(e).
 		    ///
 		    /// \note The return type of the function can be specified as a
 		    /// template parameter. For example,
 		    /// \code
 		    ///   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 Number>
 		    Number totalCost() const {
 		      Number c = 0;
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        int i = _arc_id[a];
 		        c += Number(_flow[i]) * Number(_cost[i]);
+		      }
 		      return c;
+		    }
 		#ifndef DOXYGEN
 		    Cost totalCost() const {
 		      return totalCost<Cost>();
+		    }
 		#endif
 		    /// \brief Return the flow on the given arc.
 		    ///
 		    /// This function returns the flow on the given arc.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Value flow(const Arc& a) const {
 		      return _flow[_arc_id[a]];
+		    }
 		    /// \brief Return the flow map (the primal solution).
 		    ///
 		    /// This function copies the flow value on each arc into the given
 		    /// map. The \c Value type of the algorithm must be convertible to
 		    /// the \c Value type of the map.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    template <typename FlowMap>
 		    void flowMap(FlowMap &map) const {
 		      for (ArcIt a(_graph); a != INVALID; ++a) {
 		        map.set(a, _flow[_arc_id[a]]);
+		      }
+		    }
 		    /// \brief Return the potential (dual value) of the given node.
 		    ///
 		    /// This function returns the potential (dual value) of the
 		    /// given node.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    Cost potential(const Node& n) const {
 		      return _pi[_node_id[n]];
+		    }
 		    /// \brief Return the potential map (the dual solution).
 		    ///
 		    /// This function copies the potential (dual value) of each node
 		    /// into the given map.
 		    /// The \c Cost type of the algorithm must be convertible to the
 		    /// \c Value type of the map.
 		    ///
 		    /// \pre \ref run() must be called before using this function.
 		    template <typename PotentialMap>
 		    void potentialMap(PotentialMap &map) const {
 		      for (NodeIt n(_graph); n != INVALID; ++n) {
 		        map.set(n, _pi[_node_id[n]]);
+		      }
+		    }
 		    /// @}
 		  private:
 		    // Initialize internal data structures
 		    bool init() {
 		      if (_node_num == 0) return false;
 		      // Check the sum of supply values
 		      _sum_supply = 0;
 		      for (int i = 0; i != _node_num; ++i) {
 		        _sum_supply += _supply[i];
+		      }
 		      if ( !((_stype == GEQ && _sum_supply <= 0) ||
 		             (_stype == LEQ && _sum_supply >= 0)) ) return false;
 		      // Remove non-zero lower bounds
 		      if (_have_lower) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Value c = _lower[i];
 		          if (c >= 0) {
 		            _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
 		          } else {
 		            _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
+		          }
 		          _supply[_source[i]] -= c;
 		          _supply[_target[i]] += c;
+		        }
 		      } else {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          _cap[i] = _upper[i];
+		        }
+		      }
 		      // Initialize artifical cost
 		      Cost ART_COST;
 		      if (std::numeric_limits<Cost>::is_exact) {
-		        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
+		        ART_COST = std::numeric_limits<Cost>::max() / 2 + 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;
+		      }
 		      // Initialize arc maps
 		      for (int i = 0; i != _arc_num; ++i) {
 		        _flow[i] = 0;
 		        _state[i] = STATE_LOWER;
+		      }
 		      // 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] = _node_num + 1;
 		      _last_succ[_root] = _root - 1;
 		      _supply[_root] = -_sum_supply;
-		      _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
 		      _pi[_root] = 0;
 		      // Add artificial arcs and initialize the spanning tree data structure
 		      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		        _parent[u] = _root;
 		        _pred[u] = e;
 		        _thread[u] = u + 1;
 		        _rev_thread[u + 1] = u;
 		        _succ_num[u] = 1;
 		        _last_succ[u] = u;
 		        _cost[e] = ART_COST;
 		        _cap[e] = INF;
 		        _state[e] = STATE_TREE;
 		        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
 		          _flow[e] = _supply[u];
 		          _forward[u] = true;
 		          _pi[u] = -ART_COST + _pi[_root];
 		        } else {
 		          _flow[e] = -_supply[u];
 		          _forward[u] = false;
 		          _pi[u] = ART_COST + _pi[_root];
 		      if (_sum_supply == 0) {
 		        // EQ supply constraints
 		        _search_arc_num = _arc_num;
 		        _all_arc_num = _arc_num + _node_num;
 		        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		          _parent[u] = _root;
 		          _pred[u] = e;
 		          _thread[u] = u + 1;
 		          _rev_thread[u + 1] = u;
 		          _succ_num[u] = 1;
 		          _last_succ[u] = u;
 		          _cap[e] = INF;
 		          _state[e] = STATE_TREE;
 		          if (_supply[u] >= 0) {
 		            _forward[u] = true;
 		            _pi[u] = 0;
 		            _source[e] = u;
 		            _target[e] = _root;
 		            _flow[e] = _supply[u];
 		            _cost[e] = 0;
 		          } else {
 		            _forward[u] = false;
 		            _pi[u] = ART_COST;
 		            _source[e] = _root;
 		            _target[e] = u;
 		            _flow[e] = -_supply[u];
 		            _cost[e] = ART_COST;
+		          }
+		        }
+		      }
 		      else if (_sum_supply > 0) {
 		        // LEQ supply constraints
 		        _search_arc_num = _arc_num + _node_num;
 		        int f = _arc_num + _node_num;
 		        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		          _parent[u] = _root;
 		          _thread[u] = u + 1;
 		          _rev_thread[u + 1] = u;
 		          _succ_num[u] = 1;
 		          _last_succ[u] = u;
 		          if (_supply[u] >= 0) {
 		            _forward[u] = true;
 		            _pi[u] = 0;
 		            _pred[u] = e;
 		            _source[e] = u;
 		            _target[e] = _root;
 		            _cap[e] = INF;
 		            _flow[e] = _supply[u];
 		            _cost[e] = 0;
 		            _state[e] = STATE_TREE;
 		          } else {
 		            _forward[u] = false;
 		            _pi[u] = ART_COST;
 		            _pred[u] = f;
 		            _source[f] = _root;
 		            _target[f] = u;
 		            _cap[f] = INF;
 		            _flow[f] = -_supply[u];
 		            _cost[f] = ART_COST;
 		            _state[f] = STATE_TREE;
 		            _source[e] = u;
 		            _target[e] = _root;
 		            _cap[e] = INF;
 		            _flow[e] = 0;
 		            _cost[e] = 0;
 		            _state[e] = STATE_LOWER;
 		            ++f;
+		          }
+		        }
 		        _all_arc_num = f;
+		      }
 		      else {
 		        // GEQ supply constraints
 		        _search_arc_num = _arc_num + _node_num;
 		        int f = _arc_num + _node_num;
 		        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		          _parent[u] = _root;
 		          _thread[u] = u + 1;
 		          _rev_thread[u + 1] = u;
 		          _succ_num[u] = 1;
 		          _last_succ[u] = u;
 		          if (_supply[u] <= 0) {
 		            _forward[u] = false;
 		            _pi[u] = 0;
 		            _pred[u] = e;
 		            _source[e] = _root;
 		            _target[e] = u;
 		            _cap[e] = INF;
 		            _flow[e] = -_supply[u];
 		            _cost[e] = 0;
 		            _state[e] = STATE_TREE;
 		          } else {
 		            _forward[u] = true;
 		            _pi[u] = -ART_COST;
 		            _pred[u] = f;
 		            _source[f] = u;
 		            _target[f] = _root;
 		            _cap[f] = INF;
 		            _flow[f] = _supply[u];
 		            _state[f] = STATE_TREE;
 		            _cost[f] = ART_COST;
 		            _source[e] = _root;
 		            _target[e] = u;
 		            _cap[e] = INF;
 		            _flow[e] = 0;
 		            _cost[e] = 0;
 		            _state[e] = STATE_LOWER;
 		            ++f;
+		          }
+		        }
 		        _all_arc_num = f;
+		      }
 		      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;
 		      Value 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] == INF ? INF : _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] == INF ? INF : _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) {
 		        Value 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
 		    ProblemType 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 INFEASIBLE; // avoid warning
+		    }
 		    template <typename PivotRuleImpl>
 		    ProblemType start() {
 		      PivotRuleImpl pivot(*this);
 		      // Execute the Network Simplex algorithm
 		      while (pivot.findEnteringArc()) {
 		        findJoinNode();
 		        bool change = findLeavingArc();
 		        if (delta >= INF) return UNBOUNDED;
 		        changeFlow(change);
 		        if (change) {
 		          updateTreeStructure();
 		          updatePotential();
+		        }
+		      }
 		      // Check feasibility
 		      if (_sum_supply < 0) {
 		        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
+		        }
+		      }
 		      else if (_sum_supply > 0) {
 		        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
+		        }
+		      }
 		      else {
 		        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
 		          if (_flow[e] != 0) return INFEASIBLE;
+		        }
 		      for (int e = _search_arc_num; e != _all_arc_num; ++e) {
 		        if (_flow[e] != 0) return INFEASIBLE;
+		      }
 		      // Transform the solution and the supply map to the original form
 		      if (_have_lower) {
 		        for (int i = 0; i != _arc_num; ++i) {
 		          Value c = _lower[i];
 		          if (c != 0) {
 		            _flow[i] += c;
 		            _supply[_source[i]] += c;
 		            _supply[_target[i]] -= c;
+		          }
+		        }
+		      }
 		      // Shift potentials to meet the requirements of the GEQ/LEQ type
 		      // optimality conditions
 		      if (_sum_supply == 0) {
 		        if (_stype == GEQ) {
 		          Cost max_pot = std::numeric_limits<Cost>::min();
 		          for (int i = 0; i != _node_num; ++i) {
 		            if (_pi[i] > max_pot) max_pot = _pi[i];
+		          }
 		          if (max_pot > 0) {
 		            for (int i = 0; i != _node_num; ++i)
 		              _pi[i] -= max_pot;
+		          }
 		        } else {
 		          Cost min_pot = std::numeric_limits<Cost>::max();
 		          for (int i = 0; i != _node_num; ++i) {
 		            if (_pi[i] < min_pot) min_pot = _pi[i];
+		          }
 		          if (min_pot < 0) {
 		            for (int i = 0; i != _node_num; ++i)
 		              _pi[i] -= min_pot;
+		          }
+		        }
+		      }
 		      return OPTIMAL;
+		    }
 		  }; //class NetworkSimplex
 		  ///@}
 		} //namespace lemon
 		#endif //LEMON_NETWORK_SIMPLEX_H

0 comments (0 inline)

RhodeCode

Login to your account