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 = ↦
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 = ↦
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 = ↦
3.345 + Circulation& supplyMap(const SupplyMap& map) {
3.346 + _supply = ↦
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 = ↦
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 = ↦
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 ⊃
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);