1.1 --- a/lemon/network_simplex.h Tue Mar 24 00:18:25 2009 +0100
1.2 +++ b/lemon/network_simplex.h Wed Mar 25 15:58:44 2009 +0100
1.3 @@ -22,7 +22,7 @@
1.4 /// \ingroup min_cost_flow
1.5 ///
1.6 /// \file
1.7 -/// \brief Network simplex algorithm for finding a minimum cost flow.
1.8 +/// \brief Network Simplex algorithm for finding a minimum cost flow.
1.9
1.10 #include <vector>
1.11 #include <limits>
1.12 @@ -36,80 +36,89 @@
1.13 /// \addtogroup min_cost_flow
1.14 /// @{
1.15
1.16 - /// \brief Implementation of the primal network simplex algorithm
1.17 + /// \brief Implementation of the primal Network Simplex algorithm
1.18 /// for finding a \ref min_cost_flow "minimum cost flow".
1.19 ///
1.20 - /// \ref NetworkSimplex implements the primal network simplex algorithm
1.21 + /// \ref NetworkSimplex implements the primal Network Simplex algorithm
1.22 /// for finding a \ref min_cost_flow "minimum cost flow".
1.23 ///
1.24 - /// \tparam Digraph The digraph type the algorithm runs on.
1.25 - /// \tparam LowerMap The type of the lower bound map.
1.26 - /// \tparam CapacityMap The type of the capacity (upper bound) map.
1.27 - /// \tparam CostMap The type of the cost (length) map.
1.28 - /// \tparam SupplyMap The type of the supply map.
1.29 + /// \tparam GR The digraph type the algorithm runs on.
1.30 + /// \tparam V The value type used in the algorithm.
1.31 + /// By default it is \c int.
1.32 ///
1.33 - /// \warning
1.34 - /// - Arc capacities and costs should be \e non-negative \e integers.
1.35 - /// - Supply values should be \e signed \e integers.
1.36 - /// - The value types of the maps should be convertible to each other.
1.37 - /// - \c CostMap::Value must be signed type.
1.38 + /// \warning \c V must be a signed integer type.
1.39 ///
1.40 - /// \note \ref NetworkSimplex provides five different pivot rule
1.41 - /// implementations that significantly affect the efficiency of the
1.42 - /// algorithm.
1.43 - /// By default "Block Search" pivot rule is used, which proved to be
1.44 - /// by far the most efficient according to our benchmark tests.
1.45 - /// However another pivot rule can be selected using \ref run()
1.46 - /// function with the proper parameter.
1.47 -#ifdef DOXYGEN
1.48 - template < typename Digraph,
1.49 - typename LowerMap,
1.50 - typename CapacityMap,
1.51 - typename CostMap,
1.52 - typename SupplyMap >
1.53 -
1.54 -#else
1.55 - template < typename Digraph,
1.56 - typename LowerMap = typename Digraph::template ArcMap<int>,
1.57 - typename CapacityMap = typename Digraph::template ArcMap<int>,
1.58 - typename CostMap = typename Digraph::template ArcMap<int>,
1.59 - typename SupplyMap = typename Digraph::template NodeMap<int> >
1.60 -#endif
1.61 + /// \note %NetworkSimplex provides five different pivot rule
1.62 + /// implementations. For more information see \ref PivotRule.
1.63 + template <typename GR, typename V = int>
1.64 class NetworkSimplex
1.65 {
1.66 - TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.67 + public:
1.68
1.69 - typedef typename CapacityMap::Value Capacity;
1.70 - typedef typename CostMap::Value Cost;
1.71 - typedef typename SupplyMap::Value Supply;
1.72 + /// The value type of the algorithm
1.73 + typedef V Value;
1.74 + /// The type of the flow map
1.75 + typedef typename GR::template ArcMap<Value> FlowMap;
1.76 + /// The type of the potential map
1.77 + typedef typename GR::template NodeMap<Value> PotentialMap;
1.78 +
1.79 + public:
1.80 +
1.81 + /// \brief Enum type for selecting the pivot rule.
1.82 + ///
1.83 + /// Enum type for selecting the pivot rule for the \ref run()
1.84 + /// function.
1.85 + ///
1.86 + /// \ref NetworkSimplex provides five different pivot rule
1.87 + /// implementations that significantly affect the running time
1.88 + /// of the algorithm.
1.89 + /// By default \ref BLOCK_SEARCH "Block Search" is used, which
1.90 + /// proved to be the most efficient and the most robust on various
1.91 + /// test inputs according to our benchmark tests.
1.92 + /// However another pivot rule can be selected using the \ref run()
1.93 + /// function with the proper parameter.
1.94 + enum PivotRule {
1.95 +
1.96 + /// The First Eligible pivot rule.
1.97 + /// The next eligible arc is selected in a wraparound fashion
1.98 + /// in every iteration.
1.99 + FIRST_ELIGIBLE,
1.100 +
1.101 + /// The Best Eligible pivot rule.
1.102 + /// The best eligible arc is selected in every iteration.
1.103 + BEST_ELIGIBLE,
1.104 +
1.105 + /// The Block Search pivot rule.
1.106 + /// A specified number of arcs are examined in every iteration
1.107 + /// in a wraparound fashion and the best eligible arc is selected
1.108 + /// from this block.
1.109 + BLOCK_SEARCH,
1.110 +
1.111 + /// The Candidate List pivot rule.
1.112 + /// In a major iteration a candidate list is built from eligible arcs
1.113 + /// in a wraparound fashion and in the following minor iterations
1.114 + /// the best eligible arc is selected from this list.
1.115 + CANDIDATE_LIST,
1.116 +
1.117 + /// The Altering Candidate List pivot rule.
1.118 + /// It is a modified version of the Candidate List method.
1.119 + /// It keeps only the several best eligible arcs from the former
1.120 + /// candidate list and extends this list in every iteration.
1.121 + ALTERING_LIST
1.122 + };
1.123 +
1.124 + private:
1.125 +
1.126 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.127 +
1.128 + typedef typename GR::template ArcMap<Value> ValueArcMap;
1.129 + typedef typename GR::template NodeMap<Value> ValueNodeMap;
1.130
1.131 typedef std::vector<Arc> ArcVector;
1.132 typedef std::vector<Node> NodeVector;
1.133 typedef std::vector<int> IntVector;
1.134 typedef std::vector<bool> BoolVector;
1.135 - typedef std::vector<Capacity> CapacityVector;
1.136 - typedef std::vector<Cost> CostVector;
1.137 - typedef std::vector<Supply> SupplyVector;
1.138 -
1.139 - public:
1.140 -
1.141 - /// The type of the flow map
1.142 - typedef typename Digraph::template ArcMap<Capacity> FlowMap;
1.143 - /// The type of the potential map
1.144 - typedef typename Digraph::template NodeMap<Cost> PotentialMap;
1.145 -
1.146 - public:
1.147 -
1.148 - /// Enum type for selecting the pivot rule used by \ref run()
1.149 - enum PivotRuleEnum {
1.150 - FIRST_ELIGIBLE_PIVOT,
1.151 - BEST_ELIGIBLE_PIVOT,
1.152 - BLOCK_SEARCH_PIVOT,
1.153 - CANDIDATE_LIST_PIVOT,
1.154 - ALTERING_LIST_PIVOT
1.155 - };
1.156 -
1.157 - private:
1.158 + typedef std::vector<Value> ValueVector;
1.159
1.160 // State constants for arcs
1.161 enum ArcStateEnum {
1.162 @@ -120,15 +129,19 @@
1.163
1.164 private:
1.165
1.166 - // References for the original data
1.167 - const Digraph &_graph;
1.168 - const LowerMap *_orig_lower;
1.169 - const CapacityMap &_orig_cap;
1.170 - const CostMap &_orig_cost;
1.171 - const SupplyMap *_orig_supply;
1.172 - Node _orig_source;
1.173 - Node _orig_target;
1.174 - Capacity _orig_flow_value;
1.175 + // Data related to the underlying digraph
1.176 + const GR &_graph;
1.177 + int _node_num;
1.178 + int _arc_num;
1.179 +
1.180 + // Parameters of the problem
1.181 + ValueArcMap *_plower;
1.182 + ValueArcMap *_pupper;
1.183 + ValueArcMap *_pcost;
1.184 + ValueNodeMap *_psupply;
1.185 + bool _pstsup;
1.186 + Node _psource, _ptarget;
1.187 + Value _pstflow;
1.188
1.189 // Result maps
1.190 FlowMap *_flow_map;
1.191 @@ -136,22 +149,18 @@
1.192 bool _local_flow;
1.193 bool _local_potential;
1.194
1.195 - // The number of nodes and arcs in the original graph
1.196 - int _node_num;
1.197 - int _arc_num;
1.198 -
1.199 - // Data structures for storing the graph
1.200 + // Data structures for storing the digraph
1.201 IntNodeMap _node_id;
1.202 ArcVector _arc_ref;
1.203 IntVector _source;
1.204 IntVector _target;
1.205
1.206 - // Node and arc maps
1.207 - CapacityVector _cap;
1.208 - CostVector _cost;
1.209 - CostVector _supply;
1.210 - CapacityVector _flow;
1.211 - CostVector _pi;
1.212 + // Node and arc data
1.213 + ValueVector _cap;
1.214 + ValueVector _cost;
1.215 + ValueVector _supply;
1.216 + ValueVector _flow;
1.217 + ValueVector _pi;
1.218
1.219 // Data for storing the spanning tree structure
1.220 IntVector _parent;
1.221 @@ -169,17 +178,11 @@
1.222 int in_arc, join, u_in, v_in, u_out, v_out;
1.223 int first, second, right, last;
1.224 int stem, par_stem, new_stem;
1.225 - Capacity delta;
1.226 + Value delta;
1.227
1.228 private:
1.229
1.230 - /// \brief Implementation of the "First Eligible" pivot rule for the
1.231 - /// \ref NetworkSimplex "network simplex" algorithm.
1.232 - ///
1.233 - /// This class implements the "First Eligible" pivot rule
1.234 - /// for the \ref NetworkSimplex "network simplex" algorithm.
1.235 - ///
1.236 - /// For more information see \ref NetworkSimplex::run().
1.237 + // Implementation of the First Eligible pivot rule
1.238 class FirstEligiblePivotRule
1.239 {
1.240 private:
1.241 @@ -187,9 +190,9 @@
1.242 // References to the NetworkSimplex class
1.243 const IntVector &_source;
1.244 const IntVector &_target;
1.245 - const CostVector &_cost;
1.246 + const ValueVector &_cost;
1.247 const IntVector &_state;
1.248 - const CostVector &_pi;
1.249 + const ValueVector &_pi;
1.250 int &_in_arc;
1.251 int _arc_num;
1.252
1.253 @@ -198,16 +201,16 @@
1.254
1.255 public:
1.256
1.257 - /// Constructor
1.258 + // Constructor
1.259 FirstEligiblePivotRule(NetworkSimplex &ns) :
1.260 _source(ns._source), _target(ns._target),
1.261 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.262 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
1.263 {}
1.264
1.265 - /// Find next entering arc
1.266 + // Find next entering arc
1.267 bool findEnteringArc() {
1.268 - Cost c;
1.269 + Value c;
1.270 for (int e = _next_arc; e < _arc_num; ++e) {
1.271 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.272 if (c < 0) {
1.273 @@ -230,13 +233,7 @@
1.274 }; //class FirstEligiblePivotRule
1.275
1.276
1.277 - /// \brief Implementation of the "Best Eligible" pivot rule for the
1.278 - /// \ref NetworkSimplex "network simplex" algorithm.
1.279 - ///
1.280 - /// This class implements the "Best Eligible" pivot rule
1.281 - /// for the \ref NetworkSimplex "network simplex" algorithm.
1.282 - ///
1.283 - /// For more information see \ref NetworkSimplex::run().
1.284 + // Implementation of the Best Eligible pivot rule
1.285 class BestEligiblePivotRule
1.286 {
1.287 private:
1.288 @@ -244,24 +241,24 @@
1.289 // References to the NetworkSimplex class
1.290 const IntVector &_source;
1.291 const IntVector &_target;
1.292 - const CostVector &_cost;
1.293 + const ValueVector &_cost;
1.294 const IntVector &_state;
1.295 - const CostVector &_pi;
1.296 + const ValueVector &_pi;
1.297 int &_in_arc;
1.298 int _arc_num;
1.299
1.300 public:
1.301
1.302 - /// Constructor
1.303 + // Constructor
1.304 BestEligiblePivotRule(NetworkSimplex &ns) :
1.305 _source(ns._source), _target(ns._target),
1.306 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.307 _in_arc(ns.in_arc), _arc_num(ns._arc_num)
1.308 {}
1.309
1.310 - /// Find next entering arc
1.311 + // Find next entering arc
1.312 bool findEnteringArc() {
1.313 - Cost c, min = 0;
1.314 + Value c, min = 0;
1.315 for (int e = 0; e < _arc_num; ++e) {
1.316 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.317 if (c < min) {
1.318 @@ -275,13 +272,7 @@
1.319 }; //class BestEligiblePivotRule
1.320
1.321
1.322 - /// \brief Implementation of the "Block Search" pivot rule for the
1.323 - /// \ref NetworkSimplex "network simplex" algorithm.
1.324 - ///
1.325 - /// This class implements the "Block Search" pivot rule
1.326 - /// for the \ref NetworkSimplex "network simplex" algorithm.
1.327 - ///
1.328 - /// For more information see \ref NetworkSimplex::run().
1.329 + // Implementation of the Block Search pivot rule
1.330 class BlockSearchPivotRule
1.331 {
1.332 private:
1.333 @@ -289,9 +280,9 @@
1.334 // References to the NetworkSimplex class
1.335 const IntVector &_source;
1.336 const IntVector &_target;
1.337 - const CostVector &_cost;
1.338 + const ValueVector &_cost;
1.339 const IntVector &_state;
1.340 - const CostVector &_pi;
1.341 + const ValueVector &_pi;
1.342 int &_in_arc;
1.343 int _arc_num;
1.344
1.345 @@ -301,7 +292,7 @@
1.346
1.347 public:
1.348
1.349 - /// Constructor
1.350 + // Constructor
1.351 BlockSearchPivotRule(NetworkSimplex &ns) :
1.352 _source(ns._source), _target(ns._target),
1.353 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.354 @@ -315,9 +306,9 @@
1.355 MIN_BLOCK_SIZE );
1.356 }
1.357
1.358 - /// Find next entering arc
1.359 + // Find next entering arc
1.360 bool findEnteringArc() {
1.361 - Cost c, min = 0;
1.362 + Value c, min = 0;
1.363 int cnt = _block_size;
1.364 int e, min_arc = _next_arc;
1.365 for (e = _next_arc; e < _arc_num; ++e) {
1.366 @@ -353,13 +344,7 @@
1.367 }; //class BlockSearchPivotRule
1.368
1.369
1.370 - /// \brief Implementation of the "Candidate List" pivot rule for the
1.371 - /// \ref NetworkSimplex "network simplex" algorithm.
1.372 - ///
1.373 - /// This class implements the "Candidate List" pivot rule
1.374 - /// for the \ref NetworkSimplex "network simplex" algorithm.
1.375 - ///
1.376 - /// For more information see \ref NetworkSimplex::run().
1.377 + // Implementation of the Candidate List pivot rule
1.378 class CandidateListPivotRule
1.379 {
1.380 private:
1.381 @@ -367,9 +352,9 @@
1.382 // References to the NetworkSimplex class
1.383 const IntVector &_source;
1.384 const IntVector &_target;
1.385 - const CostVector &_cost;
1.386 + const ValueVector &_cost;
1.387 const IntVector &_state;
1.388 - const CostVector &_pi;
1.389 + const ValueVector &_pi;
1.390 int &_in_arc;
1.391 int _arc_num;
1.392
1.393 @@ -403,7 +388,7 @@
1.394
1.395 /// Find next entering arc
1.396 bool findEnteringArc() {
1.397 - Cost min, c;
1.398 + Value min, c;
1.399 int e, min_arc = _next_arc;
1.400 if (_curr_length > 0 && _minor_count < _minor_limit) {
1.401 // Minor iteration: select the best eligible arc from the
1.402 @@ -464,13 +449,7 @@
1.403 }; //class CandidateListPivotRule
1.404
1.405
1.406 - /// \brief Implementation of the "Altering Candidate List" pivot rule
1.407 - /// for the \ref NetworkSimplex "network simplex" algorithm.
1.408 - ///
1.409 - /// This class implements the "Altering Candidate List" pivot rule
1.410 - /// for the \ref NetworkSimplex "network simplex" algorithm.
1.411 - ///
1.412 - /// For more information see \ref NetworkSimplex::run().
1.413 + // Implementation of the Altering Candidate List pivot rule
1.414 class AlteringListPivotRule
1.415 {
1.416 private:
1.417 @@ -478,9 +457,9 @@
1.418 // References to the NetworkSimplex class
1.419 const IntVector &_source;
1.420 const IntVector &_target;
1.421 - const CostVector &_cost;
1.422 + const ValueVector &_cost;
1.423 const IntVector &_state;
1.424 - const CostVector &_pi;
1.425 + const ValueVector &_pi;
1.426 int &_in_arc;
1.427 int _arc_num;
1.428
1.429 @@ -488,15 +467,15 @@
1.430 int _block_size, _head_length, _curr_length;
1.431 int _next_arc;
1.432 IntVector _candidates;
1.433 - CostVector _cand_cost;
1.434 + ValueVector _cand_cost;
1.435
1.436 // Functor class to compare arcs during sort of the candidate list
1.437 class SortFunc
1.438 {
1.439 private:
1.440 - const CostVector &_map;
1.441 + const ValueVector &_map;
1.442 public:
1.443 - SortFunc(const CostVector &map) : _map(map) {}
1.444 + SortFunc(const ValueVector &map) : _map(map) {}
1.445 bool operator()(int left, int right) {
1.446 return _map[left] > _map[right];
1.447 }
1.448 @@ -506,7 +485,7 @@
1.449
1.450 public:
1.451
1.452 - /// Constructor
1.453 + // Constructor
1.454 AlteringListPivotRule(NetworkSimplex &ns) :
1.455 _source(ns._source), _target(ns._target),
1.456 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.457 @@ -527,7 +506,7 @@
1.458 _curr_length = 0;
1.459 }
1.460
1.461 - /// Find next entering arc
1.462 + // Find next entering arc
1.463 bool findEnteringArc() {
1.464 // Check the current candidate list
1.465 int e;
1.466 @@ -592,95 +571,23 @@
1.467
1.468 public:
1.469
1.470 - /// \brief General constructor (with lower bounds).
1.471 + /// \brief Constructor.
1.472 ///
1.473 - /// General constructor (with lower bounds).
1.474 + /// Constructor.
1.475 ///
1.476 /// \param graph The digraph the algorithm runs on.
1.477 - /// \param lower The lower bounds of the arcs.
1.478 - /// \param capacity The capacities (upper bounds) of the arcs.
1.479 - /// \param cost The cost (length) values of the arcs.
1.480 - /// \param supply The supply values of the nodes (signed).
1.481 - NetworkSimplex( const Digraph &graph,
1.482 - const LowerMap &lower,
1.483 - const CapacityMap &capacity,
1.484 - const CostMap &cost,
1.485 - const SupplyMap &supply ) :
1.486 - _graph(graph), _orig_lower(&lower), _orig_cap(capacity),
1.487 - _orig_cost(cost), _orig_supply(&supply),
1.488 + NetworkSimplex(const GR& graph) :
1.489 + _graph(graph),
1.490 + _plower(NULL), _pupper(NULL), _pcost(NULL),
1.491 + _psupply(NULL), _pstsup(false),
1.492 _flow_map(NULL), _potential_map(NULL),
1.493 _local_flow(false), _local_potential(false),
1.494 _node_id(graph)
1.495 - {}
1.496 -
1.497 - /// \brief General constructor (without lower bounds).
1.498 - ///
1.499 - /// General constructor (without lower bounds).
1.500 - ///
1.501 - /// \param graph The digraph the algorithm runs on.
1.502 - /// \param capacity The capacities (upper bounds) of the arcs.
1.503 - /// \param cost The cost (length) values of the arcs.
1.504 - /// \param supply The supply values of the nodes (signed).
1.505 - NetworkSimplex( const Digraph &graph,
1.506 - const CapacityMap &capacity,
1.507 - const CostMap &cost,
1.508 - const SupplyMap &supply ) :
1.509 - _graph(graph), _orig_lower(NULL), _orig_cap(capacity),
1.510 - _orig_cost(cost), _orig_supply(&supply),
1.511 - _flow_map(NULL), _potential_map(NULL),
1.512 - _local_flow(false), _local_potential(false),
1.513 - _node_id(graph)
1.514 - {}
1.515 -
1.516 - /// \brief Simple constructor (with lower bounds).
1.517 - ///
1.518 - /// Simple constructor (with lower bounds).
1.519 - ///
1.520 - /// \param graph The digraph the algorithm runs on.
1.521 - /// \param lower The lower bounds of the arcs.
1.522 - /// \param capacity The capacities (upper bounds) of the arcs.
1.523 - /// \param cost The cost (length) values of the arcs.
1.524 - /// \param s The source node.
1.525 - /// \param t The target node.
1.526 - /// \param flow_value The required amount of flow from node \c s
1.527 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.528 - NetworkSimplex( const Digraph &graph,
1.529 - const LowerMap &lower,
1.530 - const CapacityMap &capacity,
1.531 - const CostMap &cost,
1.532 - Node s, Node t,
1.533 - Capacity flow_value ) :
1.534 - _graph(graph), _orig_lower(&lower), _orig_cap(capacity),
1.535 - _orig_cost(cost), _orig_supply(NULL),
1.536 - _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
1.537 - _flow_map(NULL), _potential_map(NULL),
1.538 - _local_flow(false), _local_potential(false),
1.539 - _node_id(graph)
1.540 - {}
1.541 -
1.542 - /// \brief Simple constructor (without lower bounds).
1.543 - ///
1.544 - /// Simple constructor (without lower bounds).
1.545 - ///
1.546 - /// \param graph The digraph the algorithm runs on.
1.547 - /// \param capacity The capacities (upper bounds) of the arcs.
1.548 - /// \param cost The cost (length) values of the arcs.
1.549 - /// \param s The source node.
1.550 - /// \param t The target node.
1.551 - /// \param flow_value The required amount of flow from node \c s
1.552 - /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.553 - NetworkSimplex( const Digraph &graph,
1.554 - const CapacityMap &capacity,
1.555 - const CostMap &cost,
1.556 - Node s, Node t,
1.557 - Capacity flow_value ) :
1.558 - _graph(graph), _orig_lower(NULL), _orig_cap(capacity),
1.559 - _orig_cost(cost), _orig_supply(NULL),
1.560 - _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
1.561 - _flow_map(NULL), _potential_map(NULL),
1.562 - _local_flow(false), _local_potential(false),
1.563 - _node_id(graph)
1.564 - {}
1.565 + {
1.566 + LEMON_ASSERT(std::numeric_limits<Value>::is_integer &&
1.567 + std::numeric_limits<Value>::is_signed,
1.568 + "The value type of NetworkSimplex must be a signed integer");
1.569 + }
1.570
1.571 /// Destructor.
1.572 ~NetworkSimplex() {
1.573 @@ -688,12 +595,165 @@
1.574 if (_local_potential) delete _potential_map;
1.575 }
1.576
1.577 + /// \brief Set the lower bounds on the arcs.
1.578 + ///
1.579 + /// This function sets the lower bounds on the arcs.
1.580 + /// If neither this function nor \ref boundMaps() is used before
1.581 + /// calling \ref run(), the lower bounds will be set to zero
1.582 + /// on all arcs.
1.583 + ///
1.584 + /// \param map An arc map storing the lower bounds.
1.585 + /// Its \c Value type must be convertible to the \c Value type
1.586 + /// of the algorithm.
1.587 + ///
1.588 + /// \return <tt>(*this)</tt>
1.589 + template <typename LOWER>
1.590 + NetworkSimplex& lowerMap(const LOWER& map) {
1.591 + delete _plower;
1.592 + _plower = new ValueArcMap(_graph);
1.593 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.594 + (*_plower)[a] = map[a];
1.595 + }
1.596 + return *this;
1.597 + }
1.598 +
1.599 + /// \brief Set the upper bounds (capacities) on the arcs.
1.600 + ///
1.601 + /// This function sets the upper bounds (capacities) on the arcs.
1.602 + /// If none of the functions \ref upperMap(), \ref capacityMap()
1.603 + /// and \ref boundMaps() is used before calling \ref run(),
1.604 + /// the upper bounds (capacities) will be set to
1.605 + /// \c std::numeric_limits<Value>::max() on all arcs.
1.606 + ///
1.607 + /// \param map An arc map storing the upper bounds.
1.608 + /// Its \c Value type must be convertible to the \c Value type
1.609 + /// of the algorithm.
1.610 + ///
1.611 + /// \return <tt>(*this)</tt>
1.612 + template<typename UPPER>
1.613 + NetworkSimplex& upperMap(const UPPER& map) {
1.614 + delete _pupper;
1.615 + _pupper = new ValueArcMap(_graph);
1.616 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.617 + (*_pupper)[a] = map[a];
1.618 + }
1.619 + return *this;
1.620 + }
1.621 +
1.622 + /// \brief Set the upper bounds (capacities) on the arcs.
1.623 + ///
1.624 + /// This function sets the upper bounds (capacities) on the arcs.
1.625 + /// It is just an alias for \ref upperMap().
1.626 + ///
1.627 + /// \return <tt>(*this)</tt>
1.628 + template<typename CAP>
1.629 + NetworkSimplex& capacityMap(const CAP& map) {
1.630 + return upperMap(map);
1.631 + }
1.632 +
1.633 + /// \brief Set the lower and upper bounds on the arcs.
1.634 + ///
1.635 + /// This function sets the lower and upper bounds on the arcs.
1.636 + /// If neither this function nor \ref lowerMap() is used before
1.637 + /// calling \ref run(), the lower bounds will be set to zero
1.638 + /// on all arcs.
1.639 + /// If none of the functions \ref upperMap(), \ref capacityMap()
1.640 + /// and \ref boundMaps() is used before calling \ref run(),
1.641 + /// the upper bounds (capacities) will be set to
1.642 + /// \c std::numeric_limits<Value>::max() on all arcs.
1.643 + ///
1.644 + /// \param lower An arc map storing the lower bounds.
1.645 + /// \param upper An arc map storing the upper bounds.
1.646 + ///
1.647 + /// The \c Value type of the maps must be convertible to the
1.648 + /// \c Value type of the algorithm.
1.649 + ///
1.650 + /// \note This function is just a shortcut of calling \ref lowerMap()
1.651 + /// and \ref upperMap() separately.
1.652 + ///
1.653 + /// \return <tt>(*this)</tt>
1.654 + template <typename LOWER, typename UPPER>
1.655 + NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
1.656 + return lowerMap(lower).upperMap(upper);
1.657 + }
1.658 +
1.659 + /// \brief Set the costs of the arcs.
1.660 + ///
1.661 + /// This function sets the costs of the arcs.
1.662 + /// If it is not used before calling \ref run(), the costs
1.663 + /// will be set to \c 1 on all arcs.
1.664 + ///
1.665 + /// \param map An arc map storing the costs.
1.666 + /// Its \c Value type must be convertible to the \c Value type
1.667 + /// of the algorithm.
1.668 + ///
1.669 + /// \return <tt>(*this)</tt>
1.670 + template<typename COST>
1.671 + NetworkSimplex& costMap(const COST& map) {
1.672 + delete _pcost;
1.673 + _pcost = new ValueArcMap(_graph);
1.674 + for (ArcIt a(_graph); a != INVALID; ++a) {
1.675 + (*_pcost)[a] = map[a];
1.676 + }
1.677 + return *this;
1.678 + }
1.679 +
1.680 + /// \brief Set the supply values of the nodes.
1.681 + ///
1.682 + /// This function sets the supply values of the nodes.
1.683 + /// If neither this function nor \ref stSupply() is used before
1.684 + /// calling \ref run(), the supply of each node will be set to zero.
1.685 + /// (It makes sense only if non-zero lower bounds are given.)
1.686 + ///
1.687 + /// \param map A node map storing the supply values.
1.688 + /// Its \c Value type must be convertible to the \c Value type
1.689 + /// of the algorithm.
1.690 + ///
1.691 + /// \return <tt>(*this)</tt>
1.692 + template<typename SUP>
1.693 + NetworkSimplex& supplyMap(const SUP& map) {
1.694 + delete _psupply;
1.695 + _pstsup = false;
1.696 + _psupply = new ValueNodeMap(_graph);
1.697 + for (NodeIt n(_graph); n != INVALID; ++n) {
1.698 + (*_psupply)[n] = map[n];
1.699 + }
1.700 + return *this;
1.701 + }
1.702 +
1.703 + /// \brief Set single source and target nodes and a supply value.
1.704 + ///
1.705 + /// This function sets a single source node and a single target node
1.706 + /// and the required flow value.
1.707 + /// If neither this function nor \ref supplyMap() is used before
1.708 + /// calling \ref run(), the supply of each node will be set to zero.
1.709 + /// (It makes sense only if non-zero lower bounds are given.)
1.710 + ///
1.711 + /// \param s The source node.
1.712 + /// \param t The target node.
1.713 + /// \param k The required amount of flow from node \c s to node \c t
1.714 + /// (i.e. the supply of \c s and the demand of \c t).
1.715 + ///
1.716 + /// \return <tt>(*this)</tt>
1.717 + NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
1.718 + delete _psupply;
1.719 + _psupply = NULL;
1.720 + _pstsup = true;
1.721 + _psource = s;
1.722 + _ptarget = t;
1.723 + _pstflow = k;
1.724 + return *this;
1.725 + }
1.726 +
1.727 /// \brief Set the flow map.
1.728 ///
1.729 /// This function sets the flow map.
1.730 + /// If it is not used before calling \ref run(), an instance will
1.731 + /// be allocated automatically. The destructor deallocates this
1.732 + /// automatically allocated map, of course.
1.733 ///
1.734 /// \return <tt>(*this)</tt>
1.735 - NetworkSimplex& flowMap(FlowMap &map) {
1.736 + NetworkSimplex& flowMap(FlowMap& map) {
1.737 if (_local_flow) {
1.738 delete _flow_map;
1.739 _local_flow = false;
1.740 @@ -704,10 +764,14 @@
1.741
1.742 /// \brief Set the potential map.
1.743 ///
1.744 - /// This function sets the potential map.
1.745 + /// This function sets the potential map, which is used for storing
1.746 + /// the dual solution.
1.747 + /// If it is not used before calling \ref run(), an instance will
1.748 + /// be allocated automatically. The destructor deallocates this
1.749 + /// automatically allocated map, of course.
1.750 ///
1.751 /// \return <tt>(*this)</tt>
1.752 - NetworkSimplex& potentialMap(PotentialMap &map) {
1.753 + NetworkSimplex& potentialMap(PotentialMap& map) {
1.754 if (_local_potential) {
1.755 delete _potential_map;
1.756 _local_potential = false;
1.757 @@ -716,47 +780,29 @@
1.758 return *this;
1.759 }
1.760
1.761 - /// \name Execution control
1.762 - /// The algorithm can be executed using the
1.763 - /// \ref lemon::NetworkSimplex::run() "run()" function.
1.764 + /// \name Execution Control
1.765 + /// The algorithm can be executed using \ref run().
1.766 +
1.767 /// @{
1.768
1.769 /// \brief Run the algorithm.
1.770 ///
1.771 /// This function runs the algorithm.
1.772 + /// The paramters can be specified using \ref lowerMap(),
1.773 + /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
1.774 + /// \ref costMap(), \ref supplyMap() and \ref stSupply()
1.775 + /// functions. For example,
1.776 + /// \code
1.777 + /// NetworkSimplex<ListDigraph> ns(graph);
1.778 + /// ns.boundMaps(lower, upper).costMap(cost)
1.779 + /// .supplyMap(sup).run();
1.780 + /// \endcode
1.781 ///
1.782 - /// \param pivot_rule The pivot rule that is used during the
1.783 - /// algorithm.
1.784 - ///
1.785 - /// The available pivot rules:
1.786 - ///
1.787 - /// - FIRST_ELIGIBLE_PIVOT The next eligible arc is selected in
1.788 - /// a wraparound fashion in every iteration
1.789 - /// (\ref FirstEligiblePivotRule).
1.790 - ///
1.791 - /// - BEST_ELIGIBLE_PIVOT The best eligible arc is selected in
1.792 - /// every iteration (\ref BestEligiblePivotRule).
1.793 - ///
1.794 - /// - BLOCK_SEARCH_PIVOT A specified number of arcs are examined in
1.795 - /// every iteration in a wraparound fashion and the best eligible
1.796 - /// arc is selected from this block (\ref BlockSearchPivotRule).
1.797 - ///
1.798 - /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
1.799 - /// built from eligible arcs in a wraparound fashion and in the
1.800 - /// following minor iterations the best eligible arc is selected
1.801 - /// from this list (\ref CandidateListPivotRule).
1.802 - ///
1.803 - /// - ALTERING_LIST_PIVOT It is a modified version of the
1.804 - /// "Candidate List" pivot rule. It keeps only the several best
1.805 - /// eligible arcs from the former candidate list and extends this
1.806 - /// list in every iteration (\ref AlteringListPivotRule).
1.807 - ///
1.808 - /// According to our comprehensive benchmark tests the "Block Search"
1.809 - /// pivot rule proved to be the fastest and the most robust on
1.810 - /// various test inputs. Thus it is the default option.
1.811 + /// \param pivot_rule The pivot rule that will be used during the
1.812 + /// algorithm. For more information see \ref PivotRule.
1.813 ///
1.814 /// \return \c true if a feasible flow can be found.
1.815 - bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
1.816 + bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
1.817 return init() && start(pivot_rule);
1.818 }
1.819
1.820 @@ -765,10 +811,53 @@
1.821 /// \name Query Functions
1.822 /// The results of the algorithm can be obtained using these
1.823 /// functions.\n
1.824 - /// \ref lemon::NetworkSimplex::run() "run()" must be called before
1.825 - /// using them.
1.826 + /// The \ref run() function must be called before using them.
1.827 +
1.828 /// @{
1.829
1.830 + /// \brief Return the total cost of the found flow.
1.831 + ///
1.832 + /// This function returns the total cost of the found flow.
1.833 + /// The complexity of the function is \f$ O(e) \f$.
1.834 + ///
1.835 + /// \note The return type of the function can be specified as a
1.836 + /// template parameter. For example,
1.837 + /// \code
1.838 + /// ns.totalCost<double>();
1.839 + /// \endcode
1.840 + /// It is useful if the total cost cannot be stored in the \c Value
1.841 + /// type of the algorithm, which is the default return type of the
1.842 + /// function.
1.843 + ///
1.844 + /// \pre \ref run() must be called before using this function.
1.845 + template <typename Num>
1.846 + Num totalCost() const {
1.847 + Num c = 0;
1.848 + if (_pcost) {
1.849 + for (ArcIt e(_graph); e != INVALID; ++e)
1.850 + c += (*_flow_map)[e] * (*_pcost)[e];
1.851 + } else {
1.852 + for (ArcIt e(_graph); e != INVALID; ++e)
1.853 + c += (*_flow_map)[e];
1.854 + }
1.855 + return c;
1.856 + }
1.857 +
1.858 +#ifndef DOXYGEN
1.859 + Value totalCost() const {
1.860 + return totalCost<Value>();
1.861 + }
1.862 +#endif
1.863 +
1.864 + /// \brief Return the flow on the given arc.
1.865 + ///
1.866 + /// This function returns the flow on the given arc.
1.867 + ///
1.868 + /// \pre \ref run() must be called before using this function.
1.869 + Value flow(const Arc& a) const {
1.870 + return (*_flow_map)[a];
1.871 + }
1.872 +
1.873 /// \brief Return a const reference to the flow map.
1.874 ///
1.875 /// This function returns a const reference to an arc map storing
1.876 @@ -779,48 +868,28 @@
1.877 return *_flow_map;
1.878 }
1.879
1.880 + /// \brief Return the potential (dual value) of the given node.
1.881 + ///
1.882 + /// This function returns the potential (dual value) of the
1.883 + /// given node.
1.884 + ///
1.885 + /// \pre \ref run() must be called before using this function.
1.886 + Value potential(const Node& n) const {
1.887 + return (*_potential_map)[n];
1.888 + }
1.889 +
1.890 /// \brief Return a const reference to the potential map
1.891 /// (the dual solution).
1.892 ///
1.893 /// This function returns a const reference to a node map storing
1.894 - /// the found potentials (the dual solution).
1.895 + /// the found potentials, which form the dual solution of the
1.896 + /// \ref min_cost_flow "minimum cost flow" problem.
1.897 ///
1.898 /// \pre \ref run() must be called before using this function.
1.899 const PotentialMap& potentialMap() const {
1.900 return *_potential_map;
1.901 }
1.902
1.903 - /// \brief Return the flow on the given arc.
1.904 - ///
1.905 - /// This function returns the flow on the given arc.
1.906 - ///
1.907 - /// \pre \ref run() must be called before using this function.
1.908 - Capacity flow(const Arc& arc) const {
1.909 - return (*_flow_map)[arc];
1.910 - }
1.911 -
1.912 - /// \brief Return the potential of the given node.
1.913 - ///
1.914 - /// This function returns the potential of the given node.
1.915 - ///
1.916 - /// \pre \ref run() must be called before using this function.
1.917 - Cost potential(const Node& node) const {
1.918 - return (*_potential_map)[node];
1.919 - }
1.920 -
1.921 - /// \brief Return the total cost of the found flow.
1.922 - ///
1.923 - /// This function returns the total cost of the found flow.
1.924 - /// The complexity of the function is \f$ O(e) \f$.
1.925 - ///
1.926 - /// \pre \ref run() must be called before using this function.
1.927 - Cost totalCost() const {
1.928 - Cost c = 0;
1.929 - for (ArcIt e(_graph); e != INVALID; ++e)
1.930 - c += (*_flow_map)[e] * _orig_cost[e];
1.931 - return c;
1.932 - }
1.933 -
1.934 /// @}
1.935
1.936 private:
1.937 @@ -842,6 +911,7 @@
1.938 _arc_num = countArcs(_graph);
1.939 int all_node_num = _node_num + 1;
1.940 int all_arc_num = _arc_num + _node_num;
1.941 + if (_node_num == 0) return false;
1.942
1.943 _arc_ref.resize(_arc_num);
1.944 _source.resize(all_arc_num);
1.945 @@ -864,12 +934,17 @@
1.946
1.947 // Initialize node related data
1.948 bool valid_supply = true;
1.949 - if (_orig_supply) {
1.950 - Supply sum = 0;
1.951 + if (!_pstsup && !_psupply) {
1.952 + _pstsup = true;
1.953 + _psource = _ptarget = NodeIt(_graph);
1.954 + _pstflow = 0;
1.955 + }
1.956 + if (_psupply) {
1.957 + Value sum = 0;
1.958 int i = 0;
1.959 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.960 _node_id[n] = i;
1.961 - _supply[i] = (*_orig_supply)[n];
1.962 + _supply[i] = (*_psupply)[n];
1.963 sum += _supply[i];
1.964 }
1.965 valid_supply = (sum == 0);
1.966 @@ -879,8 +954,8 @@
1.967 _node_id[n] = i;
1.968 _supply[i] = 0;
1.969 }
1.970 - _supply[_node_id[_orig_source]] = _orig_flow_value;
1.971 - _supply[_node_id[_orig_target]] = -_orig_flow_value;
1.972 + _supply[_node_id[_psource]] = _pstflow;
1.973 + _supply[_node_id[_ptarget]] = -_pstflow;
1.974 }
1.975 if (!valid_supply) return false;
1.976
1.977 @@ -904,18 +979,41 @@
1.978 }
1.979
1.980 // Initialize arc maps
1.981 - for (int i = 0; i != _arc_num; ++i) {
1.982 - Arc e = _arc_ref[i];
1.983 - _source[i] = _node_id[_graph.source(e)];
1.984 - _target[i] = _node_id[_graph.target(e)];
1.985 - _cost[i] = _orig_cost[e];
1.986 - _cap[i] = _orig_cap[e];
1.987 + if (_pupper && _pcost) {
1.988 + for (int i = 0; i != _arc_num; ++i) {
1.989 + Arc e = _arc_ref[i];
1.990 + _source[i] = _node_id[_graph.source(e)];
1.991 + _target[i] = _node_id[_graph.target(e)];
1.992 + _cap[i] = (*_pupper)[e];
1.993 + _cost[i] = (*_pcost)[e];
1.994 + }
1.995 + } else {
1.996 + for (int i = 0; i != _arc_num; ++i) {
1.997 + Arc e = _arc_ref[i];
1.998 + _source[i] = _node_id[_graph.source(e)];
1.999 + _target[i] = _node_id[_graph.target(e)];
1.1000 + }
1.1001 + if (_pupper) {
1.1002 + for (int i = 0; i != _arc_num; ++i)
1.1003 + _cap[i] = (*_pupper)[_arc_ref[i]];
1.1004 + } else {
1.1005 + Value val = std::numeric_limits<Value>::max();
1.1006 + for (int i = 0; i != _arc_num; ++i)
1.1007 + _cap[i] = val;
1.1008 + }
1.1009 + if (_pcost) {
1.1010 + for (int i = 0; i != _arc_num; ++i)
1.1011 + _cost[i] = (*_pcost)[_arc_ref[i]];
1.1012 + } else {
1.1013 + for (int i = 0; i != _arc_num; ++i)
1.1014 + _cost[i] = 1;
1.1015 + }
1.1016 }
1.1017
1.1018 // Remove non-zero lower bounds
1.1019 - if (_orig_lower) {
1.1020 + if (_plower) {
1.1021 for (int i = 0; i != _arc_num; ++i) {
1.1022 - Capacity c = (*_orig_lower)[_arc_ref[i]];
1.1023 + Value c = (*_plower)[_arc_ref[i]];
1.1024 if (c != 0) {
1.1025 _cap[i] -= c;
1.1026 _supply[_source[i]] -= c;
1.1027 @@ -925,8 +1023,8 @@
1.1028 }
1.1029
1.1030 // Add artificial arcs and initialize the spanning tree data structure
1.1031 - Cost max_cost = std::numeric_limits<Cost>::max() / 4;
1.1032 - Capacity max_cap = std::numeric_limits<Capacity>::max();
1.1033 + Value max_cap = std::numeric_limits<Value>::max();
1.1034 + Value max_cost = std::numeric_limits<Value>::max() / 4;
1.1035 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1.1036 _thread[u] = u + 1;
1.1037 _rev_thread[u + 1] = u;
1.1038 @@ -979,7 +1077,7 @@
1.1039 }
1.1040 delta = _cap[in_arc];
1.1041 int result = 0;
1.1042 - Capacity d;
1.1043 + Value d;
1.1044 int e;
1.1045
1.1046 // Search the cycle along the path form the first node to the root
1.1047 @@ -1017,7 +1115,7 @@
1.1048 void changeFlow(bool change) {
1.1049 // Augment along the cycle
1.1050 if (delta > 0) {
1.1051 - Capacity val = _state[in_arc] * delta;
1.1052 + Value val = _state[in_arc] * delta;
1.1053 _flow[in_arc] += val;
1.1054 for (int u = _source[in_arc]; u != join; u = _parent[u]) {
1.1055 _flow[_pred[u]] += _forward[u] ? -val : val;
1.1056 @@ -1158,7 +1256,7 @@
1.1057
1.1058 // Update potentials
1.1059 void updatePotential() {
1.1060 - Cost sigma = _forward[u_in] ?
1.1061 + Value sigma = _forward[u_in] ?
1.1062 _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1.1063 _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1.1064 if (_succ_num[u_in] > _node_num / 2) {
1.1065 @@ -1181,28 +1279,28 @@
1.1066 }
1.1067
1.1068 // Execute the algorithm
1.1069 - bool start(PivotRuleEnum pivot_rule) {
1.1070 + bool start(PivotRule pivot_rule) {
1.1071 // Select the pivot rule implementation
1.1072 switch (pivot_rule) {
1.1073 - case FIRST_ELIGIBLE_PIVOT:
1.1074 + case FIRST_ELIGIBLE:
1.1075 return start<FirstEligiblePivotRule>();
1.1076 - case BEST_ELIGIBLE_PIVOT:
1.1077 + case BEST_ELIGIBLE:
1.1078 return start<BestEligiblePivotRule>();
1.1079 - case BLOCK_SEARCH_PIVOT:
1.1080 + case BLOCK_SEARCH:
1.1081 return start<BlockSearchPivotRule>();
1.1082 - case CANDIDATE_LIST_PIVOT:
1.1083 + case CANDIDATE_LIST:
1.1084 return start<CandidateListPivotRule>();
1.1085 - case ALTERING_LIST_PIVOT:
1.1086 + case ALTERING_LIST:
1.1087 return start<AlteringListPivotRule>();
1.1088 }
1.1089 return false;
1.1090 }
1.1091
1.1092 - template<class PivotRuleImplementation>
1.1093 + template <typename PivotRuleImpl>
1.1094 bool start() {
1.1095 - PivotRuleImplementation pivot(*this);
1.1096 + PivotRuleImpl pivot(*this);
1.1097
1.1098 - // Execute the network simplex algorithm
1.1099 + // Execute the Network Simplex algorithm
1.1100 while (pivot.findEnteringArc()) {
1.1101 findJoinNode();
1.1102 bool change = findLeavingArc();
1.1103 @@ -1219,10 +1317,10 @@
1.1104 }
1.1105
1.1106 // Copy flow values to _flow_map
1.1107 - if (_orig_lower) {
1.1108 + if (_plower) {
1.1109 for (int i = 0; i != _arc_num; ++i) {
1.1110 Arc e = _arc_ref[i];
1.1111 - _flow_map->set(e, (*_orig_lower)[e] + _flow[i]);
1.1112 + _flow_map->set(e, (*_plower)[e] + _flow[i]);
1.1113 }
1.1114 } else {
1.1115 for (int i = 0; i != _arc_num; ++i) {