1.1 --- a/lemon/network_simplex.h Wed Apr 29 16:15:29 2009 +0100
1.2 +++ b/lemon/network_simplex.h Wed Apr 29 19:22:14 2009 +0100
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,14 +47,19 @@
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 + /// \tparam V The value type used for flow amounts, capacity bounds
1.30 /// and supply values in the algorithm. By default it is \c int.
1.31 /// \tparam C The value type used for costs and potentials in the
1.32 - /// algorithm. By default it is the same as \c F.
1.33 + /// algorithm. By default it is the same as \c V.
1.34 ///
1.35 /// \warning Both value types must be signed and all input data must
1.36 /// be integer.
1.37 @@ -65,34 +67,92 @@
1.38 /// \note %NetworkSimplex provides five different pivot rule
1.39 /// implementations, from which the most efficient one is used
1.40 /// by default. For more information see \ref PivotRule.
1.41 - template <typename GR, typename F = int, typename C = F>
1.42 + template <typename GR, typename V = int, typename C = V>
1.43 class NetworkSimplex
1.44 {
1.45 public:
1.46
1.47 - /// The flow type of the algorithm
1.48 - typedef F Flow;
1.49 - /// The cost type of the algorithm
1.50 + /// The type of the flow amounts, capacity bounds and supply values
1.51 + typedef V Value;
1.52 + /// The type of the arc costs
1.53 typedef C Cost;
1.54 -#ifdef DOXYGEN
1.55 - /// The type of the flow map
1.56 - typedef GR::ArcMap<Flow> FlowMap;
1.57 - /// The type of the potential map
1.58 - typedef GR::NodeMap<Cost> PotentialMap;
1.59 -#else
1.60 - /// The type of the flow map
1.61 - typedef typename GR::template ArcMap<Flow> FlowMap;
1.62 - /// The type of the potential map
1.63 - typedef typename GR::template NodeMap<Cost> PotentialMap;
1.64 -#endif
1.65
1.66 public:
1.67
1.68 - /// \brief Enum type for selecting the pivot rule.
1.69 + /// \brief Problem type constants for the \c run() function.
1.70 ///
1.71 - /// Enum type for selecting the pivot rule for the \ref run()
1.72 + /// Enum type containing the problem type constants that can be
1.73 + /// returned by the \ref run() function of the algorithm.
1.74 + enum ProblemType {
1.75 + /// The problem has no feasible solution (flow).
1.76 + INFEASIBLE,
1.77 + /// The problem has optimal solution (i.e. it is feasible and
1.78 + /// bounded), and the algorithm has found optimal flow and node
1.79 + /// potentials (primal and dual solutions).
1.80 + OPTIMAL,
1.81 + /// The objective function of the problem is unbounded, i.e.
1.82 + /// there is a directed cycle having negative total cost and
1.83 + /// infinite upper bound.
1.84 + UNBOUNDED
1.85 + };
1.86 +
1.87 + /// \brief Constants for selecting the type of the supply constraints.
1.88 + ///
1.89 + /// Enum type containing constants for selecting the supply type,
1.90 + /// i.e. the direction of the inequalities in the supply/demand
1.91 + /// constraints of the \ref min_cost_flow "minimum cost flow problem".
1.92 + ///
1.93 + /// The default supply type is \c GEQ, since this form is supported
1.94 + /// by other minimum cost flow algorithms and the \ref Circulation
1.95 + /// algorithm, as well.
1.96 + /// The \c LEQ problem type can be selected using the \ref supplyType()
1.97 /// function.
1.98 ///
1.99 + /// Note that the equality form is a special case of both supply types.
1.100 + enum SupplyType {
1.101 +
1.102 + /// This option means that there are <em>"greater or equal"</em>
1.103 + /// supply/demand constraints in the definition, i.e. the exact
1.104 + /// formulation of the problem is the following.
1.105 + /**
1.106 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.107 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
1.108 + sup(u) \quad \forall u\in V \f]
1.109 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.110 + */
1.111 + /// It means that the total demand must be greater or equal to the
1.112 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.113 + /// negative) and all the supplies have to be carried out from
1.114 + /// the supply nodes, but there could be demands that are not
1.115 + /// satisfied.
1.116 + GEQ,
1.117 + /// It is just an alias for the \c GEQ option.
1.118 + CARRY_SUPPLIES = GEQ,
1.119 +
1.120 + /// This option means that there are <em>"less or equal"</em>
1.121 + /// supply/demand constraints in the definition, i.e. the exact
1.122 + /// formulation of the problem is the following.
1.123 + /**
1.124 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.125 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
1.126 + sup(u) \quad \forall u\in V \f]
1.127 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.128 + */
1.129 + /// It means that the total demand must be less or equal to the
1.130 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.131 + /// positive) and all the demands have to be satisfied, but there
1.132 + /// could be supplies that are not carried out from the supply
1.133 + /// nodes.
1.134 + LEQ,
1.135 + /// It is just an alias for the \c LEQ option.
1.136 + SATISFY_DEMANDS = LEQ
1.137 + };
1.138 +
1.139 + /// \brief Constants for selecting the pivot rule.
1.140 + ///
1.141 + /// Enum type containing constants for selecting the pivot rule for
1.142 + /// the \ref run() function.
1.143 + ///
1.144 /// \ref NetworkSimplex provides five different pivot rule
1.145 /// implementations that significantly affect the running time
1.146 /// of the algorithm.
1.147 @@ -131,71 +191,15 @@
1.148 ALTERING_LIST
1.149 };
1.150
1.151 - /// \brief Enum type for selecting the problem type.
1.152 - ///
1.153 - /// Enum type for selecting the problem type, i.e. the direction of
1.154 - /// the inequalities in the supply/demand constraints of the
1.155 - /// \ref min_cost_flow "minimum cost flow problem".
1.156 - ///
1.157 - /// The default problem type is \c GEQ, since this form is supported
1.158 - /// by other minimum cost flow algorithms and the \ref Circulation
1.159 - /// algorithm as well.
1.160 - /// The \c LEQ problem type can be selected using the \ref problemType()
1.161 - /// function.
1.162 - ///
1.163 - /// Note that the equality form is a special case of both problem type.
1.164 - enum ProblemType {
1.165 -
1.166 - /// This option means that there are "<em>greater or equal</em>"
1.167 - /// constraints in the defintion, i.e. the exact formulation of the
1.168 - /// problem is the following.
1.169 - /**
1.170 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.171 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
1.172 - sup(u) \quad \forall u\in V \f]
1.173 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.174 - */
1.175 - /// It means that the total demand must be greater or equal to the
1.176 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.177 - /// negative) and all the supplies have to be carried out from
1.178 - /// the supply nodes, but there could be demands that are not
1.179 - /// satisfied.
1.180 - GEQ,
1.181 - /// It is just an alias for the \c GEQ option.
1.182 - CARRY_SUPPLIES = GEQ,
1.183 -
1.184 - /// This option means that there are "<em>less or equal</em>"
1.185 - /// constraints in the defintion, i.e. the exact formulation of the
1.186 - /// problem is the following.
1.187 - /**
1.188 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.189 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
1.190 - sup(u) \quad \forall u\in V \f]
1.191 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.192 - */
1.193 - /// It means that the total demand must be less or equal to the
1.194 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.195 - /// positive) and all the demands have to be satisfied, but there
1.196 - /// could be supplies that are not carried out from the supply
1.197 - /// nodes.
1.198 - LEQ,
1.199 - /// It is just an alias for the \c LEQ option.
1.200 - SATISFY_DEMANDS = LEQ
1.201 - };
1.202 -
1.203 private:
1.204
1.205 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.206
1.207 - typedef typename GR::template ArcMap<Flow> FlowArcMap;
1.208 - typedef typename GR::template ArcMap<Cost> CostArcMap;
1.209 - typedef typename GR::template NodeMap<Flow> FlowNodeMap;
1.210 -
1.211 typedef std::vector<Arc> ArcVector;
1.212 typedef std::vector<Node> NodeVector;
1.213 typedef std::vector<int> IntVector;
1.214 typedef std::vector<bool> BoolVector;
1.215 - typedef std::vector<Flow> FlowVector;
1.216 + typedef std::vector<Value> ValueVector;
1.217 typedef std::vector<Cost> CostVector;
1.218
1.219 // State constants for arcs
1.220 @@ -213,32 +217,23 @@
1.221 int _arc_num;
1.222
1.223 // Parameters of the problem
1.224 - FlowArcMap *_plower;
1.225 - FlowArcMap *_pupper;
1.226 - CostArcMap *_pcost;
1.227 - FlowNodeMap *_psupply;
1.228 - bool _pstsup;
1.229 - Node _psource, _ptarget;
1.230 - Flow _pstflow;
1.231 - ProblemType _ptype;
1.232 -
1.233 - // Result maps
1.234 - FlowMap *_flow_map;
1.235 - PotentialMap *_potential_map;
1.236 - bool _local_flow;
1.237 - bool _local_potential;
1.238 + bool _have_lower;
1.239 + SupplyType _stype;
1.240 + Value _sum_supply;
1.241
1.242 // Data structures for storing the digraph
1.243 IntNodeMap _node_id;
1.244 - ArcVector _arc_ref;
1.245 + IntArcMap _arc_id;
1.246 IntVector _source;
1.247 IntVector _target;
1.248
1.249 // Node and arc data
1.250 - FlowVector _cap;
1.251 + ValueVector _lower;
1.252 + ValueVector _upper;
1.253 + ValueVector _cap;
1.254 CostVector _cost;
1.255 - FlowVector _supply;
1.256 - FlowVector _flow;
1.257 + ValueVector _supply;
1.258 + ValueVector _flow;
1.259 CostVector _pi;
1.260
1.261 // Data for storing the spanning tree structure
1.262 @@ -257,7 +252,16 @@
1.263 int in_arc, join, u_in, v_in, u_out, v_out;
1.264 int first, second, right, last;
1.265 int stem, par_stem, new_stem;
1.266 - Flow delta;
1.267 + Value delta;
1.268 +
1.269 + public:
1.270 +
1.271 + /// \brief Constant for infinite upper bounds (capacities).
1.272 + ///
1.273 + /// Constant for infinite upper bounds (capacities).
1.274 + /// It is \c std::numeric_limits<Value>::infinity() if available,
1.275 + /// \c std::numeric_limits<Value>::max() otherwise.
1.276 + const Value INF;
1.277
1.278 private:
1.279
1.280 @@ -659,25 +663,68 @@
1.281 ///
1.282 /// \param graph The digraph the algorithm runs on.
1.283 NetworkSimplex(const GR& graph) :
1.284 - _graph(graph),
1.285 - _plower(NULL), _pupper(NULL), _pcost(NULL),
1.286 - _psupply(NULL), _pstsup(false), _ptype(GEQ),
1.287 - _flow_map(NULL), _potential_map(NULL),
1.288 - _local_flow(false), _local_potential(false),
1.289 - _node_id(graph)
1.290 + _graph(graph), _node_id(graph), _arc_id(graph),
1.291 + INF(std::numeric_limits<Value>::has_infinity ?
1.292 + std::numeric_limits<Value>::infinity() :
1.293 + std::numeric_limits<Value>::max())
1.294 {
1.295 - LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
1.296 - std::numeric_limits<Flow>::is_signed,
1.297 - "The flow type of NetworkSimplex must be signed integer");
1.298 - LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
1.299 - std::numeric_limits<Cost>::is_signed,
1.300 - "The cost type of NetworkSimplex must be signed integer");
1.301 - }
1.302 + // Check the value types
1.303 + LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
1.304 + "The flow type of NetworkSimplex must be signed");
1.305 + LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
1.306 + "The cost type of NetworkSimplex must be signed");
1.307 +
1.308 + // Resize vectors
1.309 + _node_num = countNodes(_graph);
1.310 + _arc_num = countArcs(_graph);
1.311 + int all_node_num = _node_num + 1;
1.312 + int all_arc_num = _arc_num + _node_num;
1.313
1.314 - /// Destructor.
1.315 - ~NetworkSimplex() {
1.316 - if (_local_flow) delete _flow_map;
1.317 - if (_local_potential) delete _potential_map;
1.318 + _source.resize(all_arc_num);
1.319 + _target.resize(all_arc_num);
1.320 +
1.321 + _lower.resize(all_arc_num);
1.322 + _upper.resize(all_arc_num);
1.323 + _cap.resize(all_arc_num);
1.324 + _cost.resize(all_arc_num);
1.325 + _supply.resize(all_node_num);
1.326 + _flow.resize(all_arc_num);
1.327 + _pi.resize(all_node_num);
1.328 +
1.329 + _parent.resize(all_node_num);
1.330 + _pred.resize(all_node_num);
1.331 + _forward.resize(all_node_num);
1.332 + _thread.resize(all_node_num);
1.333 + _rev_thread.resize(all_node_num);
1.334 + _succ_num.resize(all_node_num);
1.335 + _last_succ.resize(all_node_num);
1.336 + _state.resize(all_arc_num);
1.337 +
1.338 + // Copy the graph (store the arcs in a mixed order)
1.339 + int i = 0;
1.340 + for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.341 + _node_id[n] = i;
1.342 + }
1.343 + int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.344 + i = 0;
1.345 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.346 + _arc_id[a] = i;
1.347 + _source[i] = _node_id[_graph.source(a)];
1.348 + _target[i] = _node_id[_graph.target(a)];
1.349 + if ((i += k) >= _arc_num) i = (i % k) + 1;
1.350 + }
1.351 +
1.352 + // Initialize maps
1.353 + for (int i = 0; i != _node_num; ++i) {
1.354 + _supply[i] = 0;
1.355 + }
1.356 + for (int i = 0; i != _arc_num; ++i) {
1.357 + _lower[i] = 0;
1.358 + _upper[i] = INF;
1.359 + _cost[i] = 1;
1.360 + }
1.361 + _have_lower = false;
1.362 + _stype = GEQ;
1.363 }
1.364
1.365 /// \name Parameters
1.366 @@ -689,21 +736,19 @@
1.367 /// \brief Set the lower bounds on the arcs.
1.368 ///
1.369 /// This function sets the lower bounds on the arcs.
1.370 - /// If neither this function nor \ref boundMaps() is used before
1.371 - /// calling \ref run(), the lower bounds will be set to zero
1.372 - /// on all arcs.
1.373 + /// If it is not used before calling \ref run(), the lower bounds
1.374 + /// will be set to zero on all arcs.
1.375 ///
1.376 /// \param map An arc map storing the lower bounds.
1.377 - /// Its \c Value type must be convertible to the \c Flow type
1.378 + /// Its \c Value type must be convertible to the \c Value type
1.379 /// of the algorithm.
1.380 ///
1.381 /// \return <tt>(*this)</tt>
1.382 - template <typename LOWER>
1.383 - NetworkSimplex& lowerMap(const LOWER& map) {
1.384 - delete _plower;
1.385 - _plower = new FlowArcMap(_graph);
1.386 + template <typename LowerMap>
1.387 + NetworkSimplex& lowerMap(const LowerMap& map) {
1.388 + _have_lower = true;
1.389 for (ArcIt a(_graph); a != INVALID; ++a) {
1.390 - (*_plower)[a] = map[a];
1.391 + _lower[_arc_id[a]] = map[a];
1.392 }
1.393 return *this;
1.394 }
1.395 @@ -711,63 +756,23 @@
1.396 /// \brief Set the upper bounds (capacities) on the arcs.
1.397 ///
1.398 /// This function sets the upper bounds (capacities) on the arcs.
1.399 - /// If none of the functions \ref upperMap(), \ref capacityMap()
1.400 - /// and \ref boundMaps() is used before calling \ref run(),
1.401 - /// the upper bounds (capacities) will be set to
1.402 - /// \c std::numeric_limits<Flow>::max() on all arcs.
1.403 + /// If it is not used before calling \ref run(), the upper bounds
1.404 + /// will be set to \ref INF on all arcs (i.e. the flow value will be
1.405 + /// unbounded from above on each arc).
1.406 ///
1.407 /// \param map An arc map storing the upper bounds.
1.408 - /// Its \c Value type must be convertible to the \c Flow type
1.409 + /// Its \c Value type must be convertible to the \c Value type
1.410 /// of the algorithm.
1.411 ///
1.412 /// \return <tt>(*this)</tt>
1.413 - template<typename UPPER>
1.414 - NetworkSimplex& upperMap(const UPPER& map) {
1.415 - delete _pupper;
1.416 - _pupper = new FlowArcMap(_graph);
1.417 + template<typename UpperMap>
1.418 + NetworkSimplex& upperMap(const UpperMap& map) {
1.419 for (ArcIt a(_graph); a != INVALID; ++a) {
1.420 - (*_pupper)[a] = map[a];
1.421 + _upper[_arc_id[a]] = map[a];
1.422 }
1.423 return *this;
1.424 }
1.425
1.426 - /// \brief Set the upper bounds (capacities) on the arcs.
1.427 - ///
1.428 - /// This function sets the upper bounds (capacities) on the arcs.
1.429 - /// It is just an alias for \ref upperMap().
1.430 - ///
1.431 - /// \return <tt>(*this)</tt>
1.432 - template<typename CAP>
1.433 - NetworkSimplex& capacityMap(const CAP& map) {
1.434 - return upperMap(map);
1.435 - }
1.436 -
1.437 - /// \brief Set the lower and upper bounds on the arcs.
1.438 - ///
1.439 - /// This function sets the lower and upper bounds on the arcs.
1.440 - /// If neither this function nor \ref lowerMap() is used before
1.441 - /// calling \ref run(), the lower bounds will be set to zero
1.442 - /// on all arcs.
1.443 - /// If none of the functions \ref upperMap(), \ref capacityMap()
1.444 - /// and \ref boundMaps() is used before calling \ref run(),
1.445 - /// the upper bounds (capacities) will be set to
1.446 - /// \c std::numeric_limits<Flow>::max() on all arcs.
1.447 - ///
1.448 - /// \param lower An arc map storing the lower bounds.
1.449 - /// \param upper An arc map storing the upper bounds.
1.450 - ///
1.451 - /// The \c Value type of the maps must be convertible to the
1.452 - /// \c Flow type of the algorithm.
1.453 - ///
1.454 - /// \note This function is just a shortcut of calling \ref lowerMap()
1.455 - /// and \ref upperMap() separately.
1.456 - ///
1.457 - /// \return <tt>(*this)</tt>
1.458 - template <typename LOWER, typename UPPER>
1.459 - NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
1.460 - return lowerMap(lower).upperMap(upper);
1.461 - }
1.462 -
1.463 /// \brief Set the costs of the arcs.
1.464 ///
1.465 /// This function sets the costs of the arcs.
1.466 @@ -779,12 +784,10 @@
1.467 /// of the algorithm.
1.468 ///
1.469 /// \return <tt>(*this)</tt>
1.470 - template<typename COST>
1.471 - NetworkSimplex& costMap(const COST& map) {
1.472 - delete _pcost;
1.473 - _pcost = new CostArcMap(_graph);
1.474 + template<typename CostMap>
1.475 + NetworkSimplex& costMap(const CostMap& map) {
1.476 for (ArcIt a(_graph); a != INVALID; ++a) {
1.477 - (*_pcost)[a] = map[a];
1.478 + _cost[_arc_id[a]] = map[a];
1.479 }
1.480 return *this;
1.481 }
1.482 @@ -797,17 +800,14 @@
1.483 /// (It makes sense only if non-zero lower bounds are given.)
1.484 ///
1.485 /// \param map A node map storing the supply values.
1.486 - /// Its \c Value type must be convertible to the \c Flow type
1.487 + /// Its \c Value type must be convertible to the \c Value type
1.488 /// of the algorithm.
1.489 ///
1.490 /// \return <tt>(*this)</tt>
1.491 - template<typename SUP>
1.492 - NetworkSimplex& supplyMap(const SUP& map) {
1.493 - delete _psupply;
1.494 - _pstsup = false;
1.495 - _psupply = new FlowNodeMap(_graph);
1.496 + template<typename SupplyMap>
1.497 + NetworkSimplex& supplyMap(const SupplyMap& map) {
1.498 for (NodeIt n(_graph); n != INVALID; ++n) {
1.499 - (*_psupply)[n] = map[n];
1.500 + _supply[_node_id[n]] = map[n];
1.501 }
1.502 return *this;
1.503 }
1.504 @@ -820,71 +820,39 @@
1.505 /// calling \ref run(), the supply of each node will be set to zero.
1.506 /// (It makes sense only if non-zero lower bounds are given.)
1.507 ///
1.508 + /// Using this function has the same effect as using \ref supplyMap()
1.509 + /// with such a map in which \c k is assigned to \c s, \c -k is
1.510 + /// assigned to \c t and all other nodes have zero supply value.
1.511 + ///
1.512 /// \param s The source node.
1.513 /// \param t The target node.
1.514 /// \param k The required amount of flow from node \c s to node \c t
1.515 /// (i.e. the supply of \c s and the demand of \c t).
1.516 ///
1.517 /// \return <tt>(*this)</tt>
1.518 - NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
1.519 - delete _psupply;
1.520 - _psupply = NULL;
1.521 - _pstsup = true;
1.522 - _psource = s;
1.523 - _ptarget = t;
1.524 - _pstflow = k;
1.525 + NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
1.526 + for (int i = 0; i != _node_num; ++i) {
1.527 + _supply[i] = 0;
1.528 + }
1.529 + _supply[_node_id[s]] = k;
1.530 + _supply[_node_id[t]] = -k;
1.531 return *this;
1.532 }
1.533
1.534 - /// \brief Set the problem type.
1.535 + /// \brief Set the type of the supply constraints.
1.536 ///
1.537 - /// This function sets the problem type for the algorithm.
1.538 - /// If it is not used before calling \ref run(), the \ref GEQ problem
1.539 + /// This function sets the type of the supply/demand constraints.
1.540 + /// If it is not used before calling \ref run(), the \ref GEQ supply
1.541 /// type will be used.
1.542 ///
1.543 - /// For more information see \ref ProblemType.
1.544 + /// For more information see \ref SupplyType.
1.545 ///
1.546 /// \return <tt>(*this)</tt>
1.547 - NetworkSimplex& problemType(ProblemType problem_type) {
1.548 - _ptype = problem_type;
1.549 + NetworkSimplex& supplyType(SupplyType supply_type) {
1.550 + _stype = supply_type;
1.551 return *this;
1.552 }
1.553
1.554 - /// \brief Set the flow map.
1.555 - ///
1.556 - /// This function sets the flow map.
1.557 - /// If it is not used before calling \ref run(), an instance will
1.558 - /// be allocated automatically. The destructor deallocates this
1.559 - /// automatically allocated map, of course.
1.560 - ///
1.561 - /// \return <tt>(*this)</tt>
1.562 - NetworkSimplex& flowMap(FlowMap& map) {
1.563 - if (_local_flow) {
1.564 - delete _flow_map;
1.565 - _local_flow = false;
1.566 - }
1.567 - _flow_map = ↦
1.568 - return *this;
1.569 - }
1.570 -
1.571 - /// \brief Set the potential map.
1.572 - ///
1.573 - /// This function sets the potential map, which is used for storing
1.574 - /// the dual solution.
1.575 - /// If it is not used before calling \ref run(), an instance will
1.576 - /// be allocated automatically. The destructor deallocates this
1.577 - /// automatically allocated map, of course.
1.578 - ///
1.579 - /// \return <tt>(*this)</tt>
1.580 - NetworkSimplex& potentialMap(PotentialMap& map) {
1.581 - if (_local_potential) {
1.582 - delete _potential_map;
1.583 - _local_potential = false;
1.584 - }
1.585 - _potential_map = ↦
1.586 - return *this;
1.587 - }
1.588 -
1.589 /// @}
1.590
1.591 /// \name Execution Control
1.592 @@ -896,13 +864,12 @@
1.593 ///
1.594 /// This function runs the algorithm.
1.595 /// The paramters can be specified using functions \ref lowerMap(),
1.596 - /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
1.597 - /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.598 - /// \ref problemType(), \ref flowMap() and \ref potentialMap().
1.599 + /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.600 + /// \ref supplyType().
1.601 /// For example,
1.602 /// \code
1.603 /// NetworkSimplex<ListDigraph> ns(graph);
1.604 - /// ns.boundMaps(lower, upper).costMap(cost)
1.605 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
1.606 /// .supplyMap(sup).run();
1.607 /// \endcode
1.608 ///
1.609 @@ -910,33 +877,44 @@
1.610 /// that have been given are kept for the next call, unless
1.611 /// \ref reset() is called, thus only the modified parameters
1.612 /// have to be set again. See \ref reset() for examples.
1.613 + /// However the underlying digraph must not be modified after this
1.614 + /// class have been constructed, since it copies and extends the graph.
1.615 ///
1.616 /// \param pivot_rule The pivot rule that will be used during the
1.617 /// algorithm. For more information see \ref PivotRule.
1.618 ///
1.619 - /// \return \c true if a feasible flow can be found.
1.620 - bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
1.621 - return init() && start(pivot_rule);
1.622 + /// \return \c INFEASIBLE if no feasible flow exists,
1.623 + /// \n \c OPTIMAL if the problem has optimal solution
1.624 + /// (i.e. it is feasible and bounded), and the algorithm has found
1.625 + /// optimal flow and node potentials (primal and dual solutions),
1.626 + /// \n \c UNBOUNDED if the objective function of the problem is
1.627 + /// unbounded, i.e. there is a directed cycle having negative total
1.628 + /// cost and infinite upper bound.
1.629 + ///
1.630 + /// \see ProblemType, PivotRule
1.631 + ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
1.632 + if (!init()) return INFEASIBLE;
1.633 + return start(pivot_rule);
1.634 }
1.635
1.636 /// \brief Reset all the parameters that have been given before.
1.637 ///
1.638 /// This function resets all the paramaters that have been given
1.639 /// before using functions \ref lowerMap(), \ref upperMap(),
1.640 - /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
1.641 - /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
1.642 - /// \ref flowMap() and \ref potentialMap().
1.643 + /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
1.644 ///
1.645 /// It is useful for multiple run() calls. If this function is not
1.646 /// used, all the parameters given before are kept for the next
1.647 /// \ref run() call.
1.648 + /// However the underlying digraph must not be modified after this
1.649 + /// class have been constructed, since it copies and extends the graph.
1.650 ///
1.651 /// For example,
1.652 /// \code
1.653 /// NetworkSimplex<ListDigraph> ns(graph);
1.654 ///
1.655 /// // First run
1.656 - /// ns.lowerMap(lower).capacityMap(cap).costMap(cost)
1.657 + /// ns.lowerMap(lower).upperMap(upper).costMap(cost)
1.658 /// .supplyMap(sup).run();
1.659 ///
1.660 /// // Run again with modified cost map (reset() is not called,
1.661 @@ -947,29 +925,22 @@
1.662 /// // Run again from scratch using reset()
1.663 /// // (the lower bounds will be set to zero on all arcs)
1.664 /// ns.reset();
1.665 - /// ns.capacityMap(cap).costMap(cost)
1.666 + /// ns.upperMap(capacity).costMap(cost)
1.667 /// .supplyMap(sup).run();
1.668 /// \endcode
1.669 ///
1.670 /// \return <tt>(*this)</tt>
1.671 NetworkSimplex& reset() {
1.672 - delete _plower;
1.673 - delete _pupper;
1.674 - delete _pcost;
1.675 - delete _psupply;
1.676 - _plower = NULL;
1.677 - _pupper = NULL;
1.678 - _pcost = NULL;
1.679 - _psupply = NULL;
1.680 - _pstsup = false;
1.681 - _ptype = GEQ;
1.682 - if (_local_flow) delete _flow_map;
1.683 - if (_local_potential) delete _potential_map;
1.684 - _flow_map = NULL;
1.685 - _potential_map = NULL;
1.686 - _local_flow = false;
1.687 - _local_potential = false;
1.688 -
1.689 + for (int i = 0; i != _node_num; ++i) {
1.690 + _supply[i] = 0;
1.691 + }
1.692 + for (int i = 0; i != _arc_num; ++i) {
1.693 + _lower[i] = 0;
1.694 + _upper[i] = INF;
1.695 + _cost[i] = 1;
1.696 + }
1.697 + _have_lower = false;
1.698 + _stype = GEQ;
1.699 return *this;
1.700 }
1.701
1.702 @@ -985,7 +956,7 @@
1.703 /// \brief Return the total cost of the found flow.
1.704 ///
1.705 /// This function returns the total cost of the found flow.
1.706 - /// The complexity of the function is O(e).
1.707 + /// Its complexity is O(e).
1.708 ///
1.709 /// \note The return type of the function can be specified as a
1.710 /// template parameter. For example,
1.711 @@ -997,15 +968,12 @@
1.712 /// function.
1.713 ///
1.714 /// \pre \ref run() must be called before using this function.
1.715 - template <typename Num>
1.716 - Num totalCost() const {
1.717 - Num c = 0;
1.718 - if (_pcost) {
1.719 - for (ArcIt e(_graph); e != INVALID; ++e)
1.720 - c += (*_flow_map)[e] * (*_pcost)[e];
1.721 - } else {
1.722 - for (ArcIt e(_graph); e != INVALID; ++e)
1.723 - c += (*_flow_map)[e];
1.724 + template <typename Number>
1.725 + Number totalCost() const {
1.726 + Number c = 0;
1.727 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.728 + int i = _arc_id[a];
1.729 + c += Number(_flow[i]) * Number(_cost[i]);
1.730 }
1.731 return c;
1.732 }
1.733 @@ -1021,18 +989,22 @@
1.734 /// This function returns the flow on the given arc.
1.735 ///
1.736 /// \pre \ref run() must be called before using this function.
1.737 - Flow flow(const Arc& a) const {
1.738 - return (*_flow_map)[a];
1.739 + Value flow(const Arc& a) const {
1.740 + return _flow[_arc_id[a]];
1.741 }
1.742
1.743 - /// \brief Return a const reference to the flow map.
1.744 + /// \brief Return the flow map (the primal solution).
1.745 ///
1.746 - /// This function returns a const reference to an arc map storing
1.747 - /// the found flow.
1.748 + /// This function copies the flow value on each arc into the given
1.749 + /// map. The \c Value type of the algorithm must be convertible to
1.750 + /// the \c Value type of the map.
1.751 ///
1.752 /// \pre \ref run() must be called before using this function.
1.753 - const FlowMap& flowMap() const {
1.754 - return *_flow_map;
1.755 + template <typename FlowMap>
1.756 + void flowMap(FlowMap &map) const {
1.757 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.758 + map.set(a, _flow[_arc_id[a]]);
1.759 + }
1.760 }
1.761
1.762 /// \brief Return the potential (dual value) of the given node.
1.763 @@ -1042,19 +1014,22 @@
1.764 ///
1.765 /// \pre \ref run() must be called before using this function.
1.766 Cost potential(const Node& n) const {
1.767 - return (*_potential_map)[n];
1.768 + return _pi[_node_id[n]];
1.769 }
1.770
1.771 - /// \brief Return a const reference to the potential map
1.772 - /// (the dual solution).
1.773 + /// \brief Return the potential map (the dual solution).
1.774 ///
1.775 - /// This function returns a const reference to a node map storing
1.776 - /// the found potentials, which form the dual solution of the
1.777 - /// \ref min_cost_flow "minimum cost flow" problem.
1.778 + /// This function copies the potential (dual value) of each node
1.779 + /// into the given map.
1.780 + /// The \c Cost type of the algorithm must be convertible to the
1.781 + /// \c Value type of the map.
1.782 ///
1.783 /// \pre \ref run() must be called before using this function.
1.784 - const PotentialMap& potentialMap() const {
1.785 - return *_potential_map;
1.786 + template <typename PotentialMap>
1.787 + void potentialMap(PotentialMap &map) const {
1.788 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.789 + map.set(n, _pi[_node_id[n]]);
1.790 + }
1.791 }
1.792
1.793 /// @}
1.794 @@ -1063,245 +1038,82 @@
1.795
1.796 // Initialize internal data structures
1.797 bool init() {
1.798 - // Initialize result maps
1.799 - if (!_flow_map) {
1.800 - _flow_map = new FlowMap(_graph);
1.801 - _local_flow = true;
1.802 + if (_node_num == 0) return false;
1.803 +
1.804 + // Check the sum of supply values
1.805 + _sum_supply = 0;
1.806 + for (int i = 0; i != _node_num; ++i) {
1.807 + _sum_supply += _supply[i];
1.808 }
1.809 - if (!_potential_map) {
1.810 - _potential_map = new PotentialMap(_graph);
1.811 - _local_potential = true;
1.812 + if ( !((_stype == GEQ && _sum_supply <= 0) ||
1.813 + (_stype == LEQ && _sum_supply >= 0)) ) return false;
1.814 +
1.815 + // Remove non-zero lower bounds
1.816 + if (_have_lower) {
1.817 + for (int i = 0; i != _arc_num; ++i) {
1.818 + Value c = _lower[i];
1.819 + if (c >= 0) {
1.820 + _cap[i] = _upper[i] < INF ? _upper[i] - c : INF;
1.821 + } else {
1.822 + _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF;
1.823 + }
1.824 + _supply[_source[i]] -= c;
1.825 + _supply[_target[i]] += c;
1.826 + }
1.827 + } else {
1.828 + for (int i = 0; i != _arc_num; ++i) {
1.829 + _cap[i] = _upper[i];
1.830 + }
1.831 }
1.832
1.833 - // Initialize vectors
1.834 - _node_num = countNodes(_graph);
1.835 - _arc_num = countArcs(_graph);
1.836 - int all_node_num = _node_num + 1;
1.837 - int all_arc_num = _arc_num + _node_num;
1.838 - if (_node_num == 0) return false;
1.839 -
1.840 - _arc_ref.resize(_arc_num);
1.841 - _source.resize(all_arc_num);
1.842 - _target.resize(all_arc_num);
1.843 -
1.844 - _cap.resize(all_arc_num);
1.845 - _cost.resize(all_arc_num);
1.846 - _supply.resize(all_node_num);
1.847 - _flow.resize(all_arc_num);
1.848 - _pi.resize(all_node_num);
1.849 -
1.850 - _parent.resize(all_node_num);
1.851 - _pred.resize(all_node_num);
1.852 - _forward.resize(all_node_num);
1.853 - _thread.resize(all_node_num);
1.854 - _rev_thread.resize(all_node_num);
1.855 - _succ_num.resize(all_node_num);
1.856 - _last_succ.resize(all_node_num);
1.857 - _state.resize(all_arc_num);
1.858 -
1.859 - // Initialize node related data
1.860 - bool valid_supply = true;
1.861 - Flow sum_supply = 0;
1.862 - if (!_pstsup && !_psupply) {
1.863 - _pstsup = true;
1.864 - _psource = _ptarget = NodeIt(_graph);
1.865 - _pstflow = 0;
1.866 - }
1.867 - if (_psupply) {
1.868 - int i = 0;
1.869 - for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.870 - _node_id[n] = i;
1.871 - _supply[i] = (*_psupply)[n];
1.872 - sum_supply += _supply[i];
1.873 + // Initialize artifical cost
1.874 + Cost ART_COST;
1.875 + if (std::numeric_limits<Cost>::is_exact) {
1.876 + ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
1.877 + } else {
1.878 + ART_COST = std::numeric_limits<Cost>::min();
1.879 + for (int i = 0; i != _arc_num; ++i) {
1.880 + if (_cost[i] > ART_COST) ART_COST = _cost[i];
1.881 }
1.882 - valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
1.883 - (_ptype == LEQ && sum_supply >= 0);
1.884 - } else {
1.885 - int i = 0;
1.886 - for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.887 - _node_id[n] = i;
1.888 - _supply[i] = 0;
1.889 - }
1.890 - _supply[_node_id[_psource]] = _pstflow;
1.891 - _supply[_node_id[_ptarget]] = -_pstflow;
1.892 - }
1.893 - if (!valid_supply) return false;
1.894 -
1.895 - // Infinite capacity value
1.896 - Flow inf_cap =
1.897 - std::numeric_limits<Flow>::has_infinity ?
1.898 - std::numeric_limits<Flow>::infinity() :
1.899 - std::numeric_limits<Flow>::max();
1.900 -
1.901 - // Initialize artifical cost
1.902 - Cost art_cost;
1.903 - if (std::numeric_limits<Cost>::is_exact) {
1.904 - art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1.905 - } else {
1.906 - art_cost = std::numeric_limits<Cost>::min();
1.907 - for (int i = 0; i != _arc_num; ++i) {
1.908 - if (_cost[i] > art_cost) art_cost = _cost[i];
1.909 - }
1.910 - art_cost = (art_cost + 1) * _node_num;
1.911 + ART_COST = (ART_COST + 1) * _node_num;
1.912 }
1.913
1.914 - // Run Circulation to check if a feasible solution exists
1.915 - typedef ConstMap<Arc, Flow> ConstArcMap;
1.916 - ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
1.917 - FlowNodeMap *csup = NULL;
1.918 - bool local_csup = false;
1.919 - if (_psupply) {
1.920 - csup = _psupply;
1.921 - } else {
1.922 - csup = new FlowNodeMap(_graph, 0);
1.923 - (*csup)[_psource] = _pstflow;
1.924 - (*csup)[_ptarget] = -_pstflow;
1.925 - local_csup = true;
1.926 + // Initialize arc maps
1.927 + for (int i = 0; i != _arc_num; ++i) {
1.928 + _flow[i] = 0;
1.929 + _state[i] = STATE_LOWER;
1.930 }
1.931 - bool circ_result = false;
1.932 - if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
1.933 - // GEQ problem type
1.934 - if (_plower) {
1.935 - if (_pupper) {
1.936 - Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
1.937 - circ(_graph, *_plower, *_pupper, *csup);
1.938 - circ_result = circ.run();
1.939 - } else {
1.940 - Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
1.941 - circ(_graph, *_plower, inf_arc_map, *csup);
1.942 - circ_result = circ.run();
1.943 - }
1.944 - } else {
1.945 - if (_pupper) {
1.946 - Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
1.947 - circ(_graph, zero_arc_map, *_pupper, *csup);
1.948 - circ_result = circ.run();
1.949 - } else {
1.950 - Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
1.951 - circ(_graph, zero_arc_map, inf_arc_map, *csup);
1.952 - circ_result = circ.run();
1.953 - }
1.954 - }
1.955 - } else {
1.956 - // LEQ problem type
1.957 - typedef ReverseDigraph<const GR> RevGraph;
1.958 - typedef NegMap<FlowNodeMap> NegNodeMap;
1.959 - RevGraph rgraph(_graph);
1.960 - NegNodeMap neg_csup(*csup);
1.961 - if (_plower) {
1.962 - if (_pupper) {
1.963 - Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
1.964 - circ(rgraph, *_plower, *_pupper, neg_csup);
1.965 - circ_result = circ.run();
1.966 - } else {
1.967 - Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
1.968 - circ(rgraph, *_plower, inf_arc_map, neg_csup);
1.969 - circ_result = circ.run();
1.970 - }
1.971 - } else {
1.972 - if (_pupper) {
1.973 - Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
1.974 - circ(rgraph, zero_arc_map, *_pupper, neg_csup);
1.975 - circ_result = circ.run();
1.976 - } else {
1.977 - Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
1.978 - circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
1.979 - circ_result = circ.run();
1.980 - }
1.981 - }
1.982 - }
1.983 - if (local_csup) delete csup;
1.984 - if (!circ_result) return false;
1.985 -
1.986 +
1.987 // Set data for the artificial root node
1.988 _root = _node_num;
1.989 _parent[_root] = -1;
1.990 _pred[_root] = -1;
1.991 _thread[_root] = 0;
1.992 _rev_thread[0] = _root;
1.993 - _succ_num[_root] = all_node_num;
1.994 + _succ_num[_root] = _node_num + 1;
1.995 _last_succ[_root] = _root - 1;
1.996 - _supply[_root] = -sum_supply;
1.997 - if (sum_supply < 0) {
1.998 - _pi[_root] = -art_cost;
1.999 - } else {
1.1000 - _pi[_root] = art_cost;
1.1001 - }
1.1002 -
1.1003 - // Store the arcs in a mixed order
1.1004 - int k = std::max(int(std::sqrt(double(_arc_num))), 10);
1.1005 - int i = 0;
1.1006 - for (ArcIt e(_graph); e != INVALID; ++e) {
1.1007 - _arc_ref[i] = e;
1.1008 - if ((i += k) >= _arc_num) i = (i % k) + 1;
1.1009 - }
1.1010 -
1.1011 - // Initialize arc maps
1.1012 - if (_pupper && _pcost) {
1.1013 - for (int i = 0; i != _arc_num; ++i) {
1.1014 - Arc e = _arc_ref[i];
1.1015 - _source[i] = _node_id[_graph.source(e)];
1.1016 - _target[i] = _node_id[_graph.target(e)];
1.1017 - _cap[i] = (*_pupper)[e];
1.1018 - _cost[i] = (*_pcost)[e];
1.1019 - _flow[i] = 0;
1.1020 - _state[i] = STATE_LOWER;
1.1021 - }
1.1022 - } else {
1.1023 - for (int i = 0; i != _arc_num; ++i) {
1.1024 - Arc e = _arc_ref[i];
1.1025 - _source[i] = _node_id[_graph.source(e)];
1.1026 - _target[i] = _node_id[_graph.target(e)];
1.1027 - _flow[i] = 0;
1.1028 - _state[i] = STATE_LOWER;
1.1029 - }
1.1030 - if (_pupper) {
1.1031 - for (int i = 0; i != _arc_num; ++i)
1.1032 - _cap[i] = (*_pupper)[_arc_ref[i]];
1.1033 - } else {
1.1034 - for (int i = 0; i != _arc_num; ++i)
1.1035 - _cap[i] = inf_cap;
1.1036 - }
1.1037 - if (_pcost) {
1.1038 - for (int i = 0; i != _arc_num; ++i)
1.1039 - _cost[i] = (*_pcost)[_arc_ref[i]];
1.1040 - } else {
1.1041 - for (int i = 0; i != _arc_num; ++i)
1.1042 - _cost[i] = 1;
1.1043 - }
1.1044 - }
1.1045 -
1.1046 - // Remove non-zero lower bounds
1.1047 - if (_plower) {
1.1048 - for (int i = 0; i != _arc_num; ++i) {
1.1049 - Flow c = (*_plower)[_arc_ref[i]];
1.1050 - if (c != 0) {
1.1051 - _cap[i] -= c;
1.1052 - _supply[_source[i]] -= c;
1.1053 - _supply[_target[i]] += c;
1.1054 - }
1.1055 - }
1.1056 - }
1.1057 + _supply[_root] = -_sum_supply;
1.1058 + _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
1.1059
1.1060 // Add artificial arcs and initialize the spanning tree data structure
1.1061 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.1062 + _parent[u] = _root;
1.1063 + _pred[u] = e;
1.1064 _thread[u] = u + 1;
1.1065 _rev_thread[u + 1] = u;
1.1066 _succ_num[u] = 1;
1.1067 _last_succ[u] = u;
1.1068 - _parent[u] = _root;
1.1069 - _pred[u] = e;
1.1070 - _cost[e] = art_cost;
1.1071 - _cap[e] = inf_cap;
1.1072 + _cost[e] = ART_COST;
1.1073 + _cap[e] = INF;
1.1074 _state[e] = STATE_TREE;
1.1075 - if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
1.1076 + if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
1.1077 _flow[e] = _supply[u];
1.1078 _forward[u] = true;
1.1079 - _pi[u] = -art_cost + _pi[_root];
1.1080 + _pi[u] = -ART_COST + _pi[_root];
1.1081 } else {
1.1082 _flow[e] = -_supply[u];
1.1083 _forward[u] = false;
1.1084 - _pi[u] = art_cost + _pi[_root];
1.1085 + _pi[u] = ART_COST + _pi[_root];
1.1086 }
1.1087 }
1.1088
1.1089 @@ -1336,13 +1148,14 @@
1.1090 }
1.1091 delta = _cap[in_arc];
1.1092 int result = 0;
1.1093 - Flow d;
1.1094 + Value d;
1.1095 int e;
1.1096
1.1097 // Search the cycle along the path form the first node to the root
1.1098 for (int u = first; u != join; u = _parent[u]) {
1.1099 e = _pred[u];
1.1100 - d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
1.1101 + d = _forward[u] ?
1.1102 + _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
1.1103 if (d < delta) {
1.1104 delta = d;
1.1105 u_out = u;
1.1106 @@ -1352,7 +1165,8 @@
1.1107 // Search the cycle along the path form the second node to the root
1.1108 for (int u = second; u != join; u = _parent[u]) {
1.1109 e = _pred[u];
1.1110 - d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1.1111 + d = _forward[u] ?
1.1112 + (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
1.1113 if (d <= delta) {
1.1114 delta = d;
1.1115 u_out = u;
1.1116 @@ -1374,7 +1188,7 @@
1.1117 void changeFlow(bool change) {
1.1118 // Augment along the cycle
1.1119 if (delta > 0) {
1.1120 - Flow val = _state[in_arc] * delta;
1.1121 + Value val = _state[in_arc] * delta;
1.1122 _flow[in_arc] += val;
1.1123 for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1.1124 _flow[_pred[u]] += _forward[u] ? -val : val;
1.1125 @@ -1526,7 +1340,7 @@
1.1126 }
1.1127
1.1128 // Execute the algorithm
1.1129 - bool start(PivotRule pivot_rule) {
1.1130 + ProblemType start(PivotRule pivot_rule) {
1.1131 // Select the pivot rule implementation
1.1132 switch (pivot_rule) {
1.1133 case FIRST_ELIGIBLE:
1.1134 @@ -1540,41 +1354,55 @@
1.1135 case ALTERING_LIST:
1.1136 return start<AlteringListPivotRule>();
1.1137 }
1.1138 - return false;
1.1139 + return INFEASIBLE; // avoid warning
1.1140 }
1.1141
1.1142 template <typename PivotRuleImpl>
1.1143 - bool start() {
1.1144 + ProblemType start() {
1.1145 PivotRuleImpl pivot(*this);
1.1146
1.1147 // Execute the Network Simplex algorithm
1.1148 while (pivot.findEnteringArc()) {
1.1149 findJoinNode();
1.1150 bool change = findLeavingArc();
1.1151 + if (delta >= INF) return UNBOUNDED;
1.1152 changeFlow(change);
1.1153 if (change) {
1.1154 updateTreeStructure();
1.1155 updatePotential();
1.1156 }
1.1157 }
1.1158 -
1.1159 - // Copy flow values to _flow_map
1.1160 - if (_plower) {
1.1161 - for (int i = 0; i != _arc_num; ++i) {
1.1162 - Arc e = _arc_ref[i];
1.1163 - _flow_map->set(e, (*_plower)[e] + _flow[i]);
1.1164 - }
1.1165 - } else {
1.1166 - for (int i = 0; i != _arc_num; ++i) {
1.1167 - _flow_map->set(_arc_ref[i], _flow[i]);
1.1168 +
1.1169 + // Check feasibility
1.1170 + if (_sum_supply < 0) {
1.1171 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.1172 + if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
1.1173 }
1.1174 }
1.1175 - // Copy potential values to _potential_map
1.1176 - for (NodeIt n(_graph); n != INVALID; ++n) {
1.1177 - _potential_map->set(n, _pi[_node_id[n]]);
1.1178 + else if (_sum_supply > 0) {
1.1179 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.1180 + if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
1.1181 + }
1.1182 + }
1.1183 + else {
1.1184 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.1185 + if (_flow[e] != 0) return INFEASIBLE;
1.1186 + }
1.1187 }
1.1188
1.1189 - return true;
1.1190 + // Transform the solution and the supply map to the original form
1.1191 + if (_have_lower) {
1.1192 + for (int i = 0; i != _arc_num; ++i) {
1.1193 + Value c = _lower[i];
1.1194 + if (c != 0) {
1.1195 + _flow[i] += c;
1.1196 + _supply[_source[i]] += c;
1.1197 + _supply[_target[i]] -= c;
1.1198 + }
1.1199 + }
1.1200 + }
1.1201 +
1.1202 + return OPTIMAL;
1.1203 }
1.1204
1.1205 }; //class NetworkSimplex