1.1 --- a/doc/groups.dox Sun Apr 26 16:36:23 2009 +0100
1.2 +++ b/doc/groups.dox Wed Apr 29 03:15:24 2009 +0200
1.3 @@ -352,17 +352,17 @@
1.4 minimum total cost from a set of supply nodes to a set of demand nodes
1.5 in a network with capacity constraints (lower and upper bounds)
1.6 and arc costs.
1.7 -Formally, let \f$G=(V,A)\f$ be a digraph,
1.8 -\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
1.9 +Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
1.10 +\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
1.11 upper bounds for the flow values on the arcs, for which
1.12 -\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
1.13 -\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
1.14 -on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
1.15 +\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
1.16 +\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
1.17 +on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
1.18 signed supply values of the nodes.
1.19 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
1.20 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
1.21 \f$-sup(u)\f$ demand.
1.22 -A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
1.23 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
1.24 of the following optimization problem.
1.25
1.26 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.27 @@ -404,7 +404,7 @@
1.28
1.29 The dual solution of the minimum cost flow problem is represented by node
1.30 potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
1.31 -An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
1.32 +An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
1.33 is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
1.34 node potentials the following \e complementary \e slackness optimality
1.35 conditions hold.
1.36 @@ -413,15 +413,15 @@
1.37 - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
1.38 - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
1.39 - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
1.40 - - For all \f$u\in V\f$:
1.41 + - For all \f$u\in V\f$ nodes:
1.42 - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
1.43 then \f$\pi(u)=0\f$.
1.44
1.45 Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
1.46 -\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
1.47 +\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
1.48 \f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
1.49
1.50 -All algorithms provide dual solution (node potentials) as well
1.51 +All algorithms provide dual solution (node potentials) as well,
1.52 if an optimal flow is found.
1.53
1.54 LEMON contains several algorithms for solving minimum cost flow problems.
2.1 --- a/lemon/network_simplex.h Sun Apr 26 16:36:23 2009 +0100
2.2 +++ b/lemon/network_simplex.h Wed Apr 29 03:15:24 2009 +0200
2.3 @@ -30,9 +30,6 @@
2.4
2.5 #include <lemon/core.h>
2.6 #include <lemon/math.h>
2.7 -#include <lemon/maps.h>
2.8 -#include <lemon/circulation.h>
2.9 -#include <lemon/adaptors.h>
2.10
2.11 namespace lemon {
2.12
2.13 @@ -50,8 +47,13 @@
2.14 ///
2.15 /// In general this class is the fastest implementation available
2.16 /// in LEMON for the minimum cost flow problem.
2.17 - /// Moreover it supports both direction of the supply/demand inequality
2.18 - /// constraints. For more information see \ref ProblemType.
2.19 + /// Moreover it supports both directions of the supply/demand inequality
2.20 + /// constraints. For more information see \ref SupplyType.
2.21 + ///
2.22 + /// Most of the parameters of the problem (except for the digraph)
2.23 + /// can be given using separate functions, and the algorithm can be
2.24 + /// executed using the \ref run() function. If some parameters are not
2.25 + /// specified, then default values will be used.
2.26 ///
2.27 /// \tparam GR The digraph type the algorithm runs on.
2.28 /// \tparam F The value type used for flow amounts, capacity bounds
2.29 @@ -88,11 +90,80 @@
2.30
2.31 public:
2.32
2.33 - /// \brief Enum type for selecting the pivot rule.
2.34 + /// \brief Problem type constants for the \c run() function.
2.35 ///
2.36 - /// Enum type for selecting the pivot rule for the \ref run()
2.37 + /// Enum type containing the problem type constants that can be
2.38 + /// returned by the \ref run() function of the algorithm.
2.39 + enum ProblemType {
2.40 + /// The problem has no feasible solution (flow).
2.41 + INFEASIBLE,
2.42 + /// The problem has optimal solution (i.e. it is feasible and
2.43 + /// bounded), and the algorithm has found optimal flow and node
2.44 + /// potentials (primal and dual solutions).
2.45 + OPTIMAL,
2.46 + /// The objective function of the problem is unbounded, i.e.
2.47 + /// there is a directed cycle having negative total cost and
2.48 + /// infinite upper bound.
2.49 + UNBOUNDED
2.50 + };
2.51 +
2.52 + /// \brief Constants for selecting the type of the supply constraints.
2.53 + ///
2.54 + /// Enum type containing constants for selecting the supply type,
2.55 + /// i.e. the direction of the inequalities in the supply/demand
2.56 + /// constraints of the \ref min_cost_flow "minimum cost flow problem".
2.57 + ///
2.58 + /// The default supply type is \c GEQ, since this form is supported
2.59 + /// by other minimum cost flow algorithms and the \ref Circulation
2.60 + /// algorithm, as well.
2.61 + /// The \c LEQ problem type can be selected using the \ref supplyType()
2.62 /// function.
2.63 ///
2.64 + /// Note that the equality form is a special case of both supply types.
2.65 + enum SupplyType {
2.66 +
2.67 + /// This option means that there are <em>"greater or equal"</em>
2.68 + /// supply/demand constraints in the definition, i.e. the exact
2.69 + /// formulation of the problem is the following.
2.70 + /**
2.71 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
2.72 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
2.73 + sup(u) \quad \forall u\in V \f]
2.74 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
2.75 + */
2.76 + /// It means that the total demand must be greater or equal to the
2.77 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
2.78 + /// negative) and all the supplies have to be carried out from
2.79 + /// the supply nodes, but there could be demands that are not
2.80 + /// satisfied.
2.81 + GEQ,
2.82 + /// It is just an alias for the \c GEQ option.
2.83 + CARRY_SUPPLIES = GEQ,
2.84 +
2.85 + /// This option means that there are <em>"less or equal"</em>
2.86 + /// supply/demand constraints in the definition, i.e. the exact
2.87 + /// formulation of the problem is the following.
2.88 + /**
2.89 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
2.90 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
2.91 + sup(u) \quad \forall u\in V \f]
2.92 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
2.93 + */
2.94 + /// It means that the total demand must be less or equal to the
2.95 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
2.96 + /// positive) and all the demands have to be satisfied, but there
2.97 + /// could be supplies that are not carried out from the supply
2.98 + /// nodes.
2.99 + LEQ,
2.100 + /// It is just an alias for the \c LEQ option.
2.101 + SATISFY_DEMANDS = LEQ
2.102 + };
2.103 +
2.104 + /// \brief Constants for selecting the pivot rule.
2.105 + ///
2.106 + /// Enum type containing constants for selecting the pivot rule for
2.107 + /// the \ref run() function.
2.108 + ///
2.109 /// \ref NetworkSimplex provides five different pivot rule
2.110 /// implementations that significantly affect the running time
2.111 /// of the algorithm.
2.112 @@ -131,58 +202,6 @@
2.113 ALTERING_LIST
2.114 };
2.115
2.116 - /// \brief Enum type for selecting the problem type.
2.117 - ///
2.118 - /// Enum type for selecting the problem type, i.e. the direction of
2.119 - /// the inequalities in the supply/demand constraints of the
2.120 - /// \ref min_cost_flow "minimum cost flow problem".
2.121 - ///
2.122 - /// The default problem type is \c GEQ, since this form is supported
2.123 - /// by other minimum cost flow algorithms and the \ref Circulation
2.124 - /// algorithm as well.
2.125 - /// The \c LEQ problem type can be selected using the \ref problemType()
2.126 - /// function.
2.127 - ///
2.128 - /// Note that the equality form is a special case of both problem type.
2.129 - enum ProblemType {
2.130 -
2.131 - /// This option means that there are "<em>greater or equal</em>"
2.132 - /// constraints in the defintion, i.e. the exact formulation of the
2.133 - /// problem is the following.
2.134 - /**
2.135 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
2.136 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
2.137 - sup(u) \quad \forall u\in V \f]
2.138 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
2.139 - */
2.140 - /// It means that the total demand must be greater or equal to the
2.141 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
2.142 - /// negative) and all the supplies have to be carried out from
2.143 - /// the supply nodes, but there could be demands that are not
2.144 - /// satisfied.
2.145 - GEQ,
2.146 - /// It is just an alias for the \c GEQ option.
2.147 - CARRY_SUPPLIES = GEQ,
2.148 -
2.149 - /// This option means that there are "<em>less or equal</em>"
2.150 - /// constraints in the defintion, i.e. the exact formulation of the
2.151 - /// problem is the following.
2.152 - /**
2.153 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
2.154 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
2.155 - sup(u) \quad \forall u\in V \f]
2.156 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
2.157 - */
2.158 - /// It means that the total demand must be less or equal to the
2.159 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
2.160 - /// positive) and all the demands have to be satisfied, but there
2.161 - /// could be supplies that are not carried out from the supply
2.162 - /// nodes.
2.163 - LEQ,
2.164 - /// It is just an alias for the \c LEQ option.
2.165 - SATISFY_DEMANDS = LEQ
2.166 - };
2.167 -
2.168 private:
2.169
2.170 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
2.171 @@ -220,7 +239,9 @@
2.172 bool _pstsup;
2.173 Node _psource, _ptarget;
2.174 Flow _pstflow;
2.175 - ProblemType _ptype;
2.176 + SupplyType _stype;
2.177 +
2.178 + Flow _sum_supply;
2.179
2.180 // Result maps
2.181 FlowMap *_flow_map;
2.182 @@ -259,6 +280,15 @@
2.183 int stem, par_stem, new_stem;
2.184 Flow delta;
2.185
2.186 + public:
2.187 +
2.188 + /// \brief Constant for infinite upper bounds (capacities).
2.189 + ///
2.190 + /// Constant for infinite upper bounds (capacities).
2.191 + /// It is \c std::numeric_limits<Flow>::infinity() if available,
2.192 + /// \c std::numeric_limits<Flow>::max() otherwise.
2.193 + const Flow INF;
2.194 +
2.195 private:
2.196
2.197 // Implementation of the First Eligible pivot rule
2.198 @@ -661,17 +691,19 @@
2.199 NetworkSimplex(const GR& graph) :
2.200 _graph(graph),
2.201 _plower(NULL), _pupper(NULL), _pcost(NULL),
2.202 - _psupply(NULL), _pstsup(false), _ptype(GEQ),
2.203 + _psupply(NULL), _pstsup(false), _stype(GEQ),
2.204 _flow_map(NULL), _potential_map(NULL),
2.205 _local_flow(false), _local_potential(false),
2.206 - _node_id(graph)
2.207 + _node_id(graph),
2.208 + INF(std::numeric_limits<Flow>::has_infinity ?
2.209 + std::numeric_limits<Flow>::infinity() :
2.210 + std::numeric_limits<Flow>::max())
2.211 {
2.212 - LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
2.213 - std::numeric_limits<Flow>::is_signed,
2.214 - "The flow type of NetworkSimplex must be signed integer");
2.215 - LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
2.216 - std::numeric_limits<Cost>::is_signed,
2.217 - "The cost type of NetworkSimplex must be signed integer");
2.218 + // Check the value types
2.219 + LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,
2.220 + "The flow type of NetworkSimplex must be signed");
2.221 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
2.222 + "The cost type of NetworkSimplex must be signed");
2.223 }
2.224
2.225 /// Destructor.
2.226 @@ -689,17 +721,16 @@
2.227 /// \brief Set the lower bounds on the arcs.
2.228 ///
2.229 /// This function sets the lower bounds on the arcs.
2.230 - /// If neither this function nor \ref boundMaps() is used before
2.231 - /// calling \ref run(), the lower bounds will be set to zero
2.232 - /// on all arcs.
2.233 + /// If it is not used before calling \ref run(), the lower bounds
2.234 + /// will be set to zero on all arcs.
2.235 ///
2.236 /// \param map An arc map storing the lower bounds.
2.237 /// Its \c Value type must be convertible to the \c Flow type
2.238 /// of the algorithm.
2.239 ///
2.240 /// \return <tt>(*this)</tt>
2.241 - template <typename LOWER>
2.242 - NetworkSimplex& lowerMap(const LOWER& map) {
2.243 + template <typename LowerMap>
2.244 + NetworkSimplex& lowerMap(const LowerMap& map) {
2.245 delete _plower;
2.246 _plower = new FlowArcMap(_graph);
2.247 for (ArcIt a(_graph); a != INVALID; ++a) {
2.248 @@ -711,18 +742,17 @@
2.249 /// \brief Set the upper bounds (capacities) on the arcs.
2.250 ///
2.251 /// This function sets the upper bounds (capacities) on the arcs.
2.252 - /// If none of the functions \ref upperMap(), \ref capacityMap()
2.253 - /// and \ref boundMaps() is used before calling \ref run(),
2.254 - /// the upper bounds (capacities) will be set to
2.255 - /// \c std::numeric_limits<Flow>::max() on all arcs.
2.256 + /// If it is not used before calling \ref run(), the upper bounds
2.257 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
2.258 + /// unbounded from above on each arc).
2.259 ///
2.260 /// \param map An arc map storing the upper bounds.
2.261 /// Its \c Value type must be convertible to the \c Flow type
2.262 /// of the algorithm.
2.263 ///
2.264 /// \return <tt>(*this)</tt>
2.265 - template<typename UPPER>
2.266 - NetworkSimplex& upperMap(const UPPER& map) {
2.267 + template<typename UpperMap>
2.268 + NetworkSimplex& upperMap(const UpperMap& map) {
2.269 delete _pupper;
2.270 _pupper = new FlowArcMap(_graph);
2.271 for (ArcIt a(_graph); a != INVALID; ++a) {
2.272 @@ -731,43 +761,6 @@
2.273 return *this;
2.274 }
2.275
2.276 - /// \brief Set the upper bounds (capacities) on the arcs.
2.277 - ///
2.278 - /// This function sets the upper bounds (capacities) on the arcs.
2.279 - /// It is just an alias for \ref upperMap().
2.280 - ///
2.281 - /// \return <tt>(*this)</tt>
2.282 - template<typename CAP>
2.283 - NetworkSimplex& capacityMap(const CAP& map) {
2.284 - return upperMap(map);
2.285 - }
2.286 -
2.287 - /// \brief Set the lower and upper bounds on the arcs.
2.288 - ///
2.289 - /// This function sets the lower and upper bounds on the arcs.
2.290 - /// If neither this function nor \ref lowerMap() is used before
2.291 - /// calling \ref run(), the lower bounds will be set to zero
2.292 - /// on all arcs.
2.293 - /// If none of the functions \ref upperMap(), \ref capacityMap()
2.294 - /// and \ref boundMaps() is used before calling \ref run(),
2.295 - /// the upper bounds (capacities) will be set to
2.296 - /// \c std::numeric_limits<Flow>::max() on all arcs.
2.297 - ///
2.298 - /// \param lower An arc map storing the lower bounds.
2.299 - /// \param upper An arc map storing the upper bounds.
2.300 - ///
2.301 - /// The \c Value type of the maps must be convertible to the
2.302 - /// \c Flow type of the algorithm.
2.303 - ///
2.304 - /// \note This function is just a shortcut of calling \ref lowerMap()
2.305 - /// and \ref upperMap() separately.
2.306 - ///
2.307 - /// \return <tt>(*this)</tt>
2.308 - template <typename LOWER, typename UPPER>
2.309 - NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
2.310 - return lowerMap(lower).upperMap(upper);
2.311 - }
2.312 -
2.313 /// \brief Set the costs of the arcs.
2.314 ///
2.315 /// This function sets the costs of the arcs.
2.316 @@ -779,8 +772,8 @@
2.317 /// of the algorithm.
2.318 ///
2.319 /// \return <tt>(*this)</tt>
2.320 - template<typename COST>
2.321 - NetworkSimplex& costMap(const COST& map) {
2.322 + template<typename CostMap>
2.323 + NetworkSimplex& costMap(const CostMap& map) {
2.324 delete _pcost;
2.325 _pcost = new CostArcMap(_graph);
2.326 for (ArcIt a(_graph); a != INVALID; ++a) {
2.327 @@ -801,8 +794,8 @@
2.328 /// of the algorithm.
2.329 ///
2.330 /// \return <tt>(*this)</tt>
2.331 - template<typename SUP>
2.332 - NetworkSimplex& supplyMap(const SUP& map) {
2.333 + template<typename SupplyMap>
2.334 + NetworkSimplex& supplyMap(const SupplyMap& map) {
2.335 delete _psupply;
2.336 _pstsup = false;
2.337 _psupply = new FlowNodeMap(_graph);
2.338 @@ -820,6 +813,10 @@
2.339 /// calling \ref run(), the supply of each node will be set to zero.
2.340 /// (It makes sense only if non-zero lower bounds are given.)
2.341 ///
2.342 + /// Using this function has the same effect as using \ref supplyMap()
2.343 + /// with such a map in which \c k is assigned to \c s, \c -k is
2.344 + /// assigned to \c t and all other nodes have zero supply value.
2.345 + ///
2.346 /// \param s The source node.
2.347 /// \param t The target node.
2.348 /// \param k The required amount of flow from node \c s to node \c t
2.349 @@ -836,17 +833,17 @@
2.350 return *this;
2.351 }
2.352
2.353 - /// \brief Set the problem type.
2.354 + /// \brief Set the type of the supply constraints.
2.355 ///
2.356 - /// This function sets the problem type for the algorithm.
2.357 - /// If it is not used before calling \ref run(), the \ref GEQ problem
2.358 + /// This function sets the type of the supply/demand constraints.
2.359 + /// If it is not used before calling \ref run(), the \ref GEQ supply
2.360 /// type will be used.
2.361 ///
2.362 - /// For more information see \ref ProblemType.
2.363 + /// For more information see \ref SupplyType.
2.364 ///
2.365 /// \return <tt>(*this)</tt>
2.366 - NetworkSimplex& problemType(ProblemType problem_type) {
2.367 - _ptype = problem_type;
2.368 + NetworkSimplex& supplyType(SupplyType supply_type) {
2.369 + _stype = supply_type;
2.370 return *this;
2.371 }
2.372
2.373 @@ -896,13 +893,12 @@
2.374 ///
2.375 /// This function runs the algorithm.
2.376 /// The paramters can be specified using functions \ref lowerMap(),
2.377 - /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
2.378 - /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
2.379 - /// \ref problemType(), \ref flowMap() and \ref potentialMap().
2.380 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
2.381 + /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
2.382 /// For example,
2.383 /// \code
2.384 /// NetworkSimplex<ListDigraph> ns(graph);
2.385 - /// ns.boundMaps(lower, upper).costMap(cost)
2.386 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
2.387 /// .supplyMap(sup).run();
2.388 /// \endcode
2.389 ///
2.390 @@ -914,17 +910,25 @@
2.391 /// \param pivot_rule The pivot rule that will be used during the
2.392 /// algorithm. For more information see \ref PivotRule.
2.393 ///
2.394 - /// \return \c true if a feasible flow can be found.
2.395 - bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
2.396 - return init() && start(pivot_rule);
2.397 + /// \return \c INFEASIBLE if no feasible flow exists,
2.398 + /// \n \c OPTIMAL if the problem has optimal solution
2.399 + /// (i.e. it is feasible and bounded), and the algorithm has found
2.400 + /// optimal flow and node potentials (primal and dual solutions),
2.401 + /// \n \c UNBOUNDED if the objective function of the problem is
2.402 + /// unbounded, i.e. there is a directed cycle having negative total
2.403 + /// cost and infinite upper bound.
2.404 + ///
2.405 + /// \see ProblemType, PivotRule
2.406 + ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
2.407 + if (!init()) return INFEASIBLE;
2.408 + return start(pivot_rule);
2.409 }
2.410
2.411 /// \brief Reset all the parameters that have been given before.
2.412 ///
2.413 /// This function resets all the paramaters that have been given
2.414 /// before using functions \ref lowerMap(), \ref upperMap(),
2.415 - /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
2.416 - /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
2.417 + /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
2.418 /// \ref flowMap() and \ref potentialMap().
2.419 ///
2.420 /// It is useful for multiple run() calls. If this function is not
2.421 @@ -936,7 +940,7 @@
2.422 /// NetworkSimplex<ListDigraph> ns(graph);
2.423 ///
2.424 /// // First run
2.425 - /// ns.lowerMap(lower).capacityMap(cap).costMap(cost)
2.426 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
2.427 /// .supplyMap(sup).run();
2.428 ///
2.429 /// // Run again with modified cost map (reset() is not called,
2.430 @@ -947,7 +951,7 @@
2.431 /// // Run again from scratch using reset()
2.432 /// // (the lower bounds will be set to zero on all arcs)
2.433 /// ns.reset();
2.434 - /// ns.capacityMap(cap).costMap(cost)
2.435 + /// ns.upperMap(capacity).costMap(cost)
2.436 /// .supplyMap(sup).run();
2.437 /// \endcode
2.438 ///
2.439 @@ -962,7 +966,7 @@
2.440 _pcost = NULL;
2.441 _psupply = NULL;
2.442 _pstsup = false;
2.443 - _ptype = GEQ;
2.444 + _stype = GEQ;
2.445 if (_local_flow) delete _flow_map;
2.446 if (_local_potential) delete _potential_map;
2.447 _flow_map = NULL;
2.448 @@ -985,7 +989,7 @@
2.449 /// \brief Return the total cost of the found flow.
2.450 ///
2.451 /// This function returns the total cost of the found flow.
2.452 - /// The complexity of the function is O(e).
2.453 + /// Its complexity is O(e).
2.454 ///
2.455 /// \note The return type of the function can be specified as a
2.456 /// template parameter. For example,
2.457 @@ -997,9 +1001,9 @@
2.458 /// function.
2.459 ///
2.460 /// \pre \ref run() must be called before using this function.
2.461 - template <typename Num>
2.462 - Num totalCost() const {
2.463 - Num c = 0;
2.464 + template <typename Value>
2.465 + Value totalCost() const {
2.466 + Value c = 0;
2.467 if (_pcost) {
2.468 for (ArcIt e(_graph); e != INVALID; ++e)
2.469 c += (*_flow_map)[e] * (*_pcost)[e];
2.470 @@ -1050,7 +1054,7 @@
2.471 ///
2.472 /// This function returns a const reference to a node map storing
2.473 /// the found potentials, which form the dual solution of the
2.474 - /// \ref min_cost_flow "minimum cost flow" problem.
2.475 + /// \ref min_cost_flow "minimum cost flow problem".
2.476 ///
2.477 /// \pre \ref run() must be called before using this function.
2.478 const PotentialMap& potentialMap() const {
2.479 @@ -1101,7 +1105,7 @@
2.480
2.481 // Initialize node related data
2.482 bool valid_supply = true;
2.483 - Flow sum_supply = 0;
2.484 + _sum_supply = 0;
2.485 if (!_pstsup && !_psupply) {
2.486 _pstsup = true;
2.487 _psource = _ptarget = NodeIt(_graph);
2.488 @@ -1112,10 +1116,10 @@
2.489 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
2.490 _node_id[n] = i;
2.491 _supply[i] = (*_psupply)[n];
2.492 - sum_supply += _supply[i];
2.493 + _sum_supply += _supply[i];
2.494 }
2.495 - valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
2.496 - (_ptype == LEQ && sum_supply >= 0);
2.497 + valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
2.498 + (_stype == LEQ && _sum_supply >= 0);
2.499 } else {
2.500 int i = 0;
2.501 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
2.502 @@ -1127,92 +1131,18 @@
2.503 }
2.504 if (!valid_supply) return false;
2.505
2.506 - // Infinite capacity value
2.507 - Flow inf_cap =
2.508 - std::numeric_limits<Flow>::has_infinity ?
2.509 - std::numeric_limits<Flow>::infinity() :
2.510 - std::numeric_limits<Flow>::max();
2.511 -
2.512 // Initialize artifical cost
2.513 - Cost art_cost;
2.514 + Cost ART_COST;
2.515 if (std::numeric_limits<Cost>::is_exact) {
2.516 - art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
2.517 + ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
2.518 } else {
2.519 - art_cost = std::numeric_limits<Cost>::min();
2.520 + ART_COST = std::numeric_limits<Cost>::min();
2.521 for (int i = 0; i != _arc_num; ++i) {
2.522 - if (_cost[i] > art_cost) art_cost = _cost[i];
2.523 + if (_cost[i] > ART_COST) ART_COST = _cost[i];
2.524 }
2.525 - art_cost = (art_cost + 1) * _node_num;
2.526 + ART_COST = (ART_COST + 1) * _node_num;
2.527 }
2.528
2.529 - // Run Circulation to check if a feasible solution exists
2.530 - typedef ConstMap<Arc, Flow> ConstArcMap;
2.531 - ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
2.532 - FlowNodeMap *csup = NULL;
2.533 - bool local_csup = false;
2.534 - if (_psupply) {
2.535 - csup = _psupply;
2.536 - } else {
2.537 - csup = new FlowNodeMap(_graph, 0);
2.538 - (*csup)[_psource] = _pstflow;
2.539 - (*csup)[_ptarget] = -_pstflow;
2.540 - local_csup = true;
2.541 - }
2.542 - bool circ_result = false;
2.543 - if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
2.544 - // GEQ problem type
2.545 - if (_plower) {
2.546 - if (_pupper) {
2.547 - Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
2.548 - circ(_graph, *_plower, *_pupper, *csup);
2.549 - circ_result = circ.run();
2.550 - } else {
2.551 - Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
2.552 - circ(_graph, *_plower, inf_arc_map, *csup);
2.553 - circ_result = circ.run();
2.554 - }
2.555 - } else {
2.556 - if (_pupper) {
2.557 - Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
2.558 - circ(_graph, zero_arc_map, *_pupper, *csup);
2.559 - circ_result = circ.run();
2.560 - } else {
2.561 - Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
2.562 - circ(_graph, zero_arc_map, inf_arc_map, *csup);
2.563 - circ_result = circ.run();
2.564 - }
2.565 - }
2.566 - } else {
2.567 - // LEQ problem type
2.568 - typedef ReverseDigraph<const GR> RevGraph;
2.569 - typedef NegMap<FlowNodeMap> NegNodeMap;
2.570 - RevGraph rgraph(_graph);
2.571 - NegNodeMap neg_csup(*csup);
2.572 - if (_plower) {
2.573 - if (_pupper) {
2.574 - Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
2.575 - circ(rgraph, *_plower, *_pupper, neg_csup);
2.576 - circ_result = circ.run();
2.577 - } else {
2.578 - Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
2.579 - circ(rgraph, *_plower, inf_arc_map, neg_csup);
2.580 - circ_result = circ.run();
2.581 - }
2.582 - } else {
2.583 - if (_pupper) {
2.584 - Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
2.585 - circ(rgraph, zero_arc_map, *_pupper, neg_csup);
2.586 - circ_result = circ.run();
2.587 - } else {
2.588 - Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
2.589 - circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
2.590 - circ_result = circ.run();
2.591 - }
2.592 - }
2.593 - }
2.594 - if (local_csup) delete csup;
2.595 - if (!circ_result) return false;
2.596 -
2.597 // Set data for the artificial root node
2.598 _root = _node_num;
2.599 _parent[_root] = -1;
2.600 @@ -1221,11 +1151,11 @@
2.601 _rev_thread[0] = _root;
2.602 _succ_num[_root] = all_node_num;
2.603 _last_succ[_root] = _root - 1;
2.604 - _supply[_root] = -sum_supply;
2.605 - if (sum_supply < 0) {
2.606 - _pi[_root] = -art_cost;
2.607 + _supply[_root] = -_sum_supply;
2.608 + if (_sum_supply < 0) {
2.609 + _pi[_root] = -ART_COST;
2.610 } else {
2.611 - _pi[_root] = art_cost;
2.612 + _pi[_root] = ART_COST;
2.613 }
2.614
2.615 // Store the arcs in a mixed order
2.616 @@ -1260,7 +1190,7 @@
2.617 _cap[i] = (*_pupper)[_arc_ref[i]];
2.618 } else {
2.619 for (int i = 0; i != _arc_num; ++i)
2.620 - _cap[i] = inf_cap;
2.621 + _cap[i] = INF;
2.622 }
2.623 if (_pcost) {
2.624 for (int i = 0; i != _arc_num; ++i)
2.625 @@ -1275,8 +1205,17 @@
2.626 if (_plower) {
2.627 for (int i = 0; i != _arc_num; ++i) {
2.628 Flow c = (*_plower)[_arc_ref[i]];
2.629 - if (c != 0) {
2.630 - _cap[i] -= c;
2.631 + if (c > 0) {
2.632 + if (_cap[i] < INF) _cap[i] -= c;
2.633 + _supply[_source[i]] -= c;
2.634 + _supply[_target[i]] += c;
2.635 + }
2.636 + else if (c < 0) {
2.637 + if (_cap[i] < INF + c) {
2.638 + _cap[i] -= c;
2.639 + } else {
2.640 + _cap[i] = INF;
2.641 + }
2.642 _supply[_source[i]] -= c;
2.643 _supply[_target[i]] += c;
2.644 }
2.645 @@ -1291,17 +1230,17 @@
2.646 _last_succ[u] = u;
2.647 _parent[u] = _root;
2.648 _pred[u] = e;
2.649 - _cost[e] = art_cost;
2.650 - _cap[e] = inf_cap;
2.651 + _cost[e] = ART_COST;
2.652 + _cap[e] = INF;
2.653 _state[e] = STATE_TREE;
2.654 - if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
2.655 + if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
2.656 _flow[e] = _supply[u];
2.657 _forward[u] = true;
2.658 - _pi[u] = -art_cost + _pi[_root];
2.659 + _pi[u] = -ART_COST + _pi[_root];
2.660 } else {
2.661 _flow[e] = -_supply[u];
2.662 _forward[u] = false;
2.663 - _pi[u] = art_cost + _pi[_root];
2.664 + _pi[u] = ART_COST + _pi[_root];
2.665 }
2.666 }
2.667
2.668 @@ -1342,7 +1281,8 @@
2.669 // Search the cycle along the path form the first node to the root
2.670 for (int u = first; u != join; u = _parent[u]) {
2.671 e = _pred[u];
2.672 - d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
2.673 + d = _forward[u] ?
2.674 + _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
2.675 if (d < delta) {
2.676 delta = d;
2.677 u_out = u;
2.678 @@ -1352,7 +1292,8 @@
2.679 // Search the cycle along the path form the second node to the root
2.680 for (int u = second; u != join; u = _parent[u]) {
2.681 e = _pred[u];
2.682 - d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
2.683 + d = _forward[u] ?
2.684 + (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
2.685 if (d <= delta) {
2.686 delta = d;
2.687 u_out = u;
2.688 @@ -1526,7 +1467,7 @@
2.689 }
2.690
2.691 // Execute the algorithm
2.692 - bool start(PivotRule pivot_rule) {
2.693 + ProblemType start(PivotRule pivot_rule) {
2.694 // Select the pivot rule implementation
2.695 switch (pivot_rule) {
2.696 case FIRST_ELIGIBLE:
2.697 @@ -1540,23 +1481,41 @@
2.698 case ALTERING_LIST:
2.699 return start<AlteringListPivotRule>();
2.700 }
2.701 - return false;
2.702 + return INFEASIBLE; // avoid warning
2.703 }
2.704
2.705 template <typename PivotRuleImpl>
2.706 - bool start() {
2.707 + ProblemType start() {
2.708 PivotRuleImpl pivot(*this);
2.709
2.710 // Execute the Network Simplex algorithm
2.711 while (pivot.findEnteringArc()) {
2.712 findJoinNode();
2.713 bool change = findLeavingArc();
2.714 + if (delta >= INF) return UNBOUNDED;
2.715 changeFlow(change);
2.716 if (change) {
2.717 updateTreeStructure();
2.718 updatePotential();
2.719 }
2.720 }
2.721 +
2.722 + // Check feasibility
2.723 + if (_sum_supply < 0) {
2.724 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
2.725 + if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
2.726 + }
2.727 + }
2.728 + else if (_sum_supply > 0) {
2.729 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
2.730 + if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
2.731 + }
2.732 + }
2.733 + else {
2.734 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
2.735 + if (_flow[e] != 0) return INFEASIBLE;
2.736 + }
2.737 + }
2.738
2.739 // Copy flow values to _flow_map
2.740 if (_plower) {
2.741 @@ -1574,7 +1533,7 @@
2.742 _potential_map->set(n, _pi[_node_id[n]]);
2.743 }
2.744
2.745 - return true;
2.746 + return OPTIMAL;
2.747 }
2.748
2.749 }; //class NetworkSimplex
3.1 --- a/test/min_cost_flow_test.cc Sun Apr 26 16:36:23 2009 +0100
3.2 +++ b/test/min_cost_flow_test.cc Wed Apr 29 03:15:24 2009 +0200
3.3 @@ -18,6 +18,7 @@
3.4
3.5 #include <iostream>
3.6 #include <fstream>
3.7 +#include <limits>
3.8
3.9 #include <lemon/list_graph.h>
3.10 #include <lemon/lgf_reader.h>
3.11 @@ -33,50 +34,50 @@
3.12
3.13 char test_lgf[] =
3.14 "@nodes\n"
3.15 - "label sup1 sup2 sup3 sup4 sup5\n"
3.16 - " 1 20 27 0 20 30\n"
3.17 - " 2 -4 0 0 -8 -3\n"
3.18 - " 3 0 0 0 0 0\n"
3.19 - " 4 0 0 0 0 0\n"
3.20 - " 5 9 0 0 6 11\n"
3.21 - " 6 -6 0 0 -5 -6\n"
3.22 - " 7 0 0 0 0 0\n"
3.23 - " 8 0 0 0 0 3\n"
3.24 - " 9 3 0 0 0 0\n"
3.25 - " 10 -2 0 0 -7 -2\n"
3.26 - " 11 0 0 0 -10 0\n"
3.27 - " 12 -20 -27 0 -30 -20\n"
3.28 - "\n"
3.29 + "label sup1 sup2 sup3 sup4 sup5 sup6\n"
3.30 + " 1 20 27 0 30 20 30\n"
3.31 + " 2 -4 0 0 0 -8 -3\n"
3.32 + " 3 0 0 0 0 0 0\n"
3.33 + " 4 0 0 0 0 0 0\n"
3.34 + " 5 9 0 0 0 6 11\n"
3.35 + " 6 -6 0 0 0 -5 -6\n"
3.36 + " 7 0 0 0 0 0 0\n"
3.37 + " 8 0 0 0 0 0 3\n"
3.38 + " 9 3 0 0 0 0 0\n"
3.39 + " 10 -2 0 0 0 -7 -2\n"
3.40 + " 11 0 0 0 0 -10 0\n"
3.41 + " 12 -20 -27 0 -30 -30 -20\n"
3.42 + "\n"
3.43 "@arcs\n"
3.44 - " cost cap low1 low2\n"
3.45 - " 1 2 70 11 0 8\n"
3.46 - " 1 3 150 3 0 1\n"
3.47 - " 1 4 80 15 0 2\n"
3.48 - " 2 8 80 12 0 0\n"
3.49 - " 3 5 140 5 0 3\n"
3.50 - " 4 6 60 10 0 1\n"
3.51 - " 4 7 80 2 0 0\n"
3.52 - " 4 8 110 3 0 0\n"
3.53 - " 5 7 60 14 0 0\n"
3.54 - " 5 11 120 12 0 0\n"
3.55 - " 6 3 0 3 0 0\n"
3.56 - " 6 9 140 4 0 0\n"
3.57 - " 6 10 90 8 0 0\n"
3.58 - " 7 1 30 5 0 0\n"
3.59 - " 8 12 60 16 0 4\n"
3.60 - " 9 12 50 6 0 0\n"
3.61 - "10 12 70 13 0 5\n"
3.62 - "10 2 100 7 0 0\n"
3.63 - "10 7 60 10 0 0\n"
3.64 - "11 10 20 14 0 6\n"
3.65 - "12 11 30 10 0 0\n"
3.66 + " cost cap low1 low2 low3\n"
3.67 + " 1 2 70 11 0 8 8\n"
3.68 + " 1 3 150 3 0 1 0\n"
3.69 + " 1 4 80 15 0 2 2\n"
3.70 + " 2 8 80 12 0 0 0\n"
3.71 + " 3 5 140 5 0 3 1\n"
3.72 + " 4 6 60 10 0 1 0\n"
3.73 + " 4 7 80 2 0 0 0\n"
3.74 + " 4 8 110 3 0 0 0\n"
3.75 + " 5 7 60 14 0 0 0\n"
3.76 + " 5 11 120 12 0 0 0\n"
3.77 + " 6 3 0 3 0 0 0\n"
3.78 + " 6 9 140 4 0 0 0\n"
3.79 + " 6 10 90 8 0 0 0\n"
3.80 + " 7 1 30 5 0 0 -5\n"
3.81 + " 8 12 60 16 0 4 3\n"
3.82 + " 9 12 50 6 0 0 0\n"
3.83 + "10 12 70 13 0 5 2\n"
3.84 + "10 2 100 7 0 0 0\n"
3.85 + "10 7 60 10 0 0 -3\n"
3.86 + "11 10 20 14 0 6 -20\n"
3.87 + "12 11 30 10 0 0 -10\n"
3.88 "\n"
3.89 "@attributes\n"
3.90 "source 1\n"
3.91 "target 12\n";
3.92
3.93
3.94 -enum ProblemType {
3.95 +enum SupplyType {
3.96 EQ,
3.97 GEQ,
3.98 LEQ
3.99 @@ -98,8 +99,6 @@
3.100 b = mcf.reset()
3.101 .lowerMap(lower)
3.102 .upperMap(upper)
3.103 - .capacityMap(upper)
3.104 - .boundMaps(lower, upper)
3.105 .costMap(cost)
3.106 .supplyMap(sup)
3.107 .stSupply(n, n, k)
3.108 @@ -112,10 +111,12 @@
3.109 const typename MCF::FlowMap &fm = const_mcf.flowMap();
3.110 const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
3.111
3.112 - v = const_mcf.totalCost();
3.113 + c = const_mcf.totalCost();
3.114 double x = const_mcf.template totalCost<double>();
3.115 v = const_mcf.flow(a);
3.116 - v = const_mcf.potential(n);
3.117 + c = const_mcf.potential(n);
3.118 +
3.119 + v = const_mcf.INF;
3.120
3.121 ignore_unused_variable_warning(fm);
3.122 ignore_unused_variable_warning(pm);
3.123 @@ -137,6 +138,7 @@
3.124 const Arc &a;
3.125 const Flow &k;
3.126 Flow v;
3.127 + Cost c;
3.128 bool b;
3.129
3.130 typename MCF::FlowMap &flow;
3.131 @@ -151,7 +153,7 @@
3.132 typename SM, typename FM >
3.133 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
3.134 const SM& supply, const FM& flow,
3.135 - ProblemType type = EQ )
3.136 + SupplyType type = EQ )
3.137 {
3.138 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
3.139
3.140 @@ -208,16 +210,17 @@
3.141 // Run a minimum cost flow algorithm and check the results
3.142 template < typename MCF, typename GR,
3.143 typename LM, typename UM,
3.144 - typename CM, typename SM >
3.145 -void checkMcf( const MCF& mcf, bool mcf_result,
3.146 + typename CM, typename SM,
3.147 + typename PT >
3.148 +void checkMcf( const MCF& mcf, PT mcf_result,
3.149 const GR& gr, const LM& lower, const UM& upper,
3.150 const CM& cost, const SM& supply,
3.151 - bool result, typename CM::Value total,
3.152 + PT result, bool optimal, typename CM::Value total,
3.153 const std::string &test_id = "",
3.154 - ProblemType type = EQ )
3.155 + SupplyType type = EQ )
3.156 {
3.157 check(mcf_result == result, "Wrong result " + test_id);
3.158 - if (result) {
3.159 + if (optimal) {
3.160 check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
3.161 "The flow is not feasible " + test_id);
3.162 check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
3.163 @@ -244,8 +247,8 @@
3.164
3.165 // Read the test digraph
3.166 Digraph gr;
3.167 - Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
3.168 - Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
3.169 + Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
3.170 + Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
3.171 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
3.172 Node v, w;
3.173
3.174 @@ -255,14 +258,56 @@
3.175 .arcMap("cap", u)
3.176 .arcMap("low1", l1)
3.177 .arcMap("low2", l2)
3.178 + .arcMap("low3", l3)
3.179 .nodeMap("sup1", s1)
3.180 .nodeMap("sup2", s2)
3.181 .nodeMap("sup3", s3)
3.182 .nodeMap("sup4", s4)
3.183 .nodeMap("sup5", s5)
3.184 + .nodeMap("sup6", s6)
3.185 .node("source", v)
3.186 .node("target", w)
3.187 .run();
3.188 +
3.189 + // Build a test digraph for testing negative costs
3.190 + Digraph ngr;
3.191 + Node n1 = ngr.addNode();
3.192 + Node n2 = ngr.addNode();
3.193 + Node n3 = ngr.addNode();
3.194 + Node n4 = ngr.addNode();
3.195 + Node n5 = ngr.addNode();
3.196 + Node n6 = ngr.addNode();
3.197 + Node n7 = ngr.addNode();
3.198 +
3.199 + Arc a1 = ngr.addArc(n1, n2);
3.200 + Arc a2 = ngr.addArc(n1, n3);
3.201 + Arc a3 = ngr.addArc(n2, n4);
3.202 + Arc a4 = ngr.addArc(n3, n4);
3.203 + Arc a5 = ngr.addArc(n3, n2);
3.204 + Arc a6 = ngr.addArc(n5, n3);
3.205 + Arc a7 = ngr.addArc(n5, n6);
3.206 + Arc a8 = ngr.addArc(n6, n7);
3.207 + Arc a9 = ngr.addArc(n7, n5);
3.208 +
3.209 + Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
3.210 + ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
3.211 + Digraph::NodeMap<int> ns(ngr, 0);
3.212 +
3.213 + nl2[a7] = 1000;
3.214 + nl2[a8] = -1000;
3.215 +
3.216 + ns[n1] = 100;
3.217 + ns[n4] = -100;
3.218 +
3.219 + nc[a1] = 100;
3.220 + nc[a2] = 30;
3.221 + nc[a3] = 20;
3.222 + nc[a4] = 80;
3.223 + nc[a5] = 50;
3.224 + nc[a6] = 10;
3.225 + nc[a7] = 80;
3.226 + nc[a8] = 30;
3.227 + nc[a9] = -120;
3.228
3.229 // A. Test NetworkSimplex with the default pivot rule
3.230 {
3.231 @@ -271,63 +316,77 @@
3.232 // Check the equality form
3.233 mcf.upperMap(u).costMap(c);
3.234 checkMcf(mcf, mcf.supplyMap(s1).run(),
3.235 - gr, l1, u, c, s1, true, 5240, "#A1");
3.236 + gr, l1, u, c, s1, mcf.OPTIMAL, true, 5240, "#A1");
3.237 checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
3.238 - gr, l1, u, c, s2, true, 7620, "#A2");
3.239 + gr, l1, u, c, s2, mcf.OPTIMAL, true, 7620, "#A2");
3.240 mcf.lowerMap(l2);
3.241 checkMcf(mcf, mcf.supplyMap(s1).run(),
3.242 - gr, l2, u, c, s1, true, 5970, "#A3");
3.243 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#A3");
3.244 checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
3.245 - gr, l2, u, c, s2, true, 8010, "#A4");
3.246 + gr, l2, u, c, s2, mcf.OPTIMAL, true, 8010, "#A4");
3.247 mcf.reset();
3.248 checkMcf(mcf, mcf.supplyMap(s1).run(),
3.249 - gr, l1, cu, cc, s1, true, 74, "#A5");
3.250 + gr, l1, cu, cc, s1, mcf.OPTIMAL, true, 74, "#A5");
3.251 checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
3.252 - gr, l2, cu, cc, s2, true, 94, "#A6");
3.253 + gr, l2, cu, cc, s2, mcf.OPTIMAL, true, 94, "#A6");
3.254 mcf.reset();
3.255 checkMcf(mcf, mcf.run(),
3.256 - gr, l1, cu, cc, s3, true, 0, "#A7");
3.257 - checkMcf(mcf, mcf.boundMaps(l2, u).run(),
3.258 - gr, l2, u, cc, s3, false, 0, "#A8");
3.259 + gr, l1, cu, cc, s3, mcf.OPTIMAL, true, 0, "#A7");
3.260 + checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
3.261 + gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
3.262 + mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
3.263 + checkMcf(mcf, mcf.run(),
3.264 + gr, l3, u, c, s4, mcf.OPTIMAL, true, 6360, "#A9");
3.265
3.266 // Check the GEQ form
3.267 - mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
3.268 + mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
3.269 checkMcf(mcf, mcf.run(),
3.270 - gr, l1, u, c, s4, true, 3530, "#A9", GEQ);
3.271 - mcf.problemType(mcf.GEQ);
3.272 + gr, l1, u, c, s5, mcf.OPTIMAL, true, 3530, "#A10", GEQ);
3.273 + mcf.supplyType(mcf.GEQ);
3.274 checkMcf(mcf, mcf.lowerMap(l2).run(),
3.275 - gr, l2, u, c, s4, true, 4540, "#A10", GEQ);
3.276 - mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
3.277 + gr, l2, u, c, s5, mcf.OPTIMAL, true, 4540, "#A11", GEQ);
3.278 + mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
3.279 checkMcf(mcf, mcf.run(),
3.280 - gr, l2, u, c, s5, false, 0, "#A11", GEQ);
3.281 + gr, l2, u, c, s6, mcf.INFEASIBLE, false, 0, "#A12", GEQ);
3.282
3.283 // Check the LEQ form
3.284 - mcf.reset().problemType(mcf.LEQ);
3.285 - mcf.upperMap(u).costMap(c).supplyMap(s5);
3.286 + mcf.reset().supplyType(mcf.LEQ);
3.287 + mcf.upperMap(u).costMap(c).supplyMap(s6);
3.288 checkMcf(mcf, mcf.run(),
3.289 - gr, l1, u, c, s5, true, 5080, "#A12", LEQ);
3.290 + gr, l1, u, c, s6, mcf.OPTIMAL, true, 5080, "#A13", LEQ);
3.291 checkMcf(mcf, mcf.lowerMap(l2).run(),
3.292 - gr, l2, u, c, s5, true, 5930, "#A13", LEQ);
3.293 - mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
3.294 + gr, l2, u, c, s6, mcf.OPTIMAL, true, 5930, "#A14", LEQ);
3.295 + mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
3.296 checkMcf(mcf, mcf.run(),
3.297 - gr, l2, u, c, s4, false, 0, "#A14", LEQ);
3.298 + gr, l2, u, c, s5, mcf.INFEASIBLE, false, 0, "#A15", LEQ);
3.299 +
3.300 + // Check negative costs
3.301 + NetworkSimplex<Digraph> nmcf(ngr);
3.302 + nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
3.303 + checkMcf(nmcf, nmcf.run(),
3.304 + ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A16");
3.305 + checkMcf(nmcf, nmcf.upperMap(nu2).run(),
3.306 + ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true, -40000, "#A17");
3.307 + nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
3.308 + checkMcf(nmcf, nmcf.run(),
3.309 + ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A18");
3.310 }
3.311
3.312 // B. Test NetworkSimplex with each pivot rule
3.313 {
3.314 NetworkSimplex<Digraph> mcf(gr);
3.315 - mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
3.316 + mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
3.317
3.318 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
3.319 - gr, l2, u, c, s1, true, 5970, "#B1");
3.320 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B1");
3.321 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
3.322 - gr, l2, u, c, s1, true, 5970, "#B2");
3.323 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B2");
3.324 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
3.325 - gr, l2, u, c, s1, true, 5970, "#B3");
3.326 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B3");
3.327 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
3.328 - gr, l2, u, c, s1, true, 5970, "#B4");
3.329 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B4");
3.330 checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
3.331 - gr, l2, u, c, s1, true, 5970, "#B5");
3.332 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B5");
3.333 }
3.334
3.335 return 0;
4.1 --- a/tools/dimacs-solver.cc Sun Apr 26 16:36:23 2009 +0100
4.2 +++ b/tools/dimacs-solver.cc Wed Apr 29 03:15:24 2009 +0200
4.3 @@ -119,8 +119,8 @@
4.4
4.5 ti.restart();
4.6 NetworkSimplex<Digraph, Value> ns(g);
4.7 - ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
4.8 - if (sum_sup > 0) ns.problemType(ns.LEQ);
4.9 + ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
4.10 + if (sum_sup > 0) ns.supplyType(ns.LEQ);
4.11 if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
4.12 ti.restart();
4.13 bool res = ns.run();