1.1 --- a/lemon/network_simplex.h Sun Apr 26 16:36:23 2009 +0100
1.2 +++ b/lemon/network_simplex.h Wed Apr 29 03:15:24 2009 +0200
1.3 @@ -30,9 +30,6 @@
1.4
1.5 #include <lemon/core.h>
1.6 #include <lemon/math.h>
1.7 -#include <lemon/maps.h>
1.8 -#include <lemon/circulation.h>
1.9 -#include <lemon/adaptors.h>
1.10
1.11 namespace lemon {
1.12
1.13 @@ -50,8 +47,13 @@
1.14 ///
1.15 /// In general this class is the fastest implementation available
1.16 /// in LEMON for the minimum cost flow problem.
1.17 - /// Moreover it supports both direction of the supply/demand inequality
1.18 - /// constraints. For more information see \ref ProblemType.
1.19 + /// Moreover it supports both directions of the supply/demand inequality
1.20 + /// constraints. For more information see \ref SupplyType.
1.21 + ///
1.22 + /// Most of the parameters of the problem (except for the digraph)
1.23 + /// can be given using separate functions, and the algorithm can be
1.24 + /// executed using the \ref run() function. If some parameters are not
1.25 + /// specified, then default values will be used.
1.26 ///
1.27 /// \tparam GR The digraph type the algorithm runs on.
1.28 /// \tparam F The value type used for flow amounts, capacity bounds
1.29 @@ -88,11 +90,80 @@
1.30
1.31 public:
1.32
1.33 - /// \brief Enum type for selecting the pivot rule.
1.34 + /// \brief Problem type constants for the \c run() function.
1.35 ///
1.36 - /// Enum type for selecting the pivot rule for the \ref run()
1.37 + /// Enum type containing the problem type constants that can be
1.38 + /// returned by the \ref run() function of the algorithm.
1.39 + enum ProblemType {
1.40 + /// The problem has no feasible solution (flow).
1.41 + INFEASIBLE,
1.42 + /// The problem has optimal solution (i.e. it is feasible and
1.43 + /// bounded), and the algorithm has found optimal flow and node
1.44 + /// potentials (primal and dual solutions).
1.45 + OPTIMAL,
1.46 + /// The objective function of the problem is unbounded, i.e.
1.47 + /// there is a directed cycle having negative total cost and
1.48 + /// infinite upper bound.
1.49 + UNBOUNDED
1.50 + };
1.51 +
1.52 + /// \brief Constants for selecting the type of the supply constraints.
1.53 + ///
1.54 + /// Enum type containing constants for selecting the supply type,
1.55 + /// i.e. the direction of the inequalities in the supply/demand
1.56 + /// constraints of the \ref min_cost_flow "minimum cost flow problem".
1.57 + ///
1.58 + /// The default supply type is \c GEQ, since this form is supported
1.59 + /// by other minimum cost flow algorithms and the \ref Circulation
1.60 + /// algorithm, as well.
1.61 + /// The \c LEQ problem type can be selected using the \ref supplyType()
1.62 /// function.
1.63 ///
1.64 + /// Note that the equality form is a special case of both supply types.
1.65 + enum SupplyType {
1.66 +
1.67 + /// This option means that there are <em>"greater or equal"</em>
1.68 + /// supply/demand constraints in the definition, i.e. the exact
1.69 + /// formulation of the problem is the following.
1.70 + /**
1.71 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.72 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
1.73 + sup(u) \quad \forall u\in V \f]
1.74 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.75 + */
1.76 + /// It means that the total demand must be greater or equal to the
1.77 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.78 + /// negative) and all the supplies have to be carried out from
1.79 + /// the supply nodes, but there could be demands that are not
1.80 + /// satisfied.
1.81 + GEQ,
1.82 + /// It is just an alias for the \c GEQ option.
1.83 + CARRY_SUPPLIES = GEQ,
1.84 +
1.85 + /// This option means that there are <em>"less or equal"</em>
1.86 + /// supply/demand constraints in the definition, i.e. the exact
1.87 + /// formulation of the problem is the following.
1.88 + /**
1.89 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.90 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
1.91 + sup(u) \quad \forall u\in V \f]
1.92 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.93 + */
1.94 + /// It means that the total demand must be less or equal to the
1.95 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.96 + /// positive) and all the demands have to be satisfied, but there
1.97 + /// could be supplies that are not carried out from the supply
1.98 + /// nodes.
1.99 + LEQ,
1.100 + /// It is just an alias for the \c LEQ option.
1.101 + SATISFY_DEMANDS = LEQ
1.102 + };
1.103 +
1.104 + /// \brief Constants for selecting the pivot rule.
1.105 + ///
1.106 + /// Enum type containing constants for selecting the pivot rule for
1.107 + /// the \ref run() function.
1.108 + ///
1.109 /// \ref NetworkSimplex provides five different pivot rule
1.110 /// implementations that significantly affect the running time
1.111 /// of the algorithm.
1.112 @@ -131,58 +202,6 @@
1.113 ALTERING_LIST
1.114 };
1.115
1.116 - /// \brief Enum type for selecting the problem type.
1.117 - ///
1.118 - /// Enum type for selecting the problem type, i.e. the direction of
1.119 - /// the inequalities in the supply/demand constraints of the
1.120 - /// \ref min_cost_flow "minimum cost flow problem".
1.121 - ///
1.122 - /// The default problem type is \c GEQ, since this form is supported
1.123 - /// by other minimum cost flow algorithms and the \ref Circulation
1.124 - /// algorithm as well.
1.125 - /// The \c LEQ problem type can be selected using the \ref problemType()
1.126 - /// function.
1.127 - ///
1.128 - /// Note that the equality form is a special case of both problem type.
1.129 - enum ProblemType {
1.130 -
1.131 - /// This option means that there are "<em>greater or equal</em>"
1.132 - /// constraints in the defintion, i.e. the exact formulation of the
1.133 - /// problem is the following.
1.134 - /**
1.135 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.136 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
1.137 - sup(u) \quad \forall u\in V \f]
1.138 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.139 - */
1.140 - /// It means that the total demand must be greater or equal to the
1.141 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.142 - /// negative) and all the supplies have to be carried out from
1.143 - /// the supply nodes, but there could be demands that are not
1.144 - /// satisfied.
1.145 - GEQ,
1.146 - /// It is just an alias for the \c GEQ option.
1.147 - CARRY_SUPPLIES = GEQ,
1.148 -
1.149 - /// This option means that there are "<em>less or equal</em>"
1.150 - /// constraints in the defintion, i.e. the exact formulation of the
1.151 - /// problem is the following.
1.152 - /**
1.153 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.154 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
1.155 - sup(u) \quad \forall u\in V \f]
1.156 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.157 - */
1.158 - /// It means that the total demand must be less or equal to the
1.159 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.160 - /// positive) and all the demands have to be satisfied, but there
1.161 - /// could be supplies that are not carried out from the supply
1.162 - /// nodes.
1.163 - LEQ,
1.164 - /// It is just an alias for the \c LEQ option.
1.165 - SATISFY_DEMANDS = LEQ
1.166 - };
1.167 -
1.168 private:
1.169
1.170 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.171 @@ -220,7 +239,9 @@
1.172 bool _pstsup;
1.173 Node _psource, _ptarget;
1.174 Flow _pstflow;
1.175 - ProblemType _ptype;
1.176 + SupplyType _stype;
1.177 +
1.178 + Flow _sum_supply;
1.179
1.180 // Result maps
1.181 FlowMap *_flow_map;
1.182 @@ -259,6 +280,15 @@
1.183 int stem, par_stem, new_stem;
1.184 Flow delta;
1.185
1.186 + public:
1.187 +
1.188 + /// \brief Constant for infinite upper bounds (capacities).
1.189 + ///
1.190 + /// Constant for infinite upper bounds (capacities).
1.191 + /// It is \c std::numeric_limits<Flow>::infinity() if available,
1.192 + /// \c std::numeric_limits<Flow>::max() otherwise.
1.193 + const Flow INF;
1.194 +
1.195 private:
1.196
1.197 // Implementation of the First Eligible pivot rule
1.198 @@ -661,17 +691,19 @@
1.199 NetworkSimplex(const GR& graph) :
1.200 _graph(graph),
1.201 _plower(NULL), _pupper(NULL), _pcost(NULL),
1.202 - _psupply(NULL), _pstsup(false), _ptype(GEQ),
1.203 + _psupply(NULL), _pstsup(false), _stype(GEQ),
1.204 _flow_map(NULL), _potential_map(NULL),
1.205 _local_flow(false), _local_potential(false),
1.206 - _node_id(graph)
1.207 + _node_id(graph),
1.208 + INF(std::numeric_limits<Flow>::has_infinity ?
1.209 + std::numeric_limits<Flow>::infinity() :
1.210 + std::numeric_limits<Flow>::max())
1.211 {
1.212 - LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
1.213 - std::numeric_limits<Flow>::is_signed,
1.214 - "The flow type of NetworkSimplex must be signed integer");
1.215 - LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
1.216 - std::numeric_limits<Cost>::is_signed,
1.217 - "The cost type of NetworkSimplex must be signed integer");
1.218 + // Check the value types
1.219 + LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,
1.220 + "The flow type of NetworkSimplex must be signed");
1.221 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.222 + "The cost type of NetworkSimplex must be signed");
1.223 }
1.224
1.225 /// Destructor.
1.226 @@ -689,17 +721,16 @@
1.227 /// \brief Set the lower bounds on the arcs.
1.228 ///
1.229 /// This function sets the lower bounds on the arcs.
1.230 - /// If neither this function nor \ref boundMaps() is used before
1.231 - /// calling \ref run(), the lower bounds will be set to zero
1.232 - /// on all arcs.
1.233 + /// If it is not used before calling \ref run(), the lower bounds
1.234 + /// will be set to zero on all arcs.
1.235 ///
1.236 /// \param map An arc map storing the lower bounds.
1.237 /// Its \c Value type must be convertible to the \c Flow type
1.238 /// of the algorithm.
1.239 ///
1.240 /// \return <tt>(*this)</tt>
1.241 - template <typename LOWER>
1.242 - NetworkSimplex& lowerMap(const LOWER& map) {
1.243 + template <typename LowerMap>
1.244 + NetworkSimplex& lowerMap(const LowerMap& map) {
1.245 delete _plower;
1.246 _plower = new FlowArcMap(_graph);
1.247 for (ArcIt a(_graph); a != INVALID; ++a) {
1.248 @@ -711,18 +742,17 @@
1.249 /// \brief Set the upper bounds (capacities) on the arcs.
1.250 ///
1.251 /// This function sets the upper bounds (capacities) on the arcs.
1.252 - /// If none of the functions \ref upperMap(), \ref capacityMap()
1.253 - /// and \ref boundMaps() is used before calling \ref run(),
1.254 - /// the upper bounds (capacities) will be set to
1.255 - /// \c std::numeric_limits<Flow>::max() on all arcs.
1.256 + /// If it is not used before calling \ref run(), the upper bounds
1.257 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.258 + /// unbounded from above on each arc).
1.259 ///
1.260 /// \param map An arc map storing the upper bounds.
1.261 /// Its \c Value type must be convertible to the \c Flow type
1.262 /// of the algorithm.
1.263 ///
1.264 /// \return <tt>(*this)</tt>
1.265 - template<typename UPPER>
1.266 - NetworkSimplex& upperMap(const UPPER& map) {
1.267 + template<typename UpperMap>
1.268 + NetworkSimplex& upperMap(const UpperMap& map) {
1.269 delete _pupper;
1.270 _pupper = new FlowArcMap(_graph);
1.271 for (ArcIt a(_graph); a != INVALID; ++a) {
1.272 @@ -731,43 +761,6 @@
1.273 return *this;
1.274 }
1.275
1.276 - /// \brief Set the upper bounds (capacities) on the arcs.
1.277 - ///
1.278 - /// This function sets the upper bounds (capacities) on the arcs.
1.279 - /// It is just an alias for \ref upperMap().
1.280 - ///
1.281 - /// \return <tt>(*this)</tt>
1.282 - template<typename CAP>
1.283 - NetworkSimplex& capacityMap(const CAP& map) {
1.284 - return upperMap(map);
1.285 - }
1.286 -
1.287 - /// \brief Set the lower and upper bounds on the arcs.
1.288 - ///
1.289 - /// This function sets the lower and upper bounds on the arcs.
1.290 - /// If neither this function nor \ref lowerMap() is used before
1.291 - /// calling \ref run(), the lower bounds will be set to zero
1.292 - /// on all arcs.
1.293 - /// If none of the functions \ref upperMap(), \ref capacityMap()
1.294 - /// and \ref boundMaps() is used before calling \ref run(),
1.295 - /// the upper bounds (capacities) will be set to
1.296 - /// \c std::numeric_limits<Flow>::max() on all arcs.
1.297 - ///
1.298 - /// \param lower An arc map storing the lower bounds.
1.299 - /// \param upper An arc map storing the upper bounds.
1.300 - ///
1.301 - /// The \c Value type of the maps must be convertible to the
1.302 - /// \c Flow type of the algorithm.
1.303 - ///
1.304 - /// \note This function is just a shortcut of calling \ref lowerMap()
1.305 - /// and \ref upperMap() separately.
1.306 - ///
1.307 - /// \return <tt>(*this)</tt>
1.308 - template <typename LOWER, typename UPPER>
1.309 - NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
1.310 - return lowerMap(lower).upperMap(upper);
1.311 - }
1.312 -
1.313 /// \brief Set the costs of the arcs.
1.314 ///
1.315 /// This function sets the costs of the arcs.
1.316 @@ -779,8 +772,8 @@
1.317 /// of the algorithm.
1.318 ///
1.319 /// \return <tt>(*this)</tt>
1.320 - template<typename COST>
1.321 - NetworkSimplex& costMap(const COST& map) {
1.322 + template<typename CostMap>
1.323 + NetworkSimplex& costMap(const CostMap& map) {
1.324 delete _pcost;
1.325 _pcost = new CostArcMap(_graph);
1.326 for (ArcIt a(_graph); a != INVALID; ++a) {
1.327 @@ -801,8 +794,8 @@
1.328 /// of the algorithm.
1.329 ///
1.330 /// \return <tt>(*this)</tt>
1.331 - template<typename SUP>
1.332 - NetworkSimplex& supplyMap(const SUP& map) {
1.333 + template<typename SupplyMap>
1.334 + NetworkSimplex& supplyMap(const SupplyMap& map) {
1.335 delete _psupply;
1.336 _pstsup = false;
1.337 _psupply = new FlowNodeMap(_graph);
1.338 @@ -820,6 +813,10 @@
1.339 /// calling \ref run(), the supply of each node will be set to zero.
1.340 /// (It makes sense only if non-zero lower bounds are given.)
1.341 ///
1.342 + /// Using this function has the same effect as using \ref supplyMap()
1.343 + /// with such a map in which \c k is assigned to \c s, \c -k is
1.344 + /// assigned to \c t and all other nodes have zero supply value.
1.345 + ///
1.346 /// \param s The source node.
1.347 /// \param t The target node.
1.348 /// \param k The required amount of flow from node \c s to node \c t
1.349 @@ -836,17 +833,17 @@
1.350 return *this;
1.351 }
1.352
1.353 - /// \brief Set the problem type.
1.354 + /// \brief Set the type of the supply constraints.
1.355 ///
1.356 - /// This function sets the problem type for the algorithm.
1.357 - /// If it is not used before calling \ref run(), the \ref GEQ problem
1.358 + /// This function sets the type of the supply/demand constraints.
1.359 + /// If it is not used before calling \ref run(), the \ref GEQ supply
1.360 /// type will be used.
1.361 ///
1.362 - /// For more information see \ref ProblemType.
1.363 + /// For more information see \ref SupplyType.
1.364 ///
1.365 /// \return <tt>(*this)</tt>
1.366 - NetworkSimplex& problemType(ProblemType problem_type) {
1.367 - _ptype = problem_type;
1.368 + NetworkSimplex& supplyType(SupplyType supply_type) {
1.369 + _stype = supply_type;
1.370 return *this;
1.371 }
1.372
1.373 @@ -896,13 +893,12 @@
1.374 ///
1.375 /// This function runs the algorithm.
1.376 /// The paramters can be specified using functions \ref lowerMap(),
1.377 - /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
1.378 - /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.379 - /// \ref problemType(), \ref flowMap() and \ref potentialMap().
1.380 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.381 + /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
1.382 /// For example,
1.383 /// \code
1.384 /// NetworkSimplex<ListDigraph> ns(graph);
1.385 - /// ns.boundMaps(lower, upper).costMap(cost)
1.386 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
1.387 /// .supplyMap(sup).run();
1.388 /// \endcode
1.389 ///
1.390 @@ -914,17 +910,25 @@
1.391 /// \param pivot_rule The pivot rule that will be used during the
1.392 /// algorithm. For more information see \ref PivotRule.
1.393 ///
1.394 - /// \return \c true if a feasible flow can be found.
1.395 - bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
1.396 - return init() && start(pivot_rule);
1.397 + /// \return \c INFEASIBLE if no feasible flow exists,
1.398 + /// \n \c OPTIMAL if the problem has optimal solution
1.399 + /// (i.e. it is feasible and bounded), and the algorithm has found
1.400 + /// optimal flow and node potentials (primal and dual solutions),
1.401 + /// \n \c UNBOUNDED if the objective function of the problem is
1.402 + /// unbounded, i.e. there is a directed cycle having negative total
1.403 + /// cost and infinite upper bound.
1.404 + ///
1.405 + /// \see ProblemType, PivotRule
1.406 + ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
1.407 + if (!init()) return INFEASIBLE;
1.408 + return start(pivot_rule);
1.409 }
1.410
1.411 /// \brief Reset all the parameters that have been given before.
1.412 ///
1.413 /// This function resets all the paramaters that have been given
1.414 /// before using functions \ref lowerMap(), \ref upperMap(),
1.415 - /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
1.416 - /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
1.417 + /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
1.418 /// \ref flowMap() and \ref potentialMap().
1.419 ///
1.420 /// It is useful for multiple run() calls. If this function is not
1.421 @@ -936,7 +940,7 @@
1.422 /// NetworkSimplex<ListDigraph> ns(graph);
1.423 ///
1.424 /// // First run
1.425 - /// ns.lowerMap(lower).capacityMap(cap).costMap(cost)
1.426 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
1.427 /// .supplyMap(sup).run();
1.428 ///
1.429 /// // Run again with modified cost map (reset() is not called,
1.430 @@ -947,7 +951,7 @@
1.431 /// // Run again from scratch using reset()
1.432 /// // (the lower bounds will be set to zero on all arcs)
1.433 /// ns.reset();
1.434 - /// ns.capacityMap(cap).costMap(cost)
1.435 + /// ns.upperMap(capacity).costMap(cost)
1.436 /// .supplyMap(sup).run();
1.437 /// \endcode
1.438 ///
1.439 @@ -962,7 +966,7 @@
1.440 _pcost = NULL;
1.441 _psupply = NULL;
1.442 _pstsup = false;
1.443 - _ptype = GEQ;
1.444 + _stype = GEQ;
1.445 if (_local_flow) delete _flow_map;
1.446 if (_local_potential) delete _potential_map;
1.447 _flow_map = NULL;
1.448 @@ -985,7 +989,7 @@
1.449 /// \brief Return the total cost of the found flow.
1.450 ///
1.451 /// This function returns the total cost of the found flow.
1.452 - /// The complexity of the function is O(e).
1.453 + /// Its complexity is O(e).
1.454 ///
1.455 /// \note The return type of the function can be specified as a
1.456 /// template parameter. For example,
1.457 @@ -997,9 +1001,9 @@
1.458 /// function.
1.459 ///
1.460 /// \pre \ref run() must be called before using this function.
1.461 - template <typename Num>
1.462 - Num totalCost() const {
1.463 - Num c = 0;
1.464 + template <typename Value>
1.465 + Value totalCost() const {
1.466 + Value c = 0;
1.467 if (_pcost) {
1.468 for (ArcIt e(_graph); e != INVALID; ++e)
1.469 c += (*_flow_map)[e] * (*_pcost)[e];
1.470 @@ -1050,7 +1054,7 @@
1.471 ///
1.472 /// This function returns a const reference to a node map storing
1.473 /// the found potentials, which form the dual solution of the
1.474 - /// \ref min_cost_flow "minimum cost flow" problem.
1.475 + /// \ref min_cost_flow "minimum cost flow problem".
1.476 ///
1.477 /// \pre \ref run() must be called before using this function.
1.478 const PotentialMap& potentialMap() const {
1.479 @@ -1101,7 +1105,7 @@
1.480
1.481 // Initialize node related data
1.482 bool valid_supply = true;
1.483 - Flow sum_supply = 0;
1.484 + _sum_supply = 0;
1.485 if (!_pstsup && !_psupply) {
1.486 _pstsup = true;
1.487 _psource = _ptarget = NodeIt(_graph);
1.488 @@ -1112,10 +1116,10 @@
1.489 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.490 _node_id[n] = i;
1.491 _supply[i] = (*_psupply)[n];
1.492 - sum_supply += _supply[i];
1.493 + _sum_supply += _supply[i];
1.494 }
1.495 - valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
1.496 - (_ptype == LEQ && sum_supply >= 0);
1.497 + valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
1.498 + (_stype == LEQ && _sum_supply >= 0);
1.499 } else {
1.500 int i = 0;
1.501 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.502 @@ -1127,92 +1131,18 @@
1.503 }
1.504 if (!valid_supply) return false;
1.505
1.506 - // Infinite capacity value
1.507 - Flow inf_cap =
1.508 - std::numeric_limits<Flow>::has_infinity ?
1.509 - std::numeric_limits<Flow>::infinity() :
1.510 - std::numeric_limits<Flow>::max();
1.511 -
1.512 // Initialize artifical cost
1.513 - Cost art_cost;
1.514 + Cost ART_COST;
1.515 if (std::numeric_limits<Cost>::is_exact) {
1.516 - art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1.517 + ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
1.518 } else {
1.519 - art_cost = std::numeric_limits<Cost>::min();
1.520 + ART_COST = std::numeric_limits<Cost>::min();
1.521 for (int i = 0; i != _arc_num; ++i) {
1.522 - if (_cost[i] > art_cost) art_cost = _cost[i];
1.523 + if (_cost[i] > ART_COST) ART_COST = _cost[i];
1.524 }
1.525 - art_cost = (art_cost + 1) * _node_num;
1.526 + ART_COST = (ART_COST + 1) * _node_num;
1.527 }
1.528
1.529 - // Run Circulation to check if a feasible solution exists
1.530 - typedef ConstMap<Arc, Flow> ConstArcMap;
1.531 - ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
1.532 - FlowNodeMap *csup = NULL;
1.533 - bool local_csup = false;
1.534 - if (_psupply) {
1.535 - csup = _psupply;
1.536 - } else {
1.537 - csup = new FlowNodeMap(_graph, 0);
1.538 - (*csup)[_psource] = _pstflow;
1.539 - (*csup)[_ptarget] = -_pstflow;
1.540 - local_csup = true;
1.541 - }
1.542 - bool circ_result = false;
1.543 - if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
1.544 - // GEQ problem type
1.545 - if (_plower) {
1.546 - if (_pupper) {
1.547 - Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
1.548 - circ(_graph, *_plower, *_pupper, *csup);
1.549 - circ_result = circ.run();
1.550 - } else {
1.551 - Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
1.552 - circ(_graph, *_plower, inf_arc_map, *csup);
1.553 - circ_result = circ.run();
1.554 - }
1.555 - } else {
1.556 - if (_pupper) {
1.557 - Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
1.558 - circ(_graph, zero_arc_map, *_pupper, *csup);
1.559 - circ_result = circ.run();
1.560 - } else {
1.561 - Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
1.562 - circ(_graph, zero_arc_map, inf_arc_map, *csup);
1.563 - circ_result = circ.run();
1.564 - }
1.565 - }
1.566 - } else {
1.567 - // LEQ problem type
1.568 - typedef ReverseDigraph<const GR> RevGraph;
1.569 - typedef NegMap<FlowNodeMap> NegNodeMap;
1.570 - RevGraph rgraph(_graph);
1.571 - NegNodeMap neg_csup(*csup);
1.572 - if (_plower) {
1.573 - if (_pupper) {
1.574 - Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
1.575 - circ(rgraph, *_plower, *_pupper, neg_csup);
1.576 - circ_result = circ.run();
1.577 - } else {
1.578 - Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
1.579 - circ(rgraph, *_plower, inf_arc_map, neg_csup);
1.580 - circ_result = circ.run();
1.581 - }
1.582 - } else {
1.583 - if (_pupper) {
1.584 - Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
1.585 - circ(rgraph, zero_arc_map, *_pupper, neg_csup);
1.586 - circ_result = circ.run();
1.587 - } else {
1.588 - Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
1.589 - circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
1.590 - circ_result = circ.run();
1.591 - }
1.592 - }
1.593 - }
1.594 - if (local_csup) delete csup;
1.595 - if (!circ_result) return false;
1.596 -
1.597 // Set data for the artificial root node
1.598 _root = _node_num;
1.599 _parent[_root] = -1;
1.600 @@ -1221,11 +1151,11 @@
1.601 _rev_thread[0] = _root;
1.602 _succ_num[_root] = all_node_num;
1.603 _last_succ[_root] = _root - 1;
1.604 - _supply[_root] = -sum_supply;
1.605 - if (sum_supply < 0) {
1.606 - _pi[_root] = -art_cost;
1.607 + _supply[_root] = -_sum_supply;
1.608 + if (_sum_supply < 0) {
1.609 + _pi[_root] = -ART_COST;
1.610 } else {
1.611 - _pi[_root] = art_cost;
1.612 + _pi[_root] = ART_COST;
1.613 }
1.614
1.615 // Store the arcs in a mixed order
1.616 @@ -1260,7 +1190,7 @@
1.617 _cap[i] = (*_pupper)[_arc_ref[i]];
1.618 } else {
1.619 for (int i = 0; i != _arc_num; ++i)
1.620 - _cap[i] = inf_cap;
1.621 + _cap[i] = INF;
1.622 }
1.623 if (_pcost) {
1.624 for (int i = 0; i != _arc_num; ++i)
1.625 @@ -1275,8 +1205,17 @@
1.626 if (_plower) {
1.627 for (int i = 0; i != _arc_num; ++i) {
1.628 Flow c = (*_plower)[_arc_ref[i]];
1.629 - if (c != 0) {
1.630 - _cap[i] -= c;
1.631 + if (c > 0) {
1.632 + if (_cap[i] < INF) _cap[i] -= c;
1.633 + _supply[_source[i]] -= c;
1.634 + _supply[_target[i]] += c;
1.635 + }
1.636 + else if (c < 0) {
1.637 + if (_cap[i] < INF + c) {
1.638 + _cap[i] -= c;
1.639 + } else {
1.640 + _cap[i] = INF;
1.641 + }
1.642 _supply[_source[i]] -= c;
1.643 _supply[_target[i]] += c;
1.644 }
1.645 @@ -1291,17 +1230,17 @@
1.646 _last_succ[u] = u;
1.647 _parent[u] = _root;
1.648 _pred[u] = e;
1.649 - _cost[e] = art_cost;
1.650 - _cap[e] = inf_cap;
1.651 + _cost[e] = ART_COST;
1.652 + _cap[e] = INF;
1.653 _state[e] = STATE_TREE;
1.654 - if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
1.655 + if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
1.656 _flow[e] = _supply[u];
1.657 _forward[u] = true;
1.658 - _pi[u] = -art_cost + _pi[_root];
1.659 + _pi[u] = -ART_COST + _pi[_root];
1.660 } else {
1.661 _flow[e] = -_supply[u];
1.662 _forward[u] = false;
1.663 - _pi[u] = art_cost + _pi[_root];
1.664 + _pi[u] = ART_COST + _pi[_root];
1.665 }
1.666 }
1.667
1.668 @@ -1342,7 +1281,8 @@
1.669 // Search the cycle along the path form the first node to the root
1.670 for (int u = first; u != join; u = _parent[u]) {
1.671 e = _pred[u];
1.672 - d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
1.673 + d = _forward[u] ?
1.674 + _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
1.675 if (d < delta) {
1.676 delta = d;
1.677 u_out = u;
1.678 @@ -1352,7 +1292,8 @@
1.679 // Search the cycle along the path form the second node to the root
1.680 for (int u = second; u != join; u = _parent[u]) {
1.681 e = _pred[u];
1.682 - d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1.683 + d = _forward[u] ?
1.684 + (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
1.685 if (d <= delta) {
1.686 delta = d;
1.687 u_out = u;
1.688 @@ -1526,7 +1467,7 @@
1.689 }
1.690
1.691 // Execute the algorithm
1.692 - bool start(PivotRule pivot_rule) {
1.693 + ProblemType start(PivotRule pivot_rule) {
1.694 // Select the pivot rule implementation
1.695 switch (pivot_rule) {
1.696 case FIRST_ELIGIBLE:
1.697 @@ -1540,23 +1481,41 @@
1.698 case ALTERING_LIST:
1.699 return start<AlteringListPivotRule>();
1.700 }
1.701 - return false;
1.702 + return INFEASIBLE; // avoid warning
1.703 }
1.704
1.705 template <typename PivotRuleImpl>
1.706 - bool start() {
1.707 + ProblemType start() {
1.708 PivotRuleImpl pivot(*this);
1.709
1.710 // Execute the Network Simplex algorithm
1.711 while (pivot.findEnteringArc()) {
1.712 findJoinNode();
1.713 bool change = findLeavingArc();
1.714 + if (delta >= INF) return UNBOUNDED;
1.715 changeFlow(change);
1.716 if (change) {
1.717 updateTreeStructure();
1.718 updatePotential();
1.719 }
1.720 }
1.721 +
1.722 + // Check feasibility
1.723 + if (_sum_supply < 0) {
1.724 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.725 + if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
1.726 + }
1.727 + }
1.728 + else if (_sum_supply > 0) {
1.729 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.730 + if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
1.731 + }
1.732 + }
1.733 + else {
1.734 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.735 + if (_flow[e] != 0) return INFEASIBLE;
1.736 + }
1.737 + }
1.738
1.739 // Copy flow values to _flow_map
1.740 if (_plower) {
1.741 @@ -1574,7 +1533,7 @@
1.742 _potential_map->set(n, _pi[_node_id[n]]);
1.743 }
1.744
1.745 - return true;
1.746 + return OPTIMAL;
1.747 }
1.748
1.749 }; //class NetworkSimplex