author Alpar Juttner Tue, 21 Apr 2009 15:18:54 +0100 changeset 658 85cb3aa71cce parent 647 0ba8dfce7259 parent 657 dacc2cee2b4c child 659 0c8e5c688440 child 660 b1811c363299 child 666 ec817dfc2cb7
Merge and fix
 doc/groups.dox file | annotate | diff | comparison | revisions lemon/Makefile.am file | annotate | diff | comparison | revisions lemon/circulation.h file | annotate | diff | comparison | revisions lemon/preflow.h file | annotate | diff | comparison | revisions test/CMakeLists.txt file | annotate | diff | comparison | revisions test/Makefile.am file | annotate | diff | comparison | revisions test/circulation_test.cc file | annotate | diff | comparison | revisions tools/dimacs-solver.cc file | annotate | diff | comparison | revisions
     1.1 --- a/doc/groups.dox	Tue Apr 21 13:08:19 2009 +0100
1.2 +++ b/doc/groups.dox	Tue Apr 21 15:18:54 2009 +0100
1.3 @@ -317,15 +317,15 @@
1.4
1.5  The \e maximum \e flow \e problem is to find a flow of maximum value between
1.6  a single source and a single target. Formally, there is a \f$G=(V,A)\f$
1.7 -digraph, a \f$cap:A\rightarrow\mathbf{R}^+_0\f$ capacity function and
1.8 +digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
1.9  \f$s, t \in V\f$ source and target nodes.
1.10 -A maximum flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of the
1.11 +A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
1.12  following optimization problem.
1.13
1.14 -\f[ \max\sum_{a\in\delta_{out}(s)}f(a) - \sum_{a\in\delta_{in}(s)}f(a) \f]
1.15 -\f[ \sum_{a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)
1.16 -    \qquad \forall v\in V\setminus\{s,t\} \f]
1.17 -\f[ 0 \leq f(a) \leq cap(a) \qquad \forall a\in A \f]
1.18 +\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
1.19 +\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
1.20 +    \quad \forall u\in V\setminus\{s,t\} \f]
1.21 +\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
1.22
1.23  LEMON contains several algorithms for solving maximum flow problems:
1.24  - \ref EdmondsKarp Edmonds-Karp algorithm.
1.25 @@ -350,30 +350,98 @@
1.26
1.27  The \e minimum \e cost \e flow \e problem is to find a feasible flow of
1.28  minimum total cost from a set of supply nodes to a set of demand nodes
1.29 -in a network with capacity constraints and arc costs.
1.30 +in a network with capacity constraints (lower and upper bounds)
1.31 +and arc costs.
1.32  Formally, let \f$G=(V,A)\f$ be a digraph,
1.33  \f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
1.34 -upper bounds for the flow values on the arcs,
1.35 +upper bounds for the flow values on the arcs, for which
1.36 +\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
1.37  \f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
1.38 -on the arcs, and
1.39 -\f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values
1.40 -of the nodes.
1.41 -A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of
1.42 -the following optimization problem.
1.43 +on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
1.44 +signed supply values of the nodes.
1.45 +If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
1.46 +supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
1.47 +\f$-sup(u)\f$ demand.
1.48 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
1.49 +of the following optimization problem.
1.50
1.51 -\f[ \min\sum_{a\in A} f(a) cost(a) \f]
1.52 -\f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a) =
1.53 -    supply(v) \qquad \forall v\in V \f]
1.54 -\f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]
1.55 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.56 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
1.57 +    sup(u) \quad \forall u\in V \f]
1.58 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.59
1.60 -LEMON contains several algorithms for solving minimum cost flow problems:
1.61 - - \ref CycleCanceling Cycle-canceling algorithms.
1.62 - - \ref CapacityScaling Successive shortest path algorithm with optional
1.63 +The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
1.64 +zero or negative in order to have a feasible solution (since the sum
1.65 +of the expressions on the left-hand side of the inequalities is zero).
1.66 +It means that the total demand must be greater or equal to the total
1.67 +supply and all the supplies have to be carried out from the supply nodes,
1.68 +but there could be demands that are not satisfied.
1.69 +If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
1.70 +constraints have to be satisfied with equality, i.e. all demands
1.71 +have to be satisfied and all supplies have to be used.
1.72 +
1.73 +If you need the opposite inequalities in the supply/demand constraints
1.74 +(i.e. the total demand is less than the total supply and all the demands
1.75 +have to be satisfied while there could be supplies that are not used),
1.76 +then you could easily transform the problem to the above form by reversing
1.77 +the direction of the arcs and taking the negative of the supply values
1.78 +(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
1.79 +However \ref NetworkSimplex algorithm also supports this form directly
1.80 +for the sake of convenience.
1.81 +
1.82 +A feasible solution for this problem can be found using \ref Circulation.
1.83 +
1.84 +Note that the above formulation is actually more general than the usual
1.85 +definition of the minimum cost flow problem, in which strict equalities
1.86 +are required in the supply/demand contraints, i.e.
1.87 +
1.88 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
1.89 +    sup(u) \quad \forall u\in V. \f]
1.90 +
1.91 +However if the sum of the supply values is zero, then these two problems
1.92 +are equivalent. So if you need the equality form, you have to ensure this
1.93 +additional contraint for the algorithms.
1.94 +
1.95 +The dual solution of the minimum cost flow problem is represented by node
1.96 +potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
1.97 +An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
1.98 +is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
1.99 +node potentials the following \e complementary \e slackness optimality
1.100 +conditions hold.
1.101 +
1.102 + - For all \f$uv\in A\f$ arcs:
1.103 +   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
1.104 +   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
1.105 +   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
1.106 + - For all \f$u\in V\f$:
1.107 +   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
1.108 +     then \f$\pi(u)=0\f$.
1.109 +
1.110 +Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
1.111 +\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
1.112 +\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
1.113 +
1.114 +All algorithms provide dual solution (node potentials) as well
1.115 +if an optimal flow is found.
1.116 +
1.117 +LEMON contains several algorithms for solving minimum cost flow problems.
1.118 + - \ref NetworkSimplex Primal Network Simplex algorithm with various
1.119 +   pivot strategies.
1.120 + - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
1.121 +   cost scaling.
1.122 + - \ref CapacityScaling Successive Shortest %Path algorithm with optional
1.123     capacity scaling.
1.124 - - \ref CostScaling Push-relabel and augment-relabel algorithms based on
1.125 -   cost scaling.
1.126 - - \ref NetworkSimplex Primal network simplex algorithm with various
1.127 -   pivot strategies.
1.128 + - \ref CancelAndTighten The Cancel and Tighten algorithm.
1.129 + - \ref CycleCanceling Cycle-Canceling algorithms.
1.130 +
1.131 +Most of these implementations support the general inequality form of the
1.132 +minimum cost flow problem, but CancelAndTighten and CycleCanceling
1.133 +only support the equality form due to the primal method they use.
1.134 +
1.135 +In general NetworkSimplex is the most efficient implementation,
1.136 +but in special cases other algorithms could be faster.
1.137 +For example, if the total supply and/or capacities are rather small,
1.138 +CapacityScaling is usually the fastest algorithm (without effective scaling).
1.139  */
1.140
1.141  /**

     2.1 --- a/lemon/Makefile.am	Tue Apr 21 13:08:19 2009 +0100
2.2 +++ b/lemon/Makefile.am	Tue Apr 21 15:18:54 2009 +0100
2.3 @@ -93,6 +93,7 @@
2.4  	lemon/math.h \
2.5  	lemon/min_cost_arborescence.h \
2.7 +	lemon/network_simplex.h \
2.8  	lemon/path.h \
2.9  	lemon/preflow.h \

     3.1 --- a/lemon/circulation.h	Tue Apr 21 13:08:19 2009 +0100
3.2 +++ b/lemon/circulation.h	Tue Apr 21 15:18:54 2009 +0100
3.3 @@ -31,52 +31,52 @@
3.4    /// \brief Default traits class of Circulation class.
3.5    ///
3.6    /// Default traits class of Circulation class.
3.7 -  /// \tparam GR Digraph type.
3.8 -  /// \tparam LM Lower bound capacity map type.
3.9 -  /// \tparam UM Upper bound capacity map type.
3.10 -  /// \tparam DM Delta map type.
3.11 +  ///
3.12 +  /// \tparam GR Type of the digraph the algorithm runs on.
3.13 +  /// \tparam LM The type of the lower bound map.
3.14 +  /// \tparam UM The type of the upper bound (capacity) map.
3.15 +  /// \tparam SM The type of the supply map.
3.16    template <typename GR, typename LM,
3.17 -            typename UM, typename DM>
3.18 +            typename UM, typename SM>
3.19    struct CirculationDefaultTraits {
3.20
3.21      /// \brief The type of the digraph the algorithm runs on.
3.22      typedef GR Digraph;
3.23
3.24 -    /// \brief The type of the map that stores the circulation lower
3.25 -    /// bound.
3.26 +    /// \brief The type of the lower bound map.
3.27      ///
3.28 -    /// The type of the map that stores the circulation lower bound.
3.29 -    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
3.30 -    typedef LM LCapMap;
3.31 +    /// The type of the map that stores the lower bounds on the arcs.
3.32 +    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
3.33 +    typedef LM LowerMap;
3.34
3.35 -    /// \brief The type of the map that stores the circulation upper
3.36 -    /// bound.
3.37 +    /// \brief The type of the upper bound (capacity) map.
3.38      ///
3.39 -    /// The type of the map that stores the circulation upper bound.
3.40 -    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
3.41 -    typedef UM UCapMap;
3.42 +    /// The type of the map that stores the upper bounds (capacities)
3.43 +    /// on the arcs.
3.44 +    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
3.45 +    typedef UM UpperMap;
3.46
3.47 -    /// \brief The type of the map that stores the lower bound for
3.48 -    /// the supply of the nodes.
3.49 +    /// \brief The type of supply map.
3.50      ///
3.51 -    /// The type of the map that stores the lower bound for the supply
3.52 -    /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
3.53 -    /// concept.
3.54 -    typedef DM DeltaMap;
3.55 +    /// The type of the map that stores the signed supply values of the
3.56 +    /// nodes.
3.57 +    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
3.58 +    typedef SM SupplyMap;
3.59
3.60      /// \brief The type of the flow values.
3.61 -    typedef typename DeltaMap::Value Value;
3.62 +    typedef typename SupplyMap::Value Flow;
3.63
3.64      /// \brief The type of the map that stores the flow values.
3.65      ///
3.66      /// The type of the map that stores the flow values.
3.67 -    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
3.68 -    typedef typename Digraph::template ArcMap<Value> FlowMap;
3.69 +    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
3.70 +    /// concept.
3.71 +    typedef typename Digraph::template ArcMap<Flow> FlowMap;
3.72
3.73      /// \brief Instantiates a FlowMap.
3.74      ///
3.75      /// This function instantiates a \ref FlowMap.
3.76 -    /// \param digraph The digraph, to which we would like to define
3.77 +    /// \param digraph The digraph for which we would like to define
3.78      /// the flow map.
3.79      static FlowMap* createFlowMap(const Digraph& digraph) {
3.80        return new FlowMap(digraph);
3.81 @@ -93,7 +93,7 @@
3.82      /// \brief Instantiates an Elevator.
3.83      ///
3.84      /// This function instantiates an \ref Elevator.
3.85 -    /// \param digraph The digraph, to which we would like to define
3.86 +    /// \param digraph The digraph for which we would like to define
3.87      /// the elevator.
3.88      /// \param max_level The maximum level of the elevator.
3.89      static Elevator* createElevator(const Digraph& digraph, int max_level) {
3.90 @@ -103,7 +103,7 @@
3.91      /// \brief The tolerance used by the algorithm
3.92      ///
3.93      /// The tolerance used by the algorithm to handle inexact computation.
3.94 -    typedef lemon::Tolerance<Value> Tolerance;
3.95 +    typedef lemon::Tolerance<Flow> Tolerance;
3.96
3.97    };
3.98
3.99 @@ -111,53 +111,69 @@
3.100       \brief Push-relabel algorithm for the network circulation problem.
3.101
3.102       \ingroup max_flow
3.103 -     This class implements a push-relabel algorithm for the network
3.104 -     circulation problem.
3.105 +     This class implements a push-relabel algorithm for the \e network
3.106 +     \e circulation problem.
3.107       It is to find a feasible circulation when lower and upper bounds
3.108 -     are given for the flow values on the arcs and lower bounds
3.109 -     are given for the supply values of the nodes.
3.110 +     are given for the flow values on the arcs and lower bounds are
3.111 +     given for the difference between the outgoing and incoming flow
3.112 +     at the nodes.
3.113
3.114       The exact formulation of this problem is the following.
3.115       Let \f$G=(V,A)\f$ be a digraph,
3.116 -     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
3.117 -     \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
3.118 -     \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
3.119 -     \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
3.120 -     \geq delta(v) \quad \forall v\in V, \f]
3.121 -     \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
3.122 -     \note \f$delta(v)\f$ specifies a lower bound for the supply of node
3.123 -     \f$v\f$. It can be either positive or negative, however note that
3.124 -     \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
3.125 -     have a feasible solution.
3.126 +     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$ denote the lower and
3.127 +     upper bounds on the arcs, for which \f$0 \leq lower(uv) \leq upper(uv)\f$
3.128 +     holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
3.129 +     denotes the signed supply values of the nodes.
3.130 +     If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
3.131 +     supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
3.132 +     \f$-sup(u)\f$ demand.
3.133 +     A feasible circulation is an \f$f: A\rightarrow\mathbf{R}^+_0\f$
3.134 +     solution of the following problem.
3.135
3.136 -     \note A special case of this problem is when
3.137 -     \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
3.138 -     will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
3.139 -     Thus a feasible solution for the
3.140 -     \ref min_cost_flow "minimum cost flow" problem can be calculated
3.141 -     in this way.
3.142 +     \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
3.143 +     \geq sup(u) \quad \forall u\in V, \f]
3.144 +     \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]
3.145 +
3.146 +     The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
3.147 +     zero or negative in order to have a feasible solution (since the sum
3.148 +     of the expressions on the left-hand side of the inequalities is zero).
3.149 +     It means that the total demand must be greater or equal to the total
3.150 +     supply and all the supplies have to be carried out from the supply nodes,
3.151 +     but there could be demands that are not satisfied.
3.152 +     If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
3.153 +     constraints have to be satisfied with equality, i.e. all demands
3.154 +     have to be satisfied and all supplies have to be used.
3.155 +
3.156 +     If you need the opposite inequalities in the supply/demand constraints
3.157 +     (i.e. the total demand is less than the total supply and all the demands
3.158 +     have to be satisfied while there could be supplies that are not used),
3.159 +     then you could easily transform the problem to the above form by reversing
3.160 +     the direction of the arcs and taking the negative of the supply values
3.161 +     (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
3.162 +
3.163 +     Note that this algorithm also provides a feasible solution for the
3.164 +     \ref min_cost_flow "minimum cost flow problem".
3.165
3.166       \tparam GR The type of the digraph the algorithm runs on.
3.167 -     \tparam LM The type of the lower bound capacity map. The default
3.168 +     \tparam LM The type of the lower bound map. The default
3.169       map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
3.170 -     \tparam UM The type of the upper bound capacity map. The default
3.171 -     map type is \c LM.
3.172 -     \tparam DM The type of the map that stores the lower bound
3.173 -     for the supply of the nodes. The default map type is
3.174 +     \tparam UM The type of the upper bound (capacity) map.
3.175 +     The default map type is \c LM.
3.176 +     \tparam SM The type of the supply map. The default map type is
3.177       \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
3.178    */
3.179  #ifdef DOXYGEN
3.180  template< typename GR,
3.181            typename LM,
3.182            typename UM,
3.183 -          typename DM,
3.184 +          typename SM,
3.185            typename TR >
3.186  #else
3.187  template< typename GR,
3.188            typename LM = typename GR::template ArcMap<int>,
3.189            typename UM = LM,
3.190 -          typename DM = typename GR::template NodeMap<typename UM::Value>,
3.191 -          typename TR = CirculationDefaultTraits<GR, LM, UM, DM> >
3.192 +          typename SM = typename GR::template NodeMap<typename UM::Value>,
3.193 +          typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
3.194  #endif
3.195    class Circulation {
3.196    public:
3.197 @@ -167,15 +183,14 @@
3.198      ///The type of the digraph the algorithm runs on.
3.199      typedef typename Traits::Digraph Digraph;
3.200      ///The type of the flow values.
3.201 -    typedef typename Traits::Value Value;
3.202 +    typedef typename Traits::Flow Flow;
3.203
3.204 -    /// The type of the lower bound capacity map.
3.205 -    typedef typename Traits::LCapMap LCapMap;
3.206 -    /// The type of the upper bound capacity map.
3.207 -    typedef typename Traits::UCapMap UCapMap;
3.208 -    /// \brief The type of the map that stores the lower bound for
3.209 -    /// the supply of the nodes.
3.210 -    typedef typename Traits::DeltaMap DeltaMap;
3.211 +    ///The type of the lower bound map.
3.212 +    typedef typename Traits::LowerMap LowerMap;
3.213 +    ///The type of the upper bound (capacity) map.
3.214 +    typedef typename Traits::UpperMap UpperMap;
3.215 +    ///The type of the supply map.
3.216 +    typedef typename Traits::SupplyMap SupplyMap;
3.217      ///The type of the flow map.
3.218      typedef typename Traits::FlowMap FlowMap;
3.219
3.220 @@ -191,9 +206,9 @@
3.221      const Digraph &_g;
3.222      int _node_num;
3.223
3.224 -    const LCapMap *_lo;
3.225 -    const UCapMap *_up;
3.226 -    const DeltaMap *_delta;
3.227 +    const LowerMap *_lo;
3.228 +    const UpperMap *_up;
3.229 +    const SupplyMap *_supply;
3.230
3.231      FlowMap *_flow;
3.232      bool _local_flow;
3.233 @@ -201,7 +216,7 @@
3.234      Elevator* _level;
3.235      bool _local_level;
3.236
3.237 -    typedef typename Digraph::template NodeMap<Value> ExcessMap;
3.238 +    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
3.239      ExcessMap* _excess;
3.240
3.241      Tolerance _tol;
3.242 @@ -231,9 +246,9 @@
3.243      /// type.
3.244      template <typename T>
3.245      struct SetFlowMap
3.246 -      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
3.247 +      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
3.248                             SetFlowMapTraits<T> > {
3.249 -      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
3.250 +      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
3.251                            SetFlowMapTraits<T> > Create;
3.252      };
3.253
3.254 @@ -257,9 +272,9 @@
3.255      /// \sa SetStandardElevator
3.256      template <typename T>
3.257      struct SetElevator
3.258 -      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
3.259 +      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
3.260                             SetElevatorTraits<T> > {
3.261 -      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
3.262 +      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
3.263                            SetElevatorTraits<T> > Create;
3.264      };
3.265
3.266 @@ -285,9 +300,9 @@
3.267      /// \sa SetElevator
3.268      template <typename T>
3.269      struct SetStandardElevator
3.270 -      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
3.271 +      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
3.272                         SetStandardElevatorTraits<T> > {
3.273 -      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
3.274 +      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
3.275                        SetStandardElevatorTraits<T> > Create;
3.276      };
3.277
3.278 @@ -299,18 +314,20 @@
3.279
3.280    public:
3.281
3.282 -    /// The constructor of the class.
3.283 +    /// Constructor.
3.284
3.285      /// The constructor of the class.
3.286 -    /// \param g The digraph the algorithm runs on.
3.287 -    /// \param lo The lower bound capacity of the arcs.
3.288 -    /// \param up The upper bound capacity of the arcs.
3.289 -    /// \param delta The lower bound for the supply of the nodes.
3.290 -    Circulation(const Digraph &g,const LCapMap &lo,
3.291 -                const UCapMap &up,const DeltaMap &delta)
3.292 -      : _g(g), _node_num(),
3.293 -        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
3.294 -        _level(0), _local_level(false), _excess(0), _el() {}
3.295 +    ///
3.296 +    /// \param graph The digraph the algorithm runs on.
3.297 +    /// \param lower The lower bounds for the flow values on the arcs.
3.298 +    /// \param upper The upper bounds (capacities) for the flow values
3.299 +    /// on the arcs.
3.300 +    /// \param supply The signed supply values of the nodes.
3.301 +    Circulation(const Digraph &graph, const LowerMap &lower,
3.302 +                const UpperMap &upper, const SupplyMap &supply)
3.303 +      : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
3.304 +        _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
3.305 +        _excess(NULL) {}
3.306
3.307      /// Destructor.
3.308      ~Circulation() {
3.309 @@ -350,30 +367,30 @@
3.310
3.311    public:
3.312
3.313 -    /// Sets the lower bound capacity map.
3.314 +    /// Sets the lower bound map.
3.315
3.316 -    /// Sets the lower bound capacity map.
3.317 +    /// Sets the lower bound map.
3.318      /// \return <tt>(*this)</tt>
3.319 -    Circulation& lowerCapMap(const LCapMap& map) {
3.320 +    Circulation& lowerMap(const LowerMap& map) {
3.321        _lo = &map;
3.322        return *this;
3.323      }
3.324
3.325 -    /// Sets the upper bound capacity map.
3.326 +    /// Sets the upper bound (capacity) map.
3.327
3.328 -    /// Sets the upper bound capacity map.
3.329 +    /// Sets the upper bound (capacity) map.
3.330      /// \return <tt>(*this)</tt>
3.331 -    Circulation& upperCapMap(const LCapMap& map) {
3.332 +    Circulation& upperMap(const LowerMap& map) {
3.333        _up = &map;
3.334        return *this;
3.335      }
3.336
3.337 -    /// Sets the lower bound map for the supply of the nodes.
3.338 +    /// Sets the supply map.
3.339
3.340 -    /// Sets the lower bound map for the supply of the nodes.
3.341 +    /// Sets the supply map.
3.342      /// \return <tt>(*this)</tt>
3.343 -    Circulation& deltaMap(const DeltaMap& map) {
3.344 -      _delta = &map;
3.345 +    Circulation& supplyMap(const SupplyMap& map) {
3.346 +      _supply = &map;
3.347        return *this;
3.348      }
3.349
3.350 @@ -453,7 +470,7 @@
3.351        createStructures();
3.352
3.353        for(NodeIt n(_g);n!=INVALID;++n) {
3.354 -        (*_excess)[n] = (*_delta)[n];
3.355 +        (*_excess)[n] = (*_supply)[n];
3.356        }
3.357
3.358        for (ArcIt e(_g);e!=INVALID;++e) {
3.359 @@ -482,7 +499,7 @@
3.360        createStructures();
3.361
3.362        for(NodeIt n(_g);n!=INVALID;++n) {
3.363 -        (*_excess)[n] = (*_delta)[n];
3.364 +        (*_excess)[n] = (*_supply)[n];
3.365        }
3.366
3.367        for (ArcIt e(_g);e!=INVALID;++e) {
3.368 @@ -495,7 +512,7 @@
3.369            (*_excess)[_g.target(e)] += (*_lo)[e];
3.370            (*_excess)[_g.source(e)] -= (*_lo)[e];
3.371          } else {
3.372 -          Value fc = -(*_excess)[_g.target(e)];
3.373 +          Flow fc = -(*_excess)[_g.target(e)];
3.374            _flow->set(e, fc);
3.375            (*_excess)[_g.target(e)] = 0;
3.376            (*_excess)[_g.source(e)] -= fc;
3.377 @@ -528,11 +545,11 @@
3.378        while((act=_level->highestActive())!=INVALID) {
3.379          int actlevel=(*_level)[act];
3.380          int mlevel=_node_num;
3.381 -        Value exc=(*_excess)[act];
3.382 +        Flow exc=(*_excess)[act];
3.383
3.384          for(OutArcIt e(_g,act);e!=INVALID; ++e) {
3.385            Node v = _g.target(e);
3.386 -          Value fc=(*_up)[e]-(*_flow)[e];
3.387 +          Flow fc=(*_up)[e]-(*_flow)[e];
3.388            if(!_tol.positive(fc)) continue;
3.389            if((*_level)[v]<actlevel) {
3.390              if(!_tol.less(fc, exc)) {
3.391 @@ -556,7 +573,7 @@
3.392          }
3.393          for(InArcIt e(_g,act);e!=INVALID; ++e) {
3.394            Node v = _g.source(e);
3.395 -          Value fc=(*_flow)[e]-(*_lo)[e];
3.396 +          Flow fc=(*_flow)[e]-(*_lo)[e];
3.397            if(!_tol.positive(fc)) continue;
3.398            if((*_level)[v]<actlevel) {
3.399              if(!_tol.less(fc, exc)) {
3.400 @@ -632,7 +649,7 @@
3.401      ///
3.402      /// \pre Either \ref run() or \ref init() must be called before
3.403      /// using this function.
3.404 -    Value flow(const Arc& arc) const {
3.405 +    Flow flow(const Arc& arc) const {
3.406        return (*_flow)[arc];
3.407      }
3.408
3.409 @@ -651,8 +668,8 @@
3.410
3.411         Barrier is a set \e B of nodes for which
3.412
3.413 -       \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
3.414 -           \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
3.415 +       \f[ \sum_{uv\in A: u\in B} upper(uv) -
3.416 +           \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
3.417
3.418         holds. The existence of a set with this property prooves that a
3.419         feasible circualtion cannot exist.
3.420 @@ -715,7 +732,7 @@
3.421          if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
3.422        for(NodeIt n(_g);n!=INVALID;++n)
3.423          {
3.424 -          Value dif=-(*_delta)[n];
3.425 +          Flow dif=-(*_supply)[n];
3.426            for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
3.427            for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
3.428            if(_tol.negative(dif)) return false;
3.429 @@ -730,10 +747,10 @@
3.430      ///\sa barrierMap()
3.431      bool checkBarrier() const
3.432      {
3.433 -      Value delta=0;
3.434 +      Flow delta=0;
3.435        for(NodeIt n(_g);n!=INVALID;++n)
3.436          if(barrier(n))
3.437 -          delta-=(*_delta)[n];
3.438 +          delta-=(*_supply)[n];
3.439        for(ArcIt e(_g);e!=INVALID;++e)
3.440          {
3.441            Node s=_g.source(e);

     4.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/lemon/network_simplex.h	Tue Apr 21 15:18:54 2009 +0100
4.3 @@ -0,0 +1,1582 @@
4.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
4.5 + *
4.6 + * This file is a part of LEMON, a generic C++ optimization library.
4.7 + *
4.8 + * Copyright (C) 2003-2009
4.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
4.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
4.11 + *
4.12 + * Permission to use, modify and distribute this software is granted
4.13 + * provided that this copyright notice appears in all copies. For
4.14 + * precise terms see the accompanying LICENSE file.
4.15 + *
4.16 + * This software is provided "AS IS" with no warranty of any kind,
4.17 + * express or implied, and with no claim as to its suitability for any
4.18 + * purpose.
4.19 + *
4.20 + */
4.21 +
4.22 +#ifndef LEMON_NETWORK_SIMPLEX_H
4.23 +#define LEMON_NETWORK_SIMPLEX_H
4.24 +
4.25 +/// \ingroup min_cost_flow
4.26 +///
4.27 +/// \file
4.28 +/// \brief Network Simplex algorithm for finding a minimum cost flow.
4.29 +
4.30 +#include <vector>
4.31 +#include <limits>
4.32 +#include <algorithm>
4.33 +
4.34 +#include <lemon/core.h>
4.35 +#include <lemon/math.h>
4.36 +#include <lemon/maps.h>
4.37 +#include <lemon/circulation.h>
4.39 +
4.40 +namespace lemon {
4.41 +
4.42 +  /// \addtogroup min_cost_flow
4.43 +  /// @{
4.44 +
4.45 +  /// \brief Implementation of the primal Network Simplex algorithm
4.46 +  /// for finding a \ref min_cost_flow "minimum cost flow".
4.47 +  ///
4.48 +  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
4.49 +  /// for finding a \ref min_cost_flow "minimum cost flow".
4.50 +  /// This algorithm is a specialized version of the linear programming
4.51 +  /// simplex method directly for the minimum cost flow problem.
4.52 +  /// It is one of the most efficient solution methods.
4.53 +  ///
4.54 +  /// In general this class is the fastest implementation available
4.55 +  /// in LEMON for the minimum cost flow problem.
4.56 +  /// Moreover it supports both direction of the supply/demand inequality
4.57 +  /// constraints. For more information see \ref ProblemType.
4.58 +  ///
4.59 +  /// \tparam GR The digraph type the algorithm runs on.
4.60 +  /// \tparam F The value type used for flow amounts, capacity bounds
4.61 +  /// and supply values in the algorithm. By default it is \c int.
4.62 +  /// \tparam C The value type used for costs and potentials in the
4.63 +  /// algorithm. By default it is the same as \c F.
4.64 +  ///
4.65 +  /// \warning Both value types must be signed and all input data must
4.66 +  /// be integer.
4.67 +  ///
4.68 +  /// \note %NetworkSimplex provides five different pivot rule
4.69 +  /// implementations, from which the most efficient one is used
4.70 +  /// by default. For more information see \ref PivotRule.
4.71 +  template <typename GR, typename F = int, typename C = F>
4.72 +  class NetworkSimplex
4.73 +  {
4.74 +  public:
4.75 +
4.76 +    /// The flow type of the algorithm
4.77 +    typedef F Flow;
4.78 +    /// The cost type of the algorithm
4.79 +    typedef C Cost;
4.80 +#ifdef DOXYGEN
4.81 +    /// The type of the flow map
4.82 +    typedef GR::ArcMap<Flow> FlowMap;
4.83 +    /// The type of the potential map
4.84 +    typedef GR::NodeMap<Cost> PotentialMap;
4.85 +#else
4.86 +    /// The type of the flow map
4.87 +    typedef typename GR::template ArcMap<Flow> FlowMap;
4.88 +    /// The type of the potential map
4.89 +    typedef typename GR::template NodeMap<Cost> PotentialMap;
4.90 +#endif
4.91 +
4.92 +  public:
4.93 +
4.94 +    /// \brief Enum type for selecting the pivot rule.
4.95 +    ///
4.96 +    /// Enum type for selecting the pivot rule for the \ref run()
4.97 +    /// function.
4.98 +    ///
4.99 +    /// \ref NetworkSimplex provides five different pivot rule
4.100 +    /// implementations that significantly affect the running time
4.101 +    /// of the algorithm.
4.102 +    /// By default \ref BLOCK_SEARCH "Block Search" is used, which
4.103 +    /// proved to be the most efficient and the most robust on various
4.104 +    /// test inputs according to our benchmark tests.
4.105 +    /// However another pivot rule can be selected using the \ref run()
4.106 +    /// function with the proper parameter.
4.107 +    enum PivotRule {
4.108 +
4.109 +      /// The First Eligible pivot rule.
4.110 +      /// The next eligible arc is selected in a wraparound fashion
4.111 +      /// in every iteration.
4.112 +      FIRST_ELIGIBLE,
4.113 +
4.114 +      /// The Best Eligible pivot rule.
4.115 +      /// The best eligible arc is selected in every iteration.
4.116 +      BEST_ELIGIBLE,
4.117 +
4.118 +      /// The Block Search pivot rule.
4.119 +      /// A specified number of arcs are examined in every iteration
4.120 +      /// in a wraparound fashion and the best eligible arc is selected
4.121 +      /// from this block.
4.122 +      BLOCK_SEARCH,
4.123 +
4.124 +      /// The Candidate List pivot rule.
4.125 +      /// In a major iteration a candidate list is built from eligible arcs
4.126 +      /// in a wraparound fashion and in the following minor iterations
4.127 +      /// the best eligible arc is selected from this list.
4.128 +      CANDIDATE_LIST,
4.129 +
4.130 +      /// The Altering Candidate List pivot rule.
4.131 +      /// It is a modified version of the Candidate List method.
4.132 +      /// It keeps only the several best eligible arcs from the former
4.133 +      /// candidate list and extends this list in every iteration.
4.134 +      ALTERING_LIST
4.135 +    };
4.136 +
4.137 +    /// \brief Enum type for selecting the problem type.
4.138 +    ///
4.139 +    /// Enum type for selecting the problem type, i.e. the direction of
4.140 +    /// the inequalities in the supply/demand constraints of the
4.141 +    /// \ref min_cost_flow "minimum cost flow problem".
4.142 +    ///
4.143 +    /// The default problem type is \c GEQ, since this form is supported
4.144 +    /// by other minimum cost flow algorithms and the \ref Circulation
4.145 +    /// algorithm as well.
4.146 +    /// The \c LEQ problem type can be selected using the \ref problemType()
4.147 +    /// function.
4.148 +    ///
4.149 +    /// Note that the equality form is a special case of both problem type.
4.150 +    enum ProblemType {
4.151 +
4.152 +      /// This option means that there are "<em>greater or equal</em>"
4.153 +      /// constraints in the defintion, i.e. the exact formulation of the
4.154 +      /// problem is the following.
4.155 +      /**
4.156 +          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
4.157 +          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
4.158 +              sup(u) \quad \forall u\in V \f]
4.159 +          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
4.160 +      */
4.161 +      /// It means that the total demand must be greater or equal to the
4.162 +      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
4.163 +      /// negative) and all the supplies have to be carried out from
4.164 +      /// the supply nodes, but there could be demands that are not
4.165 +      /// satisfied.
4.166 +      GEQ,
4.167 +      /// It is just an alias for the \c GEQ option.
4.168 +      CARRY_SUPPLIES = GEQ,
4.169 +
4.170 +      /// This option means that there are "<em>less or equal</em>"
4.171 +      /// constraints in the defintion, i.e. the exact formulation of the
4.172 +      /// problem is the following.
4.173 +      /**
4.174 +          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
4.175 +          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
4.176 +              sup(u) \quad \forall u\in V \f]
4.177 +          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
4.178 +      */
4.179 +      /// It means that the total demand must be less or equal to the
4.180 +      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
4.181 +      /// positive) and all the demands have to be satisfied, but there
4.182 +      /// could be supplies that are not carried out from the supply
4.183 +      /// nodes.
4.184 +      LEQ,
4.185 +      /// It is just an alias for the \c LEQ option.
4.186 +      SATISFY_DEMANDS = LEQ
4.187 +    };
4.188 +
4.189 +  private:
4.190 +
4.191 +    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
4.192 +
4.193 +    typedef typename GR::template ArcMap<Flow> FlowArcMap;
4.194 +    typedef typename GR::template ArcMap<Cost> CostArcMap;
4.195 +    typedef typename GR::template NodeMap<Flow> FlowNodeMap;
4.196 +
4.197 +    typedef std::vector<Arc> ArcVector;
4.198 +    typedef std::vector<Node> NodeVector;
4.199 +    typedef std::vector<int> IntVector;
4.200 +    typedef std::vector<bool> BoolVector;
4.201 +    typedef std::vector<Flow> FlowVector;
4.202 +    typedef std::vector<Cost> CostVector;
4.203 +
4.204 +    // State constants for arcs
4.205 +    enum ArcStateEnum {
4.206 +      STATE_UPPER = -1,
4.207 +      STATE_TREE  =  0,
4.208 +      STATE_LOWER =  1
4.209 +    };
4.210 +
4.211 +  private:
4.212 +
4.213 +    // Data related to the underlying digraph
4.214 +    const GR &_graph;
4.215 +    int _node_num;
4.216 +    int _arc_num;
4.217 +
4.218 +    // Parameters of the problem
4.219 +    FlowArcMap *_plower;
4.220 +    FlowArcMap *_pupper;
4.221 +    CostArcMap *_pcost;
4.222 +    FlowNodeMap *_psupply;
4.223 +    bool _pstsup;
4.224 +    Node _psource, _ptarget;
4.225 +    Flow _pstflow;
4.226 +    ProblemType _ptype;
4.227 +
4.228 +    // Result maps
4.229 +    FlowMap *_flow_map;
4.230 +    PotentialMap *_potential_map;
4.231 +    bool _local_flow;
4.232 +    bool _local_potential;
4.233 +
4.234 +    // Data structures for storing the digraph
4.235 +    IntNodeMap _node_id;
4.236 +    ArcVector _arc_ref;
4.237 +    IntVector _source;
4.238 +    IntVector _target;
4.239 +
4.240 +    // Node and arc data
4.241 +    FlowVector _cap;
4.242 +    CostVector _cost;
4.243 +    FlowVector _supply;
4.244 +    FlowVector _flow;
4.245 +    CostVector _pi;
4.246 +
4.247 +    // Data for storing the spanning tree structure
4.248 +    IntVector _parent;
4.249 +    IntVector _pred;
4.250 +    IntVector _thread;
4.251 +    IntVector _rev_thread;
4.252 +    IntVector _succ_num;
4.253 +    IntVector _last_succ;
4.254 +    IntVector _dirty_revs;
4.255 +    BoolVector _forward;
4.256 +    IntVector _state;
4.257 +    int _root;
4.258 +
4.259 +    // Temporary data used in the current pivot iteration
4.260 +    int in_arc, join, u_in, v_in, u_out, v_out;
4.261 +    int first, second, right, last;
4.262 +    int stem, par_stem, new_stem;
4.263 +    Flow delta;
4.264 +
4.265 +  private:
4.266 +
4.267 +    // Implementation of the First Eligible pivot rule
4.268 +    class FirstEligiblePivotRule
4.269 +    {
4.270 +    private:
4.271 +
4.272 +      // References to the NetworkSimplex class
4.273 +      const IntVector  &_source;
4.274 +      const IntVector  &_target;
4.275 +      const CostVector &_cost;
4.276 +      const IntVector  &_state;
4.277 +      const CostVector &_pi;
4.278 +      int &_in_arc;
4.279 +      int _arc_num;
4.280 +
4.281 +      // Pivot rule data
4.282 +      int _next_arc;
4.283 +
4.284 +    public:
4.285 +
4.286 +      // Constructor
4.287 +      FirstEligiblePivotRule(NetworkSimplex &ns) :
4.288 +        _source(ns._source), _target(ns._target),
4.289 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.290 +        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
4.291 +      {}
4.292 +
4.293 +      // Find next entering arc
4.294 +      bool findEnteringArc() {
4.295 +        Cost c;
4.296 +        for (int e = _next_arc; e < _arc_num; ++e) {
4.297 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.298 +          if (c < 0) {
4.299 +            _in_arc = e;
4.300 +            _next_arc = e + 1;
4.301 +            return true;
4.302 +          }
4.303 +        }
4.304 +        for (int e = 0; e < _next_arc; ++e) {
4.305 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.306 +          if (c < 0) {
4.307 +            _in_arc = e;
4.308 +            _next_arc = e + 1;
4.309 +            return true;
4.310 +          }
4.311 +        }
4.312 +        return false;
4.313 +      }
4.314 +
4.315 +    }; //class FirstEligiblePivotRule
4.316 +
4.317 +
4.318 +    // Implementation of the Best Eligible pivot rule
4.319 +    class BestEligiblePivotRule
4.320 +    {
4.321 +    private:
4.322 +
4.323 +      // References to the NetworkSimplex class
4.324 +      const IntVector  &_source;
4.325 +      const IntVector  &_target;
4.326 +      const CostVector &_cost;
4.327 +      const IntVector  &_state;
4.328 +      const CostVector &_pi;
4.329 +      int &_in_arc;
4.330 +      int _arc_num;
4.331 +
4.332 +    public:
4.333 +
4.334 +      // Constructor
4.335 +      BestEligiblePivotRule(NetworkSimplex &ns) :
4.336 +        _source(ns._source), _target(ns._target),
4.337 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.338 +        _in_arc(ns.in_arc), _arc_num(ns._arc_num)
4.339 +      {}
4.340 +
4.341 +      // Find next entering arc
4.342 +      bool findEnteringArc() {
4.343 +        Cost c, min = 0;
4.344 +        for (int e = 0; e < _arc_num; ++e) {
4.345 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.346 +          if (c < min) {
4.347 +            min = c;
4.348 +            _in_arc = e;
4.349 +          }
4.350 +        }
4.351 +        return min < 0;
4.352 +      }
4.353 +
4.354 +    }; //class BestEligiblePivotRule
4.355 +
4.356 +
4.357 +    // Implementation of the Block Search pivot rule
4.358 +    class BlockSearchPivotRule
4.359 +    {
4.360 +    private:
4.361 +
4.362 +      // References to the NetworkSimplex class
4.363 +      const IntVector  &_source;
4.364 +      const IntVector  &_target;
4.365 +      const CostVector &_cost;
4.366 +      const IntVector  &_state;
4.367 +      const CostVector &_pi;
4.368 +      int &_in_arc;
4.369 +      int _arc_num;
4.370 +
4.371 +      // Pivot rule data
4.372 +      int _block_size;
4.373 +      int _next_arc;
4.374 +
4.375 +    public:
4.376 +
4.377 +      // Constructor
4.378 +      BlockSearchPivotRule(NetworkSimplex &ns) :
4.379 +        _source(ns._source), _target(ns._target),
4.380 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.381 +        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
4.382 +      {
4.383 +        // The main parameters of the pivot rule
4.384 +        const double BLOCK_SIZE_FACTOR = 2.0;
4.385 +        const int MIN_BLOCK_SIZE = 10;
4.386 +
4.387 +        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
4.388 +                                MIN_BLOCK_SIZE );
4.389 +      }
4.390 +
4.391 +      // Find next entering arc
4.392 +      bool findEnteringArc() {
4.393 +        Cost c, min = 0;
4.394 +        int cnt = _block_size;
4.395 +        int e, min_arc = _next_arc;
4.396 +        for (e = _next_arc; e < _arc_num; ++e) {
4.397 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.398 +          if (c < min) {
4.399 +            min = c;
4.400 +            min_arc = e;
4.401 +          }
4.402 +          if (--cnt == 0) {
4.403 +            if (min < 0) break;
4.404 +            cnt = _block_size;
4.405 +          }
4.406 +        }
4.407 +        if (min == 0 || cnt > 0) {
4.408 +          for (e = 0; e < _next_arc; ++e) {
4.409 +            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.410 +            if (c < min) {
4.411 +              min = c;
4.412 +              min_arc = e;
4.413 +            }
4.414 +            if (--cnt == 0) {
4.415 +              if (min < 0) break;
4.416 +              cnt = _block_size;
4.417 +            }
4.418 +          }
4.419 +        }
4.420 +        if (min >= 0) return false;
4.421 +        _in_arc = min_arc;
4.422 +        _next_arc = e;
4.423 +        return true;
4.424 +      }
4.425 +
4.426 +    }; //class BlockSearchPivotRule
4.427 +
4.428 +
4.429 +    // Implementation of the Candidate List pivot rule
4.430 +    class CandidateListPivotRule
4.431 +    {
4.432 +    private:
4.433 +
4.434 +      // References to the NetworkSimplex class
4.435 +      const IntVector  &_source;
4.436 +      const IntVector  &_target;
4.437 +      const CostVector &_cost;
4.438 +      const IntVector  &_state;
4.439 +      const CostVector &_pi;
4.440 +      int &_in_arc;
4.441 +      int _arc_num;
4.442 +
4.443 +      // Pivot rule data
4.444 +      IntVector _candidates;
4.445 +      int _list_length, _minor_limit;
4.446 +      int _curr_length, _minor_count;
4.447 +      int _next_arc;
4.448 +
4.449 +    public:
4.450 +
4.451 +      /// Constructor
4.452 +      CandidateListPivotRule(NetworkSimplex &ns) :
4.453 +        _source(ns._source), _target(ns._target),
4.454 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.455 +        _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
4.456 +      {
4.457 +        // The main parameters of the pivot rule
4.458 +        const double LIST_LENGTH_FACTOR = 1.0;
4.459 +        const int MIN_LIST_LENGTH = 10;
4.460 +        const double MINOR_LIMIT_FACTOR = 0.1;
4.461 +        const int MIN_MINOR_LIMIT = 3;
4.462 +
4.463 +        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_arc_num)),
4.464 +                                 MIN_LIST_LENGTH );
4.465 +        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
4.466 +                                 MIN_MINOR_LIMIT );
4.467 +        _curr_length = _minor_count = 0;
4.468 +        _candidates.resize(_list_length);
4.469 +      }
4.470 +
4.471 +      /// Find next entering arc
4.472 +      bool findEnteringArc() {
4.473 +        Cost min, c;
4.474 +        int e, min_arc = _next_arc;
4.475 +        if (_curr_length > 0 && _minor_count < _minor_limit) {
4.476 +          // Minor iteration: select the best eligible arc from the
4.477 +          // current candidate list
4.478 +          ++_minor_count;
4.479 +          min = 0;
4.480 +          for (int i = 0; i < _curr_length; ++i) {
4.481 +            e = _candidates[i];
4.482 +            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.483 +            if (c < min) {
4.484 +              min = c;
4.485 +              min_arc = e;
4.486 +            }
4.487 +            if (c >= 0) {
4.488 +              _candidates[i--] = _candidates[--_curr_length];
4.489 +            }
4.490 +          }
4.491 +          if (min < 0) {
4.492 +            _in_arc = min_arc;
4.493 +            return true;
4.494 +          }
4.495 +        }
4.496 +
4.497 +        // Major iteration: build a new candidate list
4.498 +        min = 0;
4.499 +        _curr_length = 0;
4.500 +        for (e = _next_arc; e < _arc_num; ++e) {
4.501 +          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.502 +          if (c < 0) {
4.503 +            _candidates[_curr_length++] = e;
4.504 +            if (c < min) {
4.505 +              min = c;
4.506 +              min_arc = e;
4.507 +            }
4.508 +            if (_curr_length == _list_length) break;
4.509 +          }
4.510 +        }
4.511 +        if (_curr_length < _list_length) {
4.512 +          for (e = 0; e < _next_arc; ++e) {
4.513 +            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.514 +            if (c < 0) {
4.515 +              _candidates[_curr_length++] = e;
4.516 +              if (c < min) {
4.517 +                min = c;
4.518 +                min_arc = e;
4.519 +              }
4.520 +              if (_curr_length == _list_length) break;
4.521 +            }
4.522 +          }
4.523 +        }
4.524 +        if (_curr_length == 0) return false;
4.525 +        _minor_count = 1;
4.526 +        _in_arc = min_arc;
4.527 +        _next_arc = e;
4.528 +        return true;
4.529 +      }
4.530 +
4.531 +    }; //class CandidateListPivotRule
4.532 +
4.533 +
4.534 +    // Implementation of the Altering Candidate List pivot rule
4.535 +    class AlteringListPivotRule
4.536 +    {
4.537 +    private:
4.538 +
4.539 +      // References to the NetworkSimplex class
4.540 +      const IntVector  &_source;
4.541 +      const IntVector  &_target;
4.542 +      const CostVector &_cost;
4.543 +      const IntVector  &_state;
4.544 +      const CostVector &_pi;
4.545 +      int &_in_arc;
4.546 +      int _arc_num;
4.547 +
4.548 +      // Pivot rule data
4.549 +      int _block_size, _head_length, _curr_length;
4.550 +      int _next_arc;
4.551 +      IntVector _candidates;
4.552 +      CostVector _cand_cost;
4.553 +
4.554 +      // Functor class to compare arcs during sort of the candidate list
4.555 +      class SortFunc
4.556 +      {
4.557 +      private:
4.558 +        const CostVector &_map;
4.559 +      public:
4.560 +        SortFunc(const CostVector &map) : _map(map) {}
4.561 +        bool operator()(int left, int right) {
4.562 +          return _map[left] > _map[right];
4.563 +        }
4.564 +      };
4.565 +
4.566 +      SortFunc _sort_func;
4.567 +
4.568 +    public:
4.569 +
4.570 +      // Constructor
4.571 +      AlteringListPivotRule(NetworkSimplex &ns) :
4.572 +        _source(ns._source), _target(ns._target),
4.573 +        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.574 +        _in_arc(ns.in_arc), _arc_num(ns._arc_num),
4.575 +        _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
4.576 +      {
4.577 +        // The main parameters of the pivot rule
4.578 +        const double BLOCK_SIZE_FACTOR = 1.5;
4.579 +        const int MIN_BLOCK_SIZE = 10;
4.580 +        const double HEAD_LENGTH_FACTOR = 0.1;
4.581 +        const int MIN_HEAD_LENGTH = 3;
4.582 +
4.583 +        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
4.584 +                                MIN_BLOCK_SIZE );
4.585 +        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
4.586 +                                 MIN_HEAD_LENGTH );
4.587 +        _candidates.resize(_head_length + _block_size);
4.588 +        _curr_length = 0;
4.589 +      }
4.590 +
4.591 +      // Find next entering arc
4.592 +      bool findEnteringArc() {
4.593 +        // Check the current candidate list
4.594 +        int e;
4.595 +        for (int i = 0; i < _curr_length; ++i) {
4.596 +          e = _candidates[i];
4.597 +          _cand_cost[e] = _state[e] *
4.598 +            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.599 +          if (_cand_cost[e] >= 0) {
4.600 +            _candidates[i--] = _candidates[--_curr_length];
4.601 +          }
4.602 +        }
4.603 +
4.604 +        // Extend the list
4.605 +        int cnt = _block_size;
4.606 +        int last_arc = 0;
4.607 +        int limit = _head_length;
4.608 +
4.609 +        for (int e = _next_arc; e < _arc_num; ++e) {
4.610 +          _cand_cost[e] = _state[e] *
4.611 +            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.612 +          if (_cand_cost[e] < 0) {
4.613 +            _candidates[_curr_length++] = e;
4.614 +            last_arc = e;
4.615 +          }
4.616 +          if (--cnt == 0) {
4.617 +            if (_curr_length > limit) break;
4.618 +            limit = 0;
4.619 +            cnt = _block_size;
4.620 +          }
4.621 +        }
4.622 +        if (_curr_length <= limit) {
4.623 +          for (int e = 0; e < _next_arc; ++e) {
4.624 +            _cand_cost[e] = _state[e] *
4.625 +              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.626 +            if (_cand_cost[e] < 0) {
4.627 +              _candidates[_curr_length++] = e;
4.628 +              last_arc = e;
4.629 +            }
4.630 +            if (--cnt == 0) {
4.631 +              if (_curr_length > limit) break;
4.632 +              limit = 0;
4.633 +              cnt = _block_size;
4.634 +            }
4.635 +          }
4.636 +        }
4.637 +        if (_curr_length == 0) return false;
4.638 +        _next_arc = last_arc + 1;
4.639 +
4.640 +        // Make heap of the candidate list (approximating a partial sort)
4.641 +        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
4.642 +                   _sort_func );
4.643 +
4.644 +        // Pop the first element of the heap
4.645 +        _in_arc = _candidates[0];
4.646 +        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
4.647 +                  _sort_func );
4.648 +        _curr_length = std::min(_head_length, _curr_length - 1);
4.649 +        return true;
4.650 +      }
4.651 +
4.652 +    }; //class AlteringListPivotRule
4.653 +
4.654 +  public:
4.655 +
4.656 +    /// \brief Constructor.
4.657 +    ///
4.658 +    /// The constructor of the class.
4.659 +    ///
4.660 +    /// \param graph The digraph the algorithm runs on.
4.661 +    NetworkSimplex(const GR& graph) :
4.662 +      _graph(graph),
4.663 +      _plower(NULL), _pupper(NULL), _pcost(NULL),
4.664 +      _psupply(NULL), _pstsup(false), _ptype(GEQ),
4.665 +      _flow_map(NULL), _potential_map(NULL),
4.666 +      _local_flow(false), _local_potential(false),
4.667 +      _node_id(graph)
4.668 +    {
4.669 +      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
4.670 +                   std::numeric_limits<Flow>::is_signed,
4.671 +        "The flow type of NetworkSimplex must be signed integer");
4.672 +      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
4.673 +                   std::numeric_limits<Cost>::is_signed,
4.674 +        "The cost type of NetworkSimplex must be signed integer");
4.675 +    }
4.676 +
4.677 +    /// Destructor.
4.678 +    ~NetworkSimplex() {
4.679 +      if (_local_flow) delete _flow_map;
4.680 +      if (_local_potential) delete _potential_map;
4.681 +    }
4.682 +
4.683 +    /// \name Parameters
4.684 +    /// The parameters of the algorithm can be specified using these
4.685 +    /// functions.
4.686 +
4.687 +    /// @{
4.688 +
4.689 +    /// \brief Set the lower bounds on the arcs.
4.690 +    ///
4.691 +    /// This function sets the lower bounds on the arcs.
4.692 +    /// If neither this function nor \ref boundMaps() is used before
4.693 +    /// calling \ref run(), the lower bounds will be set to zero
4.694 +    /// on all arcs.
4.695 +    ///
4.696 +    /// \param map An arc map storing the lower bounds.
4.697 +    /// Its \c Value type must be convertible to the \c Flow type
4.698 +    /// of the algorithm.
4.699 +    ///
4.700 +    /// \return <tt>(*this)</tt>
4.701 +    template <typename LOWER>
4.702 +    NetworkSimplex& lowerMap(const LOWER& map) {
4.703 +      delete _plower;
4.704 +      _plower = new FlowArcMap(_graph);
4.705 +      for (ArcIt a(_graph); a != INVALID; ++a) {
4.706 +        (*_plower)[a] = map[a];
4.707 +      }
4.708 +      return *this;
4.709 +    }
4.710 +
4.711 +    /// \brief Set the upper bounds (capacities) on the arcs.
4.712 +    ///
4.713 +    /// This function sets the upper bounds (capacities) on the arcs.
4.714 +    /// If none of the functions \ref upperMap(), \ref capacityMap()
4.715 +    /// and \ref boundMaps() is used before calling \ref run(),
4.716 +    /// the upper bounds (capacities) will be set to
4.717 +    /// \c std::numeric_limits<Flow>::max() on all arcs.
4.718 +    ///
4.719 +    /// \param map An arc map storing the upper bounds.
4.720 +    /// Its \c Value type must be convertible to the \c Flow type
4.721 +    /// of the algorithm.
4.722 +    ///
4.723 +    /// \return <tt>(*this)</tt>
4.724 +    template<typename UPPER>
4.725 +    NetworkSimplex& upperMap(const UPPER& map) {
4.726 +      delete _pupper;
4.727 +      _pupper = new FlowArcMap(_graph);
4.728 +      for (ArcIt a(_graph); a != INVALID; ++a) {
4.729 +        (*_pupper)[a] = map[a];
4.730 +      }
4.731 +      return *this;
4.732 +    }
4.733 +
4.734 +    /// \brief Set the upper bounds (capacities) on the arcs.
4.735 +    ///
4.736 +    /// This function sets the upper bounds (capacities) on the arcs.
4.737 +    /// It is just an alias for \ref upperMap().
4.738 +    ///
4.739 +    /// \return <tt>(*this)</tt>
4.740 +    template<typename CAP>
4.741 +    NetworkSimplex& capacityMap(const CAP& map) {
4.742 +      return upperMap(map);
4.743 +    }
4.744 +
4.745 +    /// \brief Set the lower and upper bounds on the arcs.
4.746 +    ///
4.747 +    /// This function sets the lower and upper bounds on the arcs.
4.748 +    /// If neither this function nor \ref lowerMap() is used before
4.749 +    /// calling \ref run(), the lower bounds will be set to zero
4.750 +    /// on all arcs.
4.751 +    /// If none of the functions \ref upperMap(), \ref capacityMap()
4.752 +    /// and \ref boundMaps() is used before calling \ref run(),
4.753 +    /// the upper bounds (capacities) will be set to
4.754 +    /// \c std::numeric_limits<Flow>::max() on all arcs.
4.755 +    ///
4.756 +    /// \param lower An arc map storing the lower bounds.
4.757 +    /// \param upper An arc map storing the upper bounds.
4.758 +    ///
4.759 +    /// The \c Value type of the maps must be convertible to the
4.760 +    /// \c Flow type of the algorithm.
4.761 +    ///
4.762 +    /// \note This function is just a shortcut of calling \ref lowerMap()
4.763 +    /// and \ref upperMap() separately.
4.764 +    ///
4.765 +    /// \return <tt>(*this)</tt>
4.766 +    template <typename LOWER, typename UPPER>
4.767 +    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
4.768 +      return lowerMap(lower).upperMap(upper);
4.769 +    }
4.770 +
4.771 +    /// \brief Set the costs of the arcs.
4.772 +    ///
4.773 +    /// This function sets the costs of the arcs.
4.774 +    /// If it is not used before calling \ref run(), the costs
4.775 +    /// will be set to \c 1 on all arcs.
4.776 +    ///
4.777 +    /// \param map An arc map storing the costs.
4.778 +    /// Its \c Value type must be convertible to the \c Cost type
4.779 +    /// of the algorithm.
4.780 +    ///
4.781 +    /// \return <tt>(*this)</tt>
4.782 +    template<typename COST>
4.783 +    NetworkSimplex& costMap(const COST& map) {
4.784 +      delete _pcost;
4.785 +      _pcost = new CostArcMap(_graph);
4.786 +      for (ArcIt a(_graph); a != INVALID; ++a) {
4.787 +        (*_pcost)[a] = map[a];
4.788 +      }
4.789 +      return *this;
4.790 +    }
4.791 +
4.792 +    /// \brief Set the supply values of the nodes.
4.793 +    ///
4.794 +    /// This function sets the supply values of the nodes.
4.795 +    /// If neither this function nor \ref stSupply() is used before
4.796 +    /// calling \ref run(), the supply of each node will be set to zero.
4.797 +    /// (It makes sense only if non-zero lower bounds are given.)
4.798 +    ///
4.799 +    /// \param map A node map storing the supply values.
4.800 +    /// Its \c Value type must be convertible to the \c Flow type
4.801 +    /// of the algorithm.
4.802 +    ///
4.803 +    /// \return <tt>(*this)</tt>
4.804 +    template<typename SUP>
4.805 +    NetworkSimplex& supplyMap(const SUP& map) {
4.806 +      delete _psupply;
4.807 +      _pstsup = false;
4.808 +      _psupply = new FlowNodeMap(_graph);
4.809 +      for (NodeIt n(_graph); n != INVALID; ++n) {
4.810 +        (*_psupply)[n] = map[n];
4.811 +      }
4.812 +      return *this;
4.813 +    }
4.814 +
4.815 +    /// \brief Set single source and target nodes and a supply value.
4.816 +    ///
4.817 +    /// This function sets a single source node and a single target node
4.818 +    /// and the required flow value.
4.819 +    /// If neither this function nor \ref supplyMap() is used before
4.820 +    /// calling \ref run(), the supply of each node will be set to zero.
4.821 +    /// (It makes sense only if non-zero lower bounds are given.)
4.822 +    ///
4.823 +    /// \param s The source node.
4.824 +    /// \param t The target node.
4.825 +    /// \param k The required amount of flow from node \c s to node \c t
4.826 +    /// (i.e. the supply of \c s and the demand of \c t).
4.827 +    ///
4.828 +    /// \return <tt>(*this)</tt>
4.829 +    NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
4.830 +      delete _psupply;
4.831 +      _psupply = NULL;
4.832 +      _pstsup = true;
4.833 +      _psource = s;
4.834 +      _ptarget = t;
4.835 +      _pstflow = k;
4.836 +      return *this;
4.837 +    }
4.838 +
4.839 +    /// \brief Set the problem type.
4.840 +    ///
4.841 +    /// This function sets the problem type for the algorithm.
4.842 +    /// If it is not used before calling \ref run(), the \ref GEQ problem
4.843 +    /// type will be used.
4.844 +    ///
4.845 +    /// For more information see \ref ProblemType.
4.846 +    ///
4.847 +    /// \return <tt>(*this)</tt>
4.848 +    NetworkSimplex& problemType(ProblemType problem_type) {
4.849 +      _ptype = problem_type;
4.850 +      return *this;
4.851 +    }
4.852 +
4.853 +    /// \brief Set the flow map.
4.854 +    ///
4.855 +    /// This function sets the flow map.
4.856 +    /// If it is not used before calling \ref run(), an instance will
4.857 +    /// be allocated automatically. The destructor deallocates this
4.858 +    /// automatically allocated map, of course.
4.859 +    ///
4.860 +    /// \return <tt>(*this)</tt>
4.861 +    NetworkSimplex& flowMap(FlowMap& map) {
4.862 +      if (_local_flow) {
4.863 +        delete _flow_map;
4.864 +        _local_flow = false;
4.865 +      }
4.866 +      _flow_map = &map;
4.867 +      return *this;
4.868 +    }
4.869 +
4.870 +    /// \brief Set the potential map.
4.871 +    ///
4.872 +    /// This function sets the potential map, which is used for storing
4.873 +    /// the dual solution.
4.874 +    /// If it is not used before calling \ref run(), an instance will
4.875 +    /// be allocated automatically. The destructor deallocates this
4.876 +    /// automatically allocated map, of course.
4.877 +    ///
4.878 +    /// \return <tt>(*this)</tt>
4.879 +    NetworkSimplex& potentialMap(PotentialMap& map) {
4.880 +      if (_local_potential) {
4.881 +        delete _potential_map;
4.882 +        _local_potential = false;
4.883 +      }
4.884 +      _potential_map = &map;
4.885 +      return *this;
4.886 +    }
4.887 +
4.888 +    /// @}
4.889 +
4.890 +    /// \name Execution Control
4.891 +    /// The algorithm can be executed using \ref run().
4.892 +
4.893 +    /// @{
4.894 +
4.895 +    /// \brief Run the algorithm.
4.896 +    ///
4.897 +    /// This function runs the algorithm.
4.898 +    /// The paramters can be specified using functions \ref lowerMap(),
4.899 +    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
4.900 +    /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
4.901 +    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
4.902 +    /// For example,
4.903 +    /// \code
4.904 +    ///   NetworkSimplex<ListDigraph> ns(graph);
4.905 +    ///   ns.boundMaps(lower, upper).costMap(cost)
4.906 +    ///     .supplyMap(sup).run();
4.907 +    /// \endcode
4.908 +    ///
4.909 +    /// This function can be called more than once. All the parameters
4.910 +    /// that have been given are kept for the next call, unless
4.911 +    /// \ref reset() is called, thus only the modified parameters
4.912 +    /// have to be set again. See \ref reset() for examples.
4.913 +    ///
4.914 +    /// \param pivot_rule The pivot rule that will be used during the
4.915 +    /// algorithm. For more information see \ref PivotRule.
4.916 +    ///
4.917 +    /// \return \c true if a feasible flow can be found.
4.918 +    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
4.919 +      return init() && start(pivot_rule);
4.920 +    }
4.921 +
4.922 +    /// \brief Reset all the parameters that have been given before.
4.923 +    ///
4.924 +    /// This function resets all the paramaters that have been given
4.925 +    /// before using functions \ref lowerMap(), \ref upperMap(),
4.926 +    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
4.927 +    /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
4.928 +    /// \ref flowMap() and \ref potentialMap().
4.929 +    ///
4.930 +    /// It is useful for multiple run() calls. If this function is not
4.931 +    /// used, all the parameters given before are kept for the next
4.932 +    /// \ref run() call.
4.933 +    ///
4.934 +    /// For example,
4.935 +    /// \code
4.936 +    ///   NetworkSimplex<ListDigraph> ns(graph);
4.937 +    ///
4.938 +    ///   // First run
4.939 +    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
4.940 +    ///     .supplyMap(sup).run();
4.941 +    ///
4.942 +    ///   // Run again with modified cost map (reset() is not called,
4.943 +    ///   // so only the cost map have to be set again)
4.944 +    ///   cost[e] += 100;
4.945 +    ///   ns.costMap(cost).run();
4.946 +    ///
4.947 +    ///   // Run again from scratch using reset()
4.948 +    ///   // (the lower bounds will be set to zero on all arcs)
4.949 +    ///   ns.reset();
4.950 +    ///   ns.capacityMap(cap).costMap(cost)
4.951 +    ///     .supplyMap(sup).run();
4.952 +    /// \endcode
4.953 +    ///
4.954 +    /// \return <tt>(*this)</tt>
4.955 +    NetworkSimplex& reset() {
4.956 +      delete _plower;
4.957 +      delete _pupper;
4.958 +      delete _pcost;
4.959 +      delete _psupply;
4.960 +      _plower = NULL;
4.961 +      _pupper = NULL;
4.962 +      _pcost = NULL;
4.963 +      _psupply = NULL;
4.964 +      _pstsup = false;
4.965 +      _ptype = GEQ;
4.966 +      if (_local_flow) delete _flow_map;
4.967 +      if (_local_potential) delete _potential_map;
4.968 +      _flow_map = NULL;
4.969 +      _potential_map = NULL;
4.970 +      _local_flow = false;
4.971 +      _local_potential = false;
4.972 +
4.973 +      return *this;
4.974 +    }
4.975 +
4.976 +    /// @}
4.977 +
4.978 +    /// \name Query Functions
4.979 +    /// The results of the algorithm can be obtained using these
4.980 +    /// functions.\n
4.981 +    /// The \ref run() function must be called before using them.
4.982 +
4.983 +    /// @{
4.984 +
4.985 +    /// \brief Return the total cost of the found flow.
4.986 +    ///
4.987 +    /// This function returns the total cost of the found flow.
4.988 +    /// The complexity of the function is O(e).
4.989 +    ///
4.990 +    /// \note The return type of the function can be specified as a
4.991 +    /// template parameter. For example,
4.992 +    /// \code
4.993 +    ///   ns.totalCost<double>();
4.994 +    /// \endcode
4.995 +    /// It is useful if the total cost cannot be stored in the \c Cost
4.996 +    /// type of the algorithm, which is the default return type of the
4.997 +    /// function.
4.998 +    ///
4.999 +    /// \pre \ref run() must be called before using this function.
4.1000 +    template <typename Num>
4.1001 +    Num totalCost() const {
4.1002 +      Num c = 0;
4.1003 +      if (_pcost) {
4.1004 +        for (ArcIt e(_graph); e != INVALID; ++e)
4.1005 +          c += (*_flow_map)[e] * (*_pcost)[e];
4.1006 +      } else {
4.1007 +        for (ArcIt e(_graph); e != INVALID; ++e)
4.1008 +          c += (*_flow_map)[e];
4.1009 +      }
4.1010 +      return c;
4.1011 +    }
4.1012 +
4.1013 +#ifndef DOXYGEN
4.1014 +    Cost totalCost() const {
4.1016 +    }
4.1017 +#endif
4.1018 +
4.1019 +    /// \brief Return the flow on the given arc.
4.1020 +    ///
4.1021 +    /// This function returns the flow on the given arc.
4.1022 +    ///
4.1023 +    /// \pre \ref run() must be called before using this function.
4.1024 +    Flow flow(const Arc& a) const {
4.1025 +      return (*_flow_map)[a];
4.1026 +    }
4.1027 +
4.1028 +    /// \brief Return a const reference to the flow map.
4.1029 +    ///
4.1030 +    /// This function returns a const reference to an arc map storing
4.1031 +    /// the found flow.
4.1032 +    ///
4.1033 +    /// \pre \ref run() must be called before using this function.
4.1034 +    const FlowMap& flowMap() const {
4.1035 +      return *_flow_map;
4.1036 +    }
4.1037 +
4.1038 +    /// \brief Return the potential (dual value) of the given node.
4.1039 +    ///
4.1040 +    /// This function returns the potential (dual value) of the
4.1041 +    /// given node.
4.1042 +    ///
4.1043 +    /// \pre \ref run() must be called before using this function.
4.1044 +    Cost potential(const Node& n) const {
4.1045 +      return (*_potential_map)[n];
4.1046 +    }
4.1047 +
4.1048 +    /// \brief Return a const reference to the potential map
4.1049 +    /// (the dual solution).
4.1050 +    ///
4.1051 +    /// This function returns a const reference to a node map storing
4.1052 +    /// the found potentials, which form the dual solution of the
4.1053 +    /// \ref min_cost_flow "minimum cost flow" problem.
4.1054 +    ///
4.1055 +    /// \pre \ref run() must be called before using this function.
4.1056 +    const PotentialMap& potentialMap() const {
4.1057 +      return *_potential_map;
4.1058 +    }
4.1059 +
4.1060 +    /// @}
4.1061 +
4.1062 +  private:
4.1063 +
4.1064 +    // Initialize internal data structures
4.1065 +    bool init() {
4.1066 +      // Initialize result maps
4.1067 +      if (!_flow_map) {
4.1068 +        _flow_map = new FlowMap(_graph);
4.1069 +        _local_flow = true;
4.1070 +      }
4.1071 +      if (!_potential_map) {
4.1072 +        _potential_map = new PotentialMap(_graph);
4.1073 +        _local_potential = true;
4.1074 +      }
4.1075 +
4.1076 +      // Initialize vectors
4.1077 +      _node_num = countNodes(_graph);
4.1078 +      _arc_num = countArcs(_graph);
4.1079 +      int all_node_num = _node_num + 1;
4.1080 +      int all_arc_num = _arc_num + _node_num;
4.1081 +      if (_node_num == 0) return false;
4.1082 +
4.1083 +      _arc_ref.resize(_arc_num);
4.1084 +      _source.resize(all_arc_num);
4.1085 +      _target.resize(all_arc_num);
4.1086 +
4.1087 +      _cap.resize(all_arc_num);
4.1088 +      _cost.resize(all_arc_num);
4.1089 +      _supply.resize(all_node_num);
4.1090 +      _flow.resize(all_arc_num);
4.1091 +      _pi.resize(all_node_num);
4.1092 +
4.1093 +      _parent.resize(all_node_num);
4.1094 +      _pred.resize(all_node_num);
4.1095 +      _forward.resize(all_node_num);
4.1098 +      _succ_num.resize(all_node_num);
4.1099 +      _last_succ.resize(all_node_num);
4.1100 +      _state.resize(all_arc_num);
4.1101 +
4.1102 +      // Initialize node related data
4.1103 +      bool valid_supply = true;
4.1104 +      Flow sum_supply = 0;
4.1105 +      if (!_pstsup && !_psupply) {
4.1106 +        _pstsup = true;
4.1107 +        _psource = _ptarget = NodeIt(_graph);
4.1108 +        _pstflow = 0;
4.1109 +      }
4.1110 +      if (_psupply) {
4.1111 +        int i = 0;
4.1112 +        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
4.1113 +          _node_id[n] = i;
4.1114 +          _supply[i] = (*_psupply)[n];
4.1115 +          sum_supply += _supply[i];
4.1116 +        }
4.1117 +        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
4.1118 +                       (_ptype == LEQ && sum_supply >= 0);
4.1119 +      } else {
4.1120 +        int i = 0;
4.1121 +        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
4.1122 +          _node_id[n] = i;
4.1123 +          _supply[i] = 0;
4.1124 +        }
4.1125 +        _supply[_node_id[_psource]] =  _pstflow;
4.1126 +        _supply[_node_id[_ptarget]] = -_pstflow;
4.1127 +      }
4.1128 +      if (!valid_supply) return false;
4.1129 +
4.1130 +      // Infinite capacity value
4.1131 +      Flow inf_cap =
4.1132 +        std::numeric_limits<Flow>::has_infinity ?
4.1133 +        std::numeric_limits<Flow>::infinity() :
4.1134 +        std::numeric_limits<Flow>::max();
4.1135 +
4.1136 +      // Initialize artifical cost
4.1137 +      Cost art_cost;
4.1138 +      if (std::numeric_limits<Cost>::is_exact) {
4.1139 +        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
4.1140 +      } else {
4.1141 +        art_cost = std::numeric_limits<Cost>::min();
4.1142 +        for (int i = 0; i != _arc_num; ++i) {
4.1143 +          if (_cost[i] > art_cost) art_cost = _cost[i];
4.1144 +        }
4.1145 +        art_cost = (art_cost + 1) * _node_num;
4.1146 +      }
4.1147 +
4.1148 +      // Run Circulation to check if a feasible solution exists
4.1149 +      typedef ConstMap<Arc, Flow> ConstArcMap;
4.1150 +      FlowNodeMap *csup = NULL;
4.1151 +      bool local_csup = false;
4.1152 +      if (_psupply) {
4.1153 +        csup = _psupply;
4.1154 +      } else {
4.1155 +        csup = new FlowNodeMap(_graph, 0);
4.1156 +        (*csup)[_psource] =  _pstflow;
4.1157 +        (*csup)[_ptarget] = -_pstflow;
4.1158 +        local_csup = true;
4.1159 +      }
4.1160 +      bool circ_result = false;
4.1161 +      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
4.1162 +        // GEQ problem type
4.1163 +        if (_plower) {
4.1164 +          if (_pupper) {
4.1165 +            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
4.1166 +              circ(_graph, *_plower, *_pupper, *csup);
4.1167 +            circ_result = circ.run();
4.1168 +          } else {
4.1169 +            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
4.1170 +              circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
4.1171 +            circ_result = circ.run();
4.1172 +          }
4.1173 +        } else {
4.1174 +          if (_pupper) {
4.1175 +            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
4.1176 +              circ(_graph, ConstArcMap(0), *_pupper, *csup);
4.1177 +            circ_result = circ.run();
4.1178 +          } else {
4.1179 +            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
4.1180 +              circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
4.1181 +            circ_result = circ.run();
4.1182 +          }
4.1183 +        }
4.1184 +      } else {
4.1185 +        // LEQ problem type
4.1186 +        typedef ReverseDigraph<const GR> RevGraph;
4.1187 +        typedef NegMap<FlowNodeMap> NegNodeMap;
4.1188 +        RevGraph rgraph(_graph);
4.1189 +        NegNodeMap neg_csup(*csup);
4.1190 +        if (_plower) {
4.1191 +          if (_pupper) {
4.1192 +            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
4.1193 +              circ(rgraph, *_plower, *_pupper, neg_csup);
4.1194 +            circ_result = circ.run();
4.1195 +          } else {
4.1196 +            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
4.1197 +              circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
4.1198 +            circ_result = circ.run();
4.1199 +          }
4.1200 +        } else {
4.1201 +          if (_pupper) {
4.1202 +            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
4.1203 +              circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
4.1204 +            circ_result = circ.run();
4.1205 +          } else {
4.1206 +            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
4.1207 +              circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
4.1208 +            circ_result = circ.run();
4.1209 +          }
4.1210 +        }
4.1211 +      }
4.1212 +      if (local_csup) delete csup;
4.1213 +      if (!circ_result) return false;
4.1214 +
4.1215 +      // Set data for the artificial root node
4.1216 +      _root = _node_num;
4.1217 +      _parent[_root] = -1;
4.1218 +      _pred[_root] = -1;
4.1219 +      _thread[_root] = 0;
4.1220 +      _rev_thread[0] = _root;
4.1221 +      _succ_num[_root] = all_node_num;
4.1222 +      _last_succ[_root] = _root - 1;
4.1223 +      _supply[_root] = -sum_supply;
4.1224 +      if (sum_supply < 0) {
4.1225 +        _pi[_root] = -art_cost;
4.1226 +      } else {
4.1227 +        _pi[_root] = art_cost;
4.1228 +      }
4.1229 +
4.1230 +      // Store the arcs in a mixed order
4.1231 +      int k = std::max(int(sqrt(_arc_num)), 10);
4.1232 +      int i = 0;
4.1233 +      for (ArcIt e(_graph); e != INVALID; ++e) {
4.1234 +        _arc_ref[i] = e;
4.1235 +        if ((i += k) >= _arc_num) i = (i % k) + 1;
4.1236 +      }
4.1237 +
4.1238 +      // Initialize arc maps
4.1239 +      if (_pupper && _pcost) {
4.1240 +        for (int i = 0; i != _arc_num; ++i) {
4.1241 +          Arc e = _arc_ref[i];
4.1242 +          _source[i] = _node_id[_graph.source(e)];
4.1243 +          _target[i] = _node_id[_graph.target(e)];
4.1244 +          _cap[i] = (*_pupper)[e];
4.1245 +          _cost[i] = (*_pcost)[e];
4.1246 +          _flow[i] = 0;
4.1247 +          _state[i] = STATE_LOWER;
4.1248 +        }
4.1249 +      } else {
4.1250 +        for (int i = 0; i != _arc_num; ++i) {
4.1251 +          Arc e = _arc_ref[i];
4.1252 +          _source[i] = _node_id[_graph.source(e)];
4.1253 +          _target[i] = _node_id[_graph.target(e)];
4.1254 +          _flow[i] = 0;
4.1255 +          _state[i] = STATE_LOWER;
4.1256 +        }
4.1257 +        if (_pupper) {
4.1258 +          for (int i = 0; i != _arc_num; ++i)
4.1259 +            _cap[i] = (*_pupper)[_arc_ref[i]];
4.1260 +        } else {
4.1261 +          for (int i = 0; i != _arc_num; ++i)
4.1262 +            _cap[i] = inf_cap;
4.1263 +        }
4.1264 +        if (_pcost) {
4.1265 +          for (int i = 0; i != _arc_num; ++i)
4.1266 +            _cost[i] = (*_pcost)[_arc_ref[i]];
4.1267 +        } else {
4.1268 +          for (int i = 0; i != _arc_num; ++i)
4.1269 +            _cost[i] = 1;
4.1270 +        }
4.1271 +      }
4.1272 +
4.1273 +      // Remove non-zero lower bounds
4.1274 +      if (_plower) {
4.1275 +        for (int i = 0; i != _arc_num; ++i) {
4.1276 +          Flow c = (*_plower)[_arc_ref[i]];
4.1277 +          if (c != 0) {
4.1278 +            _cap[i] -= c;
4.1279 +            _supply[_source[i]] -= c;
4.1280 +            _supply[_target[i]] += c;
4.1281 +          }
4.1282 +        }
4.1283 +      }
4.1284 +
4.1285 +      // Add artificial arcs and initialize the spanning tree data structure
4.1286 +      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
4.1287 +        _thread[u] = u + 1;
4.1288 +        _rev_thread[u + 1] = u;
4.1289 +        _succ_num[u] = 1;
4.1290 +        _last_succ[u] = u;
4.1291 +        _parent[u] = _root;
4.1292 +        _pred[u] = e;
4.1293 +        _cost[e] = art_cost;
4.1294 +        _cap[e] = inf_cap;
4.1295 +        _state[e] = STATE_TREE;
4.1296 +        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
4.1297 +          _flow[e] = _supply[u];
4.1298 +          _forward[u] = true;
4.1299 +          _pi[u] = -art_cost + _pi[_root];
4.1300 +        } else {
4.1301 +          _flow[e] = -_supply[u];
4.1302 +          _forward[u] = false;
4.1303 +          _pi[u] = art_cost + _pi[_root];
4.1304 +        }
4.1305 +      }
4.1306 +
4.1307 +      return true;
4.1308 +    }
4.1309 +
4.1310 +    // Find the join node
4.1311 +    void findJoinNode() {
4.1312 +      int u = _source[in_arc];
4.1313 +      int v = _target[in_arc];
4.1314 +      while (u != v) {
4.1315 +        if (_succ_num[u] < _succ_num[v]) {
4.1316 +          u = _parent[u];
4.1317 +        } else {
4.1318 +          v = _parent[v];
4.1319 +        }
4.1320 +      }
4.1321 +      join = u;
4.1322 +    }
4.1323 +
4.1324 +    // Find the leaving arc of the cycle and returns true if the
4.1325 +    // leaving arc is not the same as the entering arc
4.1326 +    bool findLeavingArc() {
4.1327 +      // Initialize first and second nodes according to the direction
4.1328 +      // of the cycle
4.1329 +      if (_state[in_arc] == STATE_LOWER) {
4.1330 +        first  = _source[in_arc];
4.1331 +        second = _target[in_arc];
4.1332 +      } else {
4.1333 +        first  = _target[in_arc];
4.1334 +        second = _source[in_arc];
4.1335 +      }
4.1336 +      delta = _cap[in_arc];
4.1337 +      int result = 0;
4.1338 +      Flow d;
4.1339 +      int e;
4.1340 +
4.1341 +      // Search the cycle along the path form the first node to the root
4.1342 +      for (int u = first; u != join; u = _parent[u]) {
4.1343 +        e = _pred[u];
4.1344 +        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
4.1345 +        if (d < delta) {
4.1346 +          delta = d;
4.1347 +          u_out = u;
4.1348 +          result = 1;
4.1349 +        }
4.1350 +      }
4.1351 +      // Search the cycle along the path form the second node to the root
4.1352 +      for (int u = second; u != join; u = _parent[u]) {
4.1353 +        e = _pred[u];
4.1354 +        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
4.1355 +        if (d <= delta) {
4.1356 +          delta = d;
4.1357 +          u_out = u;
4.1358 +          result = 2;
4.1359 +        }
4.1360 +      }
4.1361 +
4.1362 +      if (result == 1) {
4.1363 +        u_in = first;
4.1364 +        v_in = second;
4.1365 +      } else {
4.1366 +        u_in = second;
4.1367 +        v_in = first;
4.1368 +      }
4.1369 +      return result != 0;
4.1370 +    }
4.1371 +
4.1372 +    // Change _flow and _state vectors
4.1373 +    void changeFlow(bool change) {
4.1374 +      // Augment along the cycle
4.1375 +      if (delta > 0) {
4.1376 +        Flow val = _state[in_arc] * delta;
4.1377 +        _flow[in_arc] += val;
4.1378 +        for (int u = _source[in_arc]; u != join; u = _parent[u]) {
4.1379 +          _flow[_pred[u]] += _forward[u] ? -val : val;
4.1380 +        }
4.1381 +        for (int u = _target[in_arc]; u != join; u = _parent[u]) {
4.1382 +          _flow[_pred[u]] += _forward[u] ? val : -val;
4.1383 +        }
4.1384 +      }
4.1385 +      // Update the state of the entering and leaving arcs
4.1386 +      if (change) {
4.1387 +        _state[in_arc] = STATE_TREE;
4.1388 +        _state[_pred[u_out]] =
4.1389 +          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
4.1390 +      } else {
4.1391 +        _state[in_arc] = -_state[in_arc];
4.1392 +      }
4.1393 +    }
4.1394 +
4.1395 +    // Update the tree structure
4.1396 +    void updateTreeStructure() {
4.1397 +      int u, w;
4.1399 +      int old_succ_num = _succ_num[u_out];
4.1400 +      int old_last_succ = _last_succ[u_out];
4.1401 +      v_out = _parent[u_out];
4.1402 +
4.1403 +      u = _last_succ[u_in];  // the last successor of u_in
4.1404 +      right = _thread[u];    // the node after it
4.1405 +
4.1406 +      // Handle the case when old_rev_thread equals to v_in
4.1407 +      // (it also means that join and v_out coincide)
4.1408 +      if (old_rev_thread == v_in) {
4.1409 +        last = _thread[_last_succ[u_out]];
4.1410 +      } else {
4.1411 +        last = _thread[v_in];
4.1412 +      }
4.1413 +
4.1414 +      // Update _thread and _parent along the stem nodes (i.e. the nodes
4.1415 +      // between u_in and u_out, whose parent have to be changed)
4.1416 +      _thread[v_in] = stem = u_in;
4.1417 +      _dirty_revs.clear();
4.1418 +      _dirty_revs.push_back(v_in);
4.1419 +      par_stem = v_in;
4.1420 +      while (stem != u_out) {
4.1421 +        // Insert the next stem node into the thread list
4.1422 +        new_stem = _parent[stem];
4.1423 +        _thread[u] = new_stem;
4.1424 +        _dirty_revs.push_back(u);
4.1425 +
4.1426 +        // Remove the subtree of stem from the thread list
4.1427 +        w = _rev_thread[stem];
4.1428 +        _thread[w] = right;
4.1429 +        _rev_thread[right] = w;
4.1430 +
4.1431 +        // Change the parent node and shift stem nodes
4.1432 +        _parent[stem] = par_stem;
4.1433 +        par_stem = stem;
4.1434 +        stem = new_stem;
4.1435 +
4.1436 +        // Update u and right
4.1437 +        u = _last_succ[stem] == _last_succ[par_stem] ?
4.1438 +          _rev_thread[par_stem] : _last_succ[stem];
4.1439 +        right = _thread[u];
4.1440 +      }
4.1441 +      _parent[u_out] = par_stem;
4.1442 +      _thread[u] = last;
4.1443 +      _rev_thread[last] = u;
4.1444 +      _last_succ[u_out] = u;
4.1445 +
4.1446 +      // Remove the subtree of u_out from the thread list except for
4.1447 +      // the case when old_rev_thread equals to v_in
4.1448 +      // (it also means that join and v_out coincide)
4.1449 +      if (old_rev_thread != v_in) {
4.1452 +      }
4.1453 +
4.1454 +      // Update _rev_thread using the new _thread values
4.1455 +      for (int i = 0; i < int(_dirty_revs.size()); ++i) {
4.1456 +        u = _dirty_revs[i];
4.1458 +      }
4.1459 +
4.1460 +      // Update _pred, _forward, _last_succ and _succ_num for the
4.1461 +      // stem nodes from u_out to u_in
4.1462 +      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
4.1463 +      u = u_out;
4.1464 +      while (u != u_in) {
4.1465 +        w = _parent[u];
4.1466 +        _pred[u] = _pred[w];
4.1467 +        _forward[u] = !_forward[w];
4.1468 +        tmp_sc += _succ_num[u] - _succ_num[w];
4.1469 +        _succ_num[u] = tmp_sc;
4.1470 +        _last_succ[w] = tmp_ls;
4.1471 +        u = w;
4.1472 +      }
4.1473 +      _pred[u_in] = in_arc;
4.1474 +      _forward[u_in] = (u_in == _source[in_arc]);
4.1475 +      _succ_num[u_in] = old_succ_num;
4.1476 +
4.1477 +      // Set limits for updating _last_succ form v_in and v_out
4.1478 +      // towards the root
4.1479 +      int up_limit_in = -1;
4.1480 +      int up_limit_out = -1;
4.1481 +      if (_last_succ[join] == v_in) {
4.1482 +        up_limit_out = join;
4.1483 +      } else {
4.1484 +        up_limit_in = join;
4.1485 +      }
4.1486 +
4.1487 +      // Update _last_succ from v_in towards the root
4.1488 +      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
4.1489 +           u = _parent[u]) {
4.1490 +        _last_succ[u] = _last_succ[u_out];
4.1491 +      }
4.1492 +      // Update _last_succ from v_out towards the root
4.1493 +      if (join != old_rev_thread && v_in != old_rev_thread) {
4.1494 +        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
4.1495 +             u = _parent[u]) {
4.1496 +          _last_succ[u] = old_rev_thread;
4.1497 +        }
4.1498 +      } else {
4.1499 +        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
4.1500 +             u = _parent[u]) {
4.1501 +          _last_succ[u] = _last_succ[u_out];
4.1502 +        }
4.1503 +      }
4.1504 +
4.1505 +      // Update _succ_num from v_in to join
4.1506 +      for (u = v_in; u != join; u = _parent[u]) {
4.1507 +        _succ_num[u] += old_succ_num;
4.1508 +      }
4.1509 +      // Update _succ_num from v_out to join
4.1510 +      for (u = v_out; u != join; u = _parent[u]) {
4.1511 +        _succ_num[u] -= old_succ_num;
4.1512 +      }
4.1513 +    }
4.1514 +
4.1515 +    // Update potentials
4.1516 +    void updatePotential() {
4.1517 +      Cost sigma = _forward[u_in] ?
4.1518 +        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
4.1519 +        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
4.1520 +      // Update potentials in the subtree, which has been moved
4.1521 +      int end = _thread[_last_succ[u_in]];
4.1522 +      for (int u = u_in; u != end; u = _thread[u]) {
4.1523 +        _pi[u] += sigma;
4.1524 +      }
4.1525 +    }
4.1526 +
4.1527 +    // Execute the algorithm
4.1528 +    bool start(PivotRule pivot_rule) {
4.1529 +      // Select the pivot rule implementation
4.1530 +      switch (pivot_rule) {
4.1531 +        case FIRST_ELIGIBLE:
4.1532 +          return start<FirstEligiblePivotRule>();
4.1533 +        case BEST_ELIGIBLE:
4.1534 +          return start<BestEligiblePivotRule>();
4.1535 +        case BLOCK_SEARCH:
4.1536 +          return start<BlockSearchPivotRule>();
4.1537 +        case CANDIDATE_LIST:
4.1538 +          return start<CandidateListPivotRule>();
4.1539 +        case ALTERING_LIST:
4.1540 +          return start<AlteringListPivotRule>();
4.1541 +      }
4.1542 +      return false;
4.1543 +    }
4.1544 +
4.1545 +    template <typename PivotRuleImpl>
4.1546 +    bool start() {
4.1547 +      PivotRuleImpl pivot(*this);
4.1548 +
4.1549 +      // Execute the Network Simplex algorithm
4.1550 +      while (pivot.findEnteringArc()) {
4.1551 +        findJoinNode();
4.1552 +        bool change = findLeavingArc();
4.1553 +        changeFlow(change);
4.1554 +        if (change) {
4.1555 +          updateTreeStructure();
4.1556 +          updatePotential();
4.1557 +        }
4.1558 +      }
4.1559 +
4.1560 +      // Copy flow values to _flow_map
4.1561 +      if (_plower) {
4.1562 +        for (int i = 0; i != _arc_num; ++i) {
4.1563 +          Arc e = _arc_ref[i];
4.1564 +          _flow_map->set(e, (*_plower)[e] + _flow[i]);
4.1565 +        }
4.1566 +      } else {
4.1567 +        for (int i = 0; i != _arc_num; ++i) {
4.1568 +          _flow_map->set(_arc_ref[i], _flow[i]);
4.1569 +        }
4.1570 +      }
4.1571 +      // Copy potential values to _potential_map
4.1572 +      for (NodeIt n(_graph); n != INVALID; ++n) {
4.1573 +        _potential_map->set(n, _pi[_node_id[n]]);
4.1574 +      }
4.1575 +
4.1576 +      return true;
4.1577 +    }
4.1578 +
4.1579 +  }; //class NetworkSimplex
4.1580 +
4.1581 +  ///@}
4.1582 +
4.1583 +} //namespace lemon
4.1584 +
4.1585 +#endif //LEMON_NETWORK_SIMPLEX_H

     5.1 --- a/lemon/preflow.h	Tue Apr 21 13:08:19 2009 +0100
5.2 +++ b/lemon/preflow.h	Tue Apr 21 15:18:54 2009 +0100
5.3 @@ -46,18 +46,18 @@
5.4      typedef CAP CapacityMap;
5.5
5.6      /// \brief The type of the flow values.
5.7 -    typedef typename CapacityMap::Value Value;
5.8 +    typedef typename CapacityMap::Value Flow;
5.9
5.10      /// \brief The type of the map that stores the flow values.
5.11      ///
5.12      /// The type of the map that stores the flow values.
5.13      /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
5.14 -    typedef typename Digraph::template ArcMap<Value> FlowMap;
5.15 +    typedef typename Digraph::template ArcMap<Flow> FlowMap;
5.16
5.17      /// \brief Instantiates a FlowMap.
5.18      ///
5.19      /// This function instantiates a \ref FlowMap.
5.20 -    /// \param digraph The digraph, to which we would like to define
5.21 +    /// \param digraph The digraph for which we would like to define
5.22      /// the flow map.
5.23      static FlowMap* createFlowMap(const Digraph& digraph) {
5.24        return new FlowMap(digraph);
5.25 @@ -74,7 +74,7 @@
5.26      /// \brief Instantiates an Elevator.
5.27      ///
5.28      /// This function instantiates an \ref Elevator.
5.29 -    /// \param digraph The digraph, to which we would like to define
5.30 +    /// \param digraph The digraph for which we would like to define
5.31      /// the elevator.
5.32      /// \param max_level The maximum level of the elevator.
5.33      static Elevator* createElevator(const Digraph& digraph, int max_level) {
5.34 @@ -84,7 +84,7 @@
5.35      /// \brief The tolerance used by the algorithm
5.36      ///
5.37      /// The tolerance used by the algorithm to handle inexact computation.
5.38 -    typedef lemon::Tolerance<Value> Tolerance;
5.39 +    typedef lemon::Tolerance<Flow> Tolerance;
5.40
5.41    };
5.42
5.43 @@ -125,7 +125,7 @@
5.44      ///The type of the capacity map.
5.45      typedef typename Traits::CapacityMap CapacityMap;
5.46      ///The type of the flow values.
5.47 -    typedef typename Traits::Value Value;
5.48 +    typedef typename Traits::Flow Flow;
5.49
5.50      ///The type of the flow map.
5.51      typedef typename Traits::FlowMap FlowMap;
5.52 @@ -151,7 +151,7 @@
5.53      Elevator* _level;
5.54      bool _local_level;
5.55
5.56 -    typedef typename Digraph::template NodeMap<Value> ExcessMap;
5.57 +    typedef typename Digraph::template NodeMap<Flow> ExcessMap;
5.58      ExcessMap* _excess;
5.59
5.60      Tolerance _tolerance;
5.61 @@ -470,7 +470,7 @@
5.62        }
5.63
5.64        for (NodeIt n(_graph); n != INVALID; ++n) {
5.65 -        Value excess = 0;
5.66 +        Flow excess = 0;
5.67          for (InArcIt e(_graph, n); e != INVALID; ++e) {
5.68            excess += (*_flow)[e];
5.69          }
5.70 @@ -519,7 +519,7 @@
5.71        _level->initFinish();
5.72
5.73        for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
5.74 -        Value rem = (*_capacity)[e] - (*_flow)[e];
5.75 +        Flow rem = (*_capacity)[e] - (*_flow)[e];
5.76          if (_tolerance.positive(rem)) {
5.77            Node u = _graph.target(e);
5.78            if ((*_level)[u] == _level->maxLevel()) continue;
5.79 @@ -531,7 +531,7 @@
5.80          }
5.81        }
5.82        for (InArcIt e(_graph, _source); e != INVALID; ++e) {
5.83 -        Value rem = (*_flow)[e];
5.84 +        Flow rem = (*_flow)[e];
5.85          if (_tolerance.positive(rem)) {
5.86            Node v = _graph.source(e);
5.87            if ((*_level)[v] == _level->maxLevel()) continue;
5.88 @@ -564,11 +564,11 @@
5.89          int num = _node_num;
5.90
5.91          while (num > 0 && n != INVALID) {
5.92 -          Value excess = (*_excess)[n];
5.93 +          Flow excess = (*_excess)[n];
5.94            int new_level = _level->maxLevel();
5.95
5.96            for (OutArcIt e(_graph, n); e != INVALID; ++e) {
5.97 -            Value rem = (*_capacity)[e] - (*_flow)[e];
5.98 +            Flow rem = (*_capacity)[e] - (*_flow)[e];
5.99              if (!_tolerance.positive(rem)) continue;
5.100              Node v = _graph.target(e);
5.101              if ((*_level)[v] < level) {
5.102 @@ -591,7 +591,7 @@
5.103            }
5.104
5.105            for (InArcIt e(_graph, n); e != INVALID; ++e) {
5.106 -            Value rem = (*_flow)[e];
5.107 +            Flow rem = (*_flow)[e];
5.108              if (!_tolerance.positive(rem)) continue;
5.109              Node v = _graph.source(e);
5.110              if ((*_level)[v] < level) {
5.111 @@ -637,11 +637,11 @@
5.112
5.113          num = _node_num * 20;
5.114          while (num > 0 && n != INVALID) {
5.115 -          Value excess = (*_excess)[n];
5.116 +          Flow excess = (*_excess)[n];
5.117            int new_level = _level->maxLevel();
5.118
5.119            for (OutArcIt e(_graph, n); e != INVALID; ++e) {
5.120 -            Value rem = (*_capacity)[e] - (*_flow)[e];
5.121 +            Flow rem = (*_capacity)[e] - (*_flow)[e];
5.122              if (!_tolerance.positive(rem)) continue;
5.123              Node v = _graph.target(e);
5.124              if ((*_level)[v] < level) {
5.125 @@ -664,7 +664,7 @@
5.126            }
5.127
5.128            for (InArcIt e(_graph, n); e != INVALID; ++e) {
5.129 -            Value rem = (*_flow)[e];
5.130 +            Flow rem = (*_flow)[e];
5.131              if (!_tolerance.positive(rem)) continue;
5.132              Node v = _graph.source(e);
5.133              if ((*_level)[v] < level) {
5.134 @@ -778,12 +778,12 @@
5.135
5.136        Node n;
5.137        while ((n = _level->highestActive()) != INVALID) {
5.138 -        Value excess = (*_excess)[n];
5.139 +        Flow excess = (*_excess)[n];
5.140          int level = _level->highestActiveLevel();
5.141          int new_level = _level->maxLevel();
5.142
5.143          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
5.144 -          Value rem = (*_capacity)[e] - (*_flow)[e];
5.145 +          Flow rem = (*_capacity)[e] - (*_flow)[e];
5.146            if (!_tolerance.positive(rem)) continue;
5.147            Node v = _graph.target(e);
5.148            if ((*_level)[v] < level) {
5.149 @@ -806,7 +806,7 @@
5.150          }
5.151
5.152          for (InArcIt e(_graph, n); e != INVALID; ++e) {
5.153 -          Value rem = (*_flow)[e];
5.154 +          Flow rem = (*_flow)[e];
5.155            if (!_tolerance.positive(rem)) continue;
5.156            Node v = _graph.source(e);
5.157            if ((*_level)[v] < level) {
5.158 @@ -897,7 +897,7 @@
5.159      ///
5.160      /// \pre Either \ref run() or \ref init() must be called before
5.161      /// using this function.
5.162 -    Value flowValue() const {
5.163 +    Flow flowValue() const {
5.164        return (*_excess)[_target];
5.165      }
5.166
5.167 @@ -908,7 +908,7 @@
5.168      ///
5.169      /// \pre Either \ref run() or \ref init() must be called before
5.170      /// using this function.
5.171 -    Value flow(const Arc& arc) const {
5.172 +    Flow flow(const Arc& arc) const {
5.173        return (*_flow)[arc];
5.174      }
5.175

     6.1 --- a/test/CMakeLists.txt	Tue Apr 21 13:08:19 2009 +0100
6.2 +++ b/test/CMakeLists.txt	Tue Apr 21 15:18:54 2009 +0100
6.3 @@ -31,6 +31,7 @@
6.4    maps_test
6.5    matching_test
6.6    min_cost_arborescence_test
6.7 +  min_cost_flow_test
6.8    path_test
6.9    preflow_test

     7.1 --- a/test/Makefile.am	Tue Apr 21 13:08:19 2009 +0100
7.2 +++ b/test/Makefile.am	Tue Apr 21 15:18:54 2009 +0100
7.3 @@ -27,6 +27,7 @@
7.4  	test/maps_test \
7.5  	test/matching_test \
7.6  	test/min_cost_arborescence_test \
7.7 +	test/min_cost_flow_test \
7.8  	test/path_test \
7.9  	test/preflow_test \
7.11 @@ -72,6 +73,7 @@
7.12  test_mip_test_SOURCES = test/mip_test.cc
7.13  test_matching_test_SOURCES = test/matching_test.cc
7.14  test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
7.15 +test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
7.16  test_path_test_SOURCES = test/path_test.cc
7.17  test_preflow_test_SOURCES = test/preflow_test.cc

     8.1 --- a/test/circulation_test.cc	Tue Apr 21 13:08:19 2009 +0100
8.2 +++ b/test/circulation_test.cc	Tue Apr 21 15:18:54 2009 +0100
8.3 @@ -57,7 +57,7 @@
8.4    typedef Digraph::Node Node;
8.5    typedef Digraph::Arc Arc;
8.6    typedef concepts::ReadMap<Arc,VType> CapMap;
8.7 -  typedef concepts::ReadMap<Node,VType> DeltaMap;
8.8 +  typedef concepts::ReadMap<Node,VType> SupplyMap;
8.9    typedef concepts::ReadWriteMap<Arc,VType> FlowMap;
8.10    typedef concepts::WriteMap<Node,bool> BarrierMap;
8.11
8.12 @@ -68,24 +68,24 @@
8.13    Node n;
8.14    Arc a;
8.15    CapMap lcap, ucap;
8.16 -  DeltaMap delta;
8.17 +  SupplyMap supply;
8.18    FlowMap flow;
8.19    BarrierMap bar;
8.20    VType v;
8.21    bool b;
8.22
8.23 -  typedef Circulation<Digraph, CapMap, CapMap, DeltaMap>
8.24 +  typedef Circulation<Digraph, CapMap, CapMap, SupplyMap>
8.25              ::SetFlowMap<FlowMap>
8.26              ::SetElevator<Elev>
8.28              ::Create CirculationType;
8.29 -  CirculationType circ_test(g, lcap, ucap, delta);
8.30 +  CirculationType circ_test(g, lcap, ucap, supply);
8.31    const CirculationType& const_circ_test = circ_test;
8.32
8.33    circ_test
8.34 -    .lowerCapMap(lcap)
8.35 -    .upperCapMap(ucap)
8.36 -    .deltaMap(delta)
8.37 +    .lowerMap(lcap)
8.38 +    .upperMap(ucap)
8.39 +    .supplyMap(supply)
8.40      .flowMap(flow);
8.41
8.42    circ_test.init();

     9.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
9.2 +++ b/test/min_cost_flow_test.cc	Tue Apr 21 15:18:54 2009 +0100
9.3 @@ -0,0 +1,339 @@
9.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
9.5 + *
9.6 + * This file is a part of LEMON, a generic C++ optimization library.
9.7 + *
9.8 + * Copyright (C) 2003-2009
9.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
9.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
9.11 + *
9.12 + * Permission to use, modify and distribute this software is granted
9.13 + * provided that this copyright notice appears in all copies. For
9.14 + * precise terms see the accompanying LICENSE file.
9.15 + *
9.16 + * This software is provided "AS IS" with no warranty of any kind,
9.17 + * express or implied, and with no claim as to its suitability for any
9.18 + * purpose.
9.19 + *
9.20 + */
9.21 +
9.22 +#include <iostream>
9.23 +#include <fstream>
9.24 +
9.25 +#include <lemon/list_graph.h>
9.27 +
9.28 +#include <lemon/network_simplex.h>
9.29 +
9.30 +#include <lemon/concepts/digraph.h>
9.31 +#include <lemon/concept_check.h>
9.32 +
9.33 +#include "test_tools.h"
9.34 +
9.35 +using namespace lemon;
9.36 +
9.37 +char test_lgf[] =
9.38 +  "@nodes\n"
9.39 +  "label  sup1 sup2 sup3 sup4 sup5\n"
9.40 +  "    1    20   27    0   20   30\n"
9.41 +  "    2    -4    0    0   -8   -3\n"
9.42 +  "    3     0    0    0    0    0\n"
9.43 +  "    4     0    0    0    0    0\n"
9.44 +  "    5     9    0    0    6   11\n"
9.45 +  "    6    -6    0    0   -5   -6\n"
9.46 +  "    7     0    0    0    0    0\n"
9.47 +  "    8     0    0    0    0    3\n"
9.48 +  "    9     3    0    0    0    0\n"
9.49 +  "   10    -2    0    0   -7   -2\n"
9.50 +  "   11     0    0    0  -10    0\n"
9.51 +  "   12   -20  -27    0  -30  -20\n"
9.52 +  "\n"
9.53 +  "@arcs\n"
9.54 +  "       cost  cap low1 low2\n"
9.55 +  " 1  2    70   11    0    8\n"
9.56 +  " 1  3   150    3    0    1\n"
9.57 +  " 1  4    80   15    0    2\n"
9.58 +  " 2  8    80   12    0    0\n"
9.59 +  " 3  5   140    5    0    3\n"
9.60 +  " 4  6    60   10    0    1\n"
9.61 +  " 4  7    80    2    0    0\n"
9.62 +  " 4  8   110    3    0    0\n"
9.63 +  " 5  7    60   14    0    0\n"
9.64 +  " 5 11   120   12    0    0\n"
9.65 +  " 6  3     0    3    0    0\n"
9.66 +  " 6  9   140    4    0    0\n"
9.67 +  " 6 10    90    8    0    0\n"
9.68 +  " 7  1    30    5    0    0\n"
9.69 +  " 8 12    60   16    0    4\n"
9.70 +  " 9 12    50    6    0    0\n"
9.71 +  "10 12    70   13    0    5\n"
9.72 +  "10  2   100    7    0    0\n"
9.73 +  "10  7    60   10    0    0\n"
9.74 +  "11 10    20   14    0    6\n"
9.75 +  "12 11    30   10    0    0\n"
9.76 +  "\n"
9.77 +  "@attributes\n"
9.78 +  "source 1\n"
9.79 +  "target 12\n";
9.80 +
9.81 +
9.82 +enum ProblemType {
9.83 +  EQ,
9.84 +  GEQ,
9.85 +  LEQ
9.86 +};
9.87 +
9.88 +// Check the interface of an MCF algorithm
9.89 +template <typename GR, typename Flow, typename Cost>
9.90 +class McfClassConcept
9.91 +{
9.92 +public:
9.93 +
9.94 +  template <typename MCF>
9.95 +  struct Constraints {
9.96 +    void constraints() {
9.97 +      checkConcept<concepts::Digraph, GR>();
9.98 +
9.99 +      MCF mcf(g);
9.100 +
9.101 +      b = mcf.reset()
9.102 +             .lowerMap(lower)
9.103 +             .upperMap(upper)
9.104 +             .capacityMap(upper)
9.105 +             .boundMaps(lower, upper)
9.106 +             .costMap(cost)
9.107 +             .supplyMap(sup)
9.108 +             .stSupply(n, n, k)
9.109 +             .flowMap(flow)
9.110 +             .potentialMap(pot)
9.111 +             .run();
9.112 +
9.113 +      const MCF& const_mcf = mcf;
9.114 +
9.115 +      const typename MCF::FlowMap &fm = const_mcf.flowMap();
9.116 +      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
9.117 +
9.118 +      v = const_mcf.totalCost();
9.119 +      double x = const_mcf.template totalCost<double>();
9.120 +      v = const_mcf.flow(a);
9.121 +      v = const_mcf.potential(n);
9.122 +
9.123 +      ignore_unused_variable_warning(fm);
9.124 +      ignore_unused_variable_warning(pm);
9.125 +      ignore_unused_variable_warning(x);
9.126 +    }
9.127 +
9.128 +    typedef typename GR::Node Node;
9.129 +    typedef typename GR::Arc Arc;
9.130 +    typedef concepts::ReadMap<Node, Flow> NM;
9.131 +    typedef concepts::ReadMap<Arc, Flow> FAM;
9.132 +    typedef concepts::ReadMap<Arc, Cost> CAM;
9.133 +
9.134 +    const GR &g;
9.135 +    const FAM &lower;
9.136 +    const FAM &upper;
9.137 +    const CAM &cost;
9.138 +    const NM &sup;
9.139 +    const Node &n;
9.140 +    const Arc &a;
9.141 +    const Flow &k;
9.142 +    Flow v;
9.143 +    bool b;
9.144 +
9.145 +    typename MCF::FlowMap &flow;
9.146 +    typename MCF::PotentialMap &pot;
9.147 +  };
9.148 +
9.149 +};
9.150 +
9.151 +
9.152 +// Check the feasibility of the given flow (primal soluiton)
9.153 +template < typename GR, typename LM, typename UM,
9.154 +           typename SM, typename FM >
9.155 +bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
9.156 +                const SM& supply, const FM& flow,
9.157 +                ProblemType type = EQ )
9.158 +{
9.159 +  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
9.160 +
9.161 +  for (ArcIt e(gr); e != INVALID; ++e) {
9.162 +    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
9.163 +  }
9.164 +
9.165 +  for (NodeIt n(gr); n != INVALID; ++n) {
9.166 +    typename SM::Value sum = 0;
9.167 +    for (OutArcIt e(gr, n); e != INVALID; ++e)
9.168 +      sum += flow[e];
9.169 +    for (InArcIt e(gr, n); e != INVALID; ++e)
9.170 +      sum -= flow[e];
9.171 +    bool b = (type ==  EQ && sum == supply[n]) ||
9.172 +             (type == GEQ && sum >= supply[n]) ||
9.173 +             (type == LEQ && sum <= supply[n]);
9.174 +    if (!b) return false;
9.175 +  }
9.176 +
9.177 +  return true;
9.178 +}
9.179 +
9.180 +// Check the feasibility of the given potentials (dual soluiton)
9.181 +// using the "Complementary Slackness" optimality condition
9.182 +template < typename GR, typename LM, typename UM,
9.183 +           typename CM, typename SM, typename FM, typename PM >
9.184 +bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
9.185 +                     const CM& cost, const SM& supply, const FM& flow,
9.186 +                     const PM& pi )
9.187 +{
9.188 +  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
9.189 +
9.190 +  bool opt = true;
9.191 +  for (ArcIt e(gr); opt && e != INVALID; ++e) {
9.192 +    typename CM::Value red_cost =
9.193 +      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
9.194 +    opt = red_cost == 0 ||
9.195 +          (red_cost > 0 && flow[e] == lower[e]) ||
9.196 +          (red_cost < 0 && flow[e] == upper[e]);
9.197 +  }
9.198 +
9.199 +  for (NodeIt n(gr); opt && n != INVALID; ++n) {
9.200 +    typename SM::Value sum = 0;
9.201 +    for (OutArcIt e(gr, n); e != INVALID; ++e)
9.202 +      sum += flow[e];
9.203 +    for (InArcIt e(gr, n); e != INVALID; ++e)
9.204 +      sum -= flow[e];
9.205 +    opt = (sum == supply[n]) || (pi[n] == 0);
9.206 +  }
9.207 +
9.208 +  return opt;
9.209 +}
9.210 +
9.211 +// Run a minimum cost flow algorithm and check the results
9.212 +template < typename MCF, typename GR,
9.213 +           typename LM, typename UM,
9.214 +           typename CM, typename SM >
9.215 +void checkMcf( const MCF& mcf, bool mcf_result,
9.216 +               const GR& gr, const LM& lower, const UM& upper,
9.217 +               const CM& cost, const SM& supply,
9.218 +               bool result, typename CM::Value total,
9.219 +               const std::string &test_id = "",
9.220 +               ProblemType type = EQ )
9.221 +{
9.222 +  check(mcf_result == result, "Wrong result " + test_id);
9.223 +  if (result) {
9.224 +    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
9.225 +          "The flow is not feasible " + test_id);
9.226 +    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
9.227 +    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
9.228 +                         mcf.potentialMap()),
9.229 +          "Wrong potentials " + test_id);
9.230 +  }
9.231 +}
9.232 +
9.233 +int main()
9.234 +{
9.235 +  // Check the interfaces
9.236 +  {
9.237 +    typedef int Flow;
9.238 +    typedef int Cost;
9.239 +    // TODO: This typedef should be enabled if the standard maps are
9.240 +    // reference maps in the graph concepts (See #190).
9.241 +/**/
9.242 +    //typedef concepts::Digraph GR;
9.243 +    typedef ListDigraph GR;
9.244 +/**/
9.245 +    checkConcept< McfClassConcept<GR, Flow, Cost>,
9.246 +                  NetworkSimplex<GR, Flow, Cost> >();
9.247 +  }
9.248 +
9.249 +  // Run various MCF tests
9.250 +  typedef ListDigraph Digraph;
9.251 +  DIGRAPH_TYPEDEFS(ListDigraph);
9.252 +
9.253 +  // Read the test digraph
9.254 +  Digraph gr;
9.255 +  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
9.256 +  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
9.257 +  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
9.258 +  Node v, w;
9.259 +
9.260 +  std::istringstream input(test_lgf);
9.261 +  DigraphReader<Digraph>(gr, input)
9.262 +    .arcMap("cost", c)
9.263 +    .arcMap("cap", u)
9.264 +    .arcMap("low1", l1)
9.265 +    .arcMap("low2", l2)
9.266 +    .nodeMap("sup1", s1)
9.267 +    .nodeMap("sup2", s2)
9.268 +    .nodeMap("sup3", s3)
9.269 +    .nodeMap("sup4", s4)
9.270 +    .nodeMap("sup5", s5)
9.271 +    .node("source", v)
9.272 +    .node("target", w)
9.273 +    .run();
9.274 +
9.275 +  // A. Test NetworkSimplex with the default pivot rule
9.276 +  {
9.277 +    NetworkSimplex<Digraph> mcf(gr);
9.278 +
9.279 +    // Check the equality form
9.280 +    mcf.upperMap(u).costMap(c);
9.281 +    checkMcf(mcf, mcf.supplyMap(s1).run(),
9.282 +             gr, l1, u, c, s1, true,  5240, "#A1");
9.283 +    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
9.284 +             gr, l1, u, c, s2, true,  7620, "#A2");
9.285 +    mcf.lowerMap(l2);
9.286 +    checkMcf(mcf, mcf.supplyMap(s1).run(),
9.287 +             gr, l2, u, c, s1, true,  5970, "#A3");
9.288 +    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
9.289 +             gr, l2, u, c, s2, true,  8010, "#A4");
9.290 +    mcf.reset();
9.291 +    checkMcf(mcf, mcf.supplyMap(s1).run(),
9.292 +             gr, l1, cu, cc, s1, true,  74, "#A5");
9.293 +    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
9.294 +             gr, l2, cu, cc, s2, true,  94, "#A6");
9.295 +    mcf.reset();
9.296 +    checkMcf(mcf, mcf.run(),
9.297 +             gr, l1, cu, cc, s3, true,   0, "#A7");
9.298 +    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
9.299 +             gr, l2, u, cc, s3, false,   0, "#A8");
9.300 +
9.301 +    // Check the GEQ form
9.302 +    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
9.303 +    checkMcf(mcf, mcf.run(),
9.304 +             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
9.305 +    mcf.problemType(mcf.GEQ);
9.306 +    checkMcf(mcf, mcf.lowerMap(l2).run(),
9.307 +             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
9.308 +    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
9.309 +    checkMcf(mcf, mcf.run(),
9.310 +             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
9.311 +
9.312 +    // Check the LEQ form
9.313 +    mcf.reset().problemType(mcf.LEQ);
9.314 +    mcf.upperMap(u).costMap(c).supplyMap(s5);
9.315 +    checkMcf(mcf, mcf.run(),
9.316 +             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
9.317 +    checkMcf(mcf, mcf.lowerMap(l2).run(),
9.318 +             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
9.319 +    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
9.320 +    checkMcf(mcf, mcf.run(),
9.321 +             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
9.322 +  }
9.323 +
9.324 +  // B. Test NetworkSimplex with each pivot rule
9.325 +  {
9.326 +    NetworkSimplex<Digraph> mcf(gr);
9.327 +    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
9.328 +
9.329 +    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
9.330 +             gr, l2, u, c, s1, true,  5970, "#B1");
9.331 +    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
9.332 +             gr, l2, u, c, s1, true,  5970, "#B2");
9.333 +    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
9.334 +             gr, l2, u, c, s1, true,  5970, "#B3");
9.335 +    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
9.336 +             gr, l2, u, c, s1, true,  5970, "#B4");
9.337 +    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
9.338 +             gr, l2, u, c, s1, true,  5970, "#B5");
9.339 +  }
9.340 +
9.341 +  return 0;
9.342 +}

    10.1 --- a/tools/dimacs-solver.cc	Tue Apr 21 13:08:19 2009 +0100
10.2 +++ b/tools/dimacs-solver.cc	Tue Apr 21 15:18:54 2009 +0100
10.3 @@ -43,6 +43,7 @@
10.4  #include <lemon/dijkstra.h>
10.5  #include <lemon/preflow.h>
10.6  #include <lemon/matching.h>
10.7 +#include <lemon/network_simplex.h>
10.8
10.9  using namespace lemon;
10.10  typedef SmartDigraph Digraph;
10.11 @@ -90,6 +91,28 @@
10.12    if(report) std::cerr << "\nMax flow value: " << pre.flowValue() << '\n';
10.13  }
10.14
10.15 +template<class Value>
10.16 +void solve_min(ArgParser &ap, std::istream &is, std::ostream &,
10.17 +               DimacsDescriptor &desc)
10.18 +{
10.19 +  bool report = !ap.given("q");
10.20 +  Digraph g;
10.21 +  Digraph::ArcMap<Value> lower(g), cap(g), cost(g);
10.22 +  Digraph::NodeMap<Value> sup(g);
10.23 +  Timer ti;
10.24 +  ti.restart();
10.25 +  readDimacsMin(is, g, lower, cap, cost, sup, 0, desc);
10.26 +  if (report) std::cerr << "Read the file: " << ti << '\n';
10.27 +  ti.restart();
10.28 +  NetworkSimplex<Digraph, Value> ns(g);
10.29 +  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
10.30 +  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
10.31 +  ti.restart();
10.32 +  ns.run();
10.33 +  if (report) std::cerr << "Run NetworkSimplex: " << ti << '\n';
10.34 +  if (report) std::cerr << "\nMin flow cost: " << ns.totalCost() << '\n';
10.35 +}
10.36 +
10.37  void solve_mat(ArgParser &ap, std::istream &is, std::ostream &,
10.38                DimacsDescriptor &desc)
10.39  {
10.40 @@ -128,8 +151,7 @@
10.41    switch(desc.type)
10.42      {
10.43      case DimacsDescriptor::MIN:
10.44 -      std::cerr <<
10.45 -        "\n\n Sorry, the min. cost flow solver is not yet available.\n";
10.46 +      solve_min<Value>(ap,is,os,desc);
10.47        break;
10.48      case DimacsDescriptor::MAX:
10.49        solve_max<Value>(ap,is,os,infty,desc);