Merge and fix
authorAlpar Juttner <alpar@cs.elte.hu>
Tue, 21 Apr 2009 15:18:54 +0100
changeset 60385cb3aa71cce
parent 592 0ba8dfce7259
parent 602 dacc2cee2b4c
child 604 0c8e5c688440
child 605 b1811c363299
child 611 ec817dfc2cb7
Merge and fix
doc/groups.dox
lemon/Makefile.am
lemon/circulation.h
lemon/preflow.h
test/CMakeLists.txt
test/Makefile.am
test/circulation_test.cc
tools/dimacs-solver.cc
     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.6  	lemon/nauty_reader.h \
     2.7 +	lemon/network_simplex.h \
     2.8  	lemon/path.h \
     2.9  	lemon/preflow.h \
    2.10  	lemon/radix_sort.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.38 +#include <lemon/adaptors.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.1015 +      return totalCost<Cost>();
  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.1096 +      _thread.resize(all_node_num);
  4.1097 +      _rev_thread.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.1398 +      int old_rev_thread = _rev_thread[u_out];
  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.1450 +        _thread[old_rev_thread] = right;
  4.1451 +        _rev_thread[right] = old_rev_thread;
  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.1457 +        _rev_thread[_thread[u]] = u;
  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
    6.10    radix_sort_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.10  	test/radix_sort_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
    7.18  test_radix_sort_test_SOURCES = test/radix_sort_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.27              ::SetStandardElevator<LinkedElev>
    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.26 +#include <lemon/lgf_reader.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);