1.1 --- a/lemon/network_simplex.h Fri Sep 19 15:14:41 2008 +0000
1.2 +++ b/lemon/network_simplex.h Sun Oct 05 13:36:43 2008 +0000
1.3 @@ -32,16 +32,18 @@
1.4 #include <lemon/smart_graph.h>
1.5 #include <lemon/math.h>
1.6
1.7 +#define _DEBUG_
1.8 +
1.9 namespace lemon {
1.10
1.11 /// \addtogroup min_cost_flow
1.12 /// @{
1.13
1.14 - /// \brief Implementation of the network simplex algorithm for
1.15 - /// finding a minimum cost flow.
1.16 + /// \brief Implementation of the primal network simplex algorithm
1.17 + /// for finding a minimum cost flow.
1.18 ///
1.19 - /// \ref NetworkSimplex implements the network simplex algorithm for
1.20 - /// finding a minimum cost flow.
1.21 + /// \ref NetworkSimplex implements the primal network simplex algorithm
1.22 + /// for finding a minimum cost flow.
1.23 ///
1.24 /// \tparam Graph The directed graph type the algorithm runs on.
1.25 /// \tparam LowerMap The type of the lower bound map.
1.26 @@ -55,16 +57,15 @@
1.27 /// - The value types of the maps should be convertible to each other.
1.28 /// - \c CostMap::Value must be signed type.
1.29 ///
1.30 - /// \note \ref NetworkSimplex provides six different pivot rule
1.31 + /// \note \ref NetworkSimplex provides five different pivot rule
1.32 /// implementations that significantly affect the efficiency of the
1.33 /// algorithm.
1.34 - /// By default a combined pivot rule is used, which is the fastest
1.35 - /// implementation according to our benchmark tests.
1.36 - /// Another pivot rule can be selected using \ref run() function
1.37 - /// with the proper parameter.
1.38 + /// By default "Block Search" pivot rule is used, which proved to be
1.39 + /// by far the most efficient according to our benchmark tests.
1.40 + /// However another pivot rule can be selected using \ref run()
1.41 + /// function with the proper parameter.
1.42 ///
1.43 /// \author Peter Kovacs
1.44 -
1.45 template < typename Graph,
1.46 typename LowerMap = typename Graph::template EdgeMap<int>,
1.47 typename CapacityMap = typename Graph::template EdgeMap<int>,
1.48 @@ -89,10 +90,13 @@
1.49 typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
1.50 typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
1.51 typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
1.52 + typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;
1.53
1.54 typedef typename Graph::template NodeMap<Node> NodeRefMap;
1.55 typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
1.56
1.57 + typedef std::vector<Edge> EdgeVector;
1.58 +
1.59 public:
1.60
1.61 /// The type of the flow map.
1.62 @@ -107,9 +111,8 @@
1.63 FIRST_ELIGIBLE_PIVOT,
1.64 BEST_ELIGIBLE_PIVOT,
1.65 BLOCK_SEARCH_PIVOT,
1.66 - LIMITED_SEARCH_PIVOT,
1.67 CANDIDATE_LIST_PIVOT,
1.68 - COMBINED_PIVOT
1.69 + ALTERING_LIST_PIVOT
1.70 };
1.71
1.72 private:
1.73 @@ -148,32 +151,40 @@
1.74 ///
1.75 /// This class implements the "First Eligible" pivot rule
1.76 /// for the \ref NetworkSimplex "network simplex" algorithm.
1.77 + ///
1.78 + /// For more information see \ref NetworkSimplex::run().
1.79 class FirstEligiblePivotRule
1.80 {
1.81 private:
1.82
1.83 + // References to the NetworkSimplex class
1.84 NetworkSimplex &_ns;
1.85 - EdgeIt _next_edge;
1.86 + EdgeVector &_edges;
1.87 +
1.88 + int _next_edge;
1.89
1.90 public:
1.91
1.92 - /// Constructor.
1.93 - FirstEligiblePivotRule(NetworkSimplex &ns) :
1.94 - _ns(ns), _next_edge(ns._graph) {}
1.95 + /// Constructor
1.96 + FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.97 + _ns(ns), _edges(edges), _next_edge(0) {}
1.98
1.99 - /// Finds the next entering edge.
1.100 - bool findEnteringEdge() {
1.101 - for (EdgeIt e = _next_edge; e != INVALID; ++e) {
1.102 + /// Find next entering edge
1.103 + inline bool findEnteringEdge() {
1.104 + Edge e;
1.105 + for (int i = _next_edge; i < int(_edges.size()); ++i) {
1.106 + e = _edges[i];
1.107 if (_ns._state[e] * _ns._red_cost[e] < 0) {
1.108 _ns._in_edge = e;
1.109 - _next_edge = ++e;
1.110 + _next_edge = i + 1;
1.111 return true;
1.112 }
1.113 }
1.114 - for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
1.115 + for (int i = 0; i < _next_edge; ++i) {
1.116 + e = _edges[i];
1.117 if (_ns._state[e] * _ns._red_cost[e] < 0) {
1.118 _ns._in_edge = e;
1.119 - _next_edge = ++e;
1.120 + _next_edge = i + 1;
1.121 return true;
1.122 }
1.123 }
1.124 @@ -186,21 +197,28 @@
1.125 ///
1.126 /// This class implements the "Best Eligible" pivot rule
1.127 /// for the \ref NetworkSimplex "network simplex" algorithm.
1.128 + ///
1.129 + /// For more information see \ref NetworkSimplex::run().
1.130 class BestEligiblePivotRule
1.131 {
1.132 private:
1.133
1.134 + // References to the NetworkSimplex class
1.135 NetworkSimplex &_ns;
1.136 + EdgeVector &_edges;
1.137
1.138 public:
1.139
1.140 - /// Constructor.
1.141 - BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
1.142 + /// Constructor
1.143 + BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.144 + _ns(ns), _edges(edges) {}
1.145
1.146 - /// Finds the next entering edge.
1.147 - bool findEnteringEdge() {
1.148 + /// Find next entering edge
1.149 + inline bool findEnteringEdge() {
1.150 Cost min = 0;
1.151 - for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
1.152 + Edge e;
1.153 + for (int i = 0; i < int(_edges.size()); ++i) {
1.154 + e = _edges[i];
1.155 if (_ns._state[e] * _ns._red_cost[e] < min) {
1.156 min = _ns._state[e] * _ns._red_cost[e];
1.157 _ns._in_edge = e;
1.158 @@ -215,206 +233,291 @@
1.159 ///
1.160 /// This class implements the "Block Search" pivot rule
1.161 /// for the \ref NetworkSimplex "network simplex" algorithm.
1.162 + ///
1.163 + /// For more information see \ref NetworkSimplex::run().
1.164 class BlockSearchPivotRule
1.165 {
1.166 private:
1.167
1.168 + // References to the NetworkSimplex class
1.169 NetworkSimplex &_ns;
1.170 - EdgeIt _next_edge, _min_edge;
1.171 + EdgeVector &_edges;
1.172 +
1.173 int _block_size;
1.174 -
1.175 - static const int MIN_BLOCK_SIZE = 10;
1.176 + int _next_edge, _min_edge;
1.177
1.178 public:
1.179
1.180 - /// Constructor.
1.181 - BlockSearchPivotRule(NetworkSimplex &ns) :
1.182 - _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
1.183 + /// Constructor
1.184 + BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.185 + _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
1.186 {
1.187 - _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
1.188 - if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
1.189 + // The main parameters of the pivot rule
1.190 + const double BLOCK_SIZE_FACTOR = 2.0;
1.191 + const int MIN_BLOCK_SIZE = 10;
1.192 +
1.193 + _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
1.194 + MIN_BLOCK_SIZE );
1.195 }
1.196
1.197 - /// Finds the next entering edge.
1.198 - bool findEnteringEdge() {
1.199 + /// Find next entering edge
1.200 + inline bool findEnteringEdge() {
1.201 Cost curr, min = 0;
1.202 - int cnt = 0;
1.203 - for (EdgeIt e = _next_edge; e != INVALID; ++e) {
1.204 + Edge e;
1.205 + int cnt = _block_size;
1.206 + int i;
1.207 + for (i = _next_edge; i < int(_edges.size()); ++i) {
1.208 + e = _edges[i];
1.209 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.210 min = curr;
1.211 - _min_edge = e;
1.212 + _min_edge = i;
1.213 }
1.214 - if (++cnt == _block_size) {
1.215 + if (--cnt == 0) {
1.216 if (min < 0) break;
1.217 - cnt = 0;
1.218 + cnt = _block_size;
1.219 }
1.220 }
1.221 - if (min == 0) {
1.222 - for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
1.223 + if (min == 0 || cnt > 0) {
1.224 + for (i = 0; i < _next_edge; ++i) {
1.225 + e = _edges[i];
1.226 if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.227 min = curr;
1.228 - _min_edge = e;
1.229 + _min_edge = i;
1.230 }
1.231 - if (++cnt == _block_size) {
1.232 + if (--cnt == 0) {
1.233 if (min < 0) break;
1.234 - cnt = 0;
1.235 + cnt = _block_size;
1.236 }
1.237 }
1.238 }
1.239 - _ns._in_edge = _min_edge;
1.240 - _next_edge = ++_min_edge;
1.241 - return min < 0;
1.242 + if (min >= 0) return false;
1.243 + _ns._in_edge = _edges[_min_edge];
1.244 + _next_edge = i;
1.245 + return true;
1.246 }
1.247 }; //class BlockSearchPivotRule
1.248
1.249 - /// \brief Implementation of the "Limited Search" pivot rule for the
1.250 - /// \ref NetworkSimplex "network simplex" algorithm.
1.251 - ///
1.252 - /// This class implements the "Limited Search" pivot rule
1.253 - /// for the \ref NetworkSimplex "network simplex" algorithm.
1.254 - class LimitedSearchPivotRule
1.255 - {
1.256 - private:
1.257 -
1.258 - NetworkSimplex &_ns;
1.259 - EdgeIt _next_edge, _min_edge;
1.260 - int _sample_size;
1.261 -
1.262 - static const int SAMPLE_SIZE_FACTOR = 15;
1.263 - static const int MIN_SAMPLE_SIZE = 10;
1.264 -
1.265 - public:
1.266 -
1.267 - /// Constructor.
1.268 - LimitedSearchPivotRule(NetworkSimplex &ns) :
1.269 - _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
1.270 - {
1.271 - _sample_size = countEdges(_ns._graph) *
1.272 - SAMPLE_SIZE_FACTOR / 10000;
1.273 - if (_sample_size < MIN_SAMPLE_SIZE)
1.274 - _sample_size = MIN_SAMPLE_SIZE;
1.275 - }
1.276 -
1.277 - /// Finds the next entering edge.
1.278 - bool findEnteringEdge() {
1.279 - Cost curr, min = 0;
1.280 - int cnt = 0;
1.281 - for (EdgeIt e = _next_edge; e != INVALID; ++e) {
1.282 - if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.283 - min = curr;
1.284 - _min_edge = e;
1.285 - }
1.286 - if (curr < 0 && ++cnt == _sample_size) break;
1.287 - }
1.288 - if (min == 0) {
1.289 - for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
1.290 - if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.291 - min = curr;
1.292 - _min_edge = e;
1.293 - }
1.294 - if (curr < 0 && ++cnt == _sample_size) break;
1.295 - }
1.296 - }
1.297 - _ns._in_edge = _min_edge;
1.298 - _next_edge = ++_min_edge;
1.299 - return min < 0;
1.300 - }
1.301 - }; //class LimitedSearchPivotRule
1.302 -
1.303 /// \brief Implementation of the "Candidate List" pivot rule for the
1.304 /// \ref NetworkSimplex "network simplex" algorithm.
1.305 ///
1.306 /// This class implements the "Candidate List" pivot rule
1.307 /// for the \ref NetworkSimplex "network simplex" algorithm.
1.308 + ///
1.309 + /// For more information see \ref NetworkSimplex::run().
1.310 class CandidateListPivotRule
1.311 {
1.312 private:
1.313
1.314 + // References to the NetworkSimplex class
1.315 NetworkSimplex &_ns;
1.316 + EdgeVector &_edges;
1.317
1.318 - // The list of candidate edges.
1.319 - std::vector<Edge> _candidates;
1.320 - // The maximum length of the edge list.
1.321 - int _list_length;
1.322 - // The maximum number of minor iterations between two major
1.323 - // itarations.
1.324 - int _minor_limit;
1.325 -
1.326 - int _minor_count;
1.327 - EdgeIt _next_edge;
1.328 -
1.329 - static const int LIST_LENGTH_FACTOR = 20;
1.330 - static const int MINOR_LIMIT_FACTOR = 10;
1.331 - static const int MIN_LIST_LENGTH = 10;
1.332 - static const int MIN_MINOR_LIMIT = 2;
1.333 + EdgeVector _candidates;
1.334 + int _list_length, _minor_limit;
1.335 + int _curr_length, _minor_count;
1.336 + int _next_edge, _min_edge;
1.337
1.338 public:
1.339
1.340 - /// Constructor.
1.341 - CandidateListPivotRule(NetworkSimplex &ns) :
1.342 - _ns(ns), _next_edge(ns._graph)
1.343 + /// Constructor
1.344 + CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.345 + _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
1.346 {
1.347 - int edge_num = countEdges(_ns._graph);
1.348 - _minor_count = 0;
1.349 - _list_length = edge_num * LIST_LENGTH_FACTOR / 10000;
1.350 - if (_list_length < MIN_LIST_LENGTH)
1.351 - _list_length = MIN_LIST_LENGTH;
1.352 - _minor_limit = _list_length * MINOR_LIMIT_FACTOR / 100;
1.353 - if (_minor_limit < MIN_MINOR_LIMIT)
1.354 - _minor_limit = MIN_MINOR_LIMIT;
1.355 + // The main parameters of the pivot rule
1.356 + const double LIST_LENGTH_FACTOR = 1.0;
1.357 + const int MIN_LIST_LENGTH = 10;
1.358 + const double MINOR_LIMIT_FACTOR = 0.1;
1.359 + const int MIN_MINOR_LIMIT = 3;
1.360 +
1.361 + _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())),
1.362 + MIN_LIST_LENGTH );
1.363 + _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
1.364 + MIN_MINOR_LIMIT );
1.365 + _curr_length = _minor_count = 0;
1.366 + _candidates.resize(_list_length);
1.367 }
1.368
1.369 - /// Finds the next entering edge.
1.370 - bool findEnteringEdge() {
1.371 + /// Find next entering edge
1.372 + inline bool findEnteringEdge() {
1.373 Cost min, curr;
1.374 - if (_minor_count < _minor_limit && _candidates.size() > 0) {
1.375 - // Minor iteration
1.376 + if (_curr_length > 0 && _minor_count < _minor_limit) {
1.377 + // Minor iteration: selecting the best eligible edge from
1.378 + // the current candidate list
1.379 ++_minor_count;
1.380 Edge e;
1.381 min = 0;
1.382 - for (int i = 0; i < int(_candidates.size()); ++i) {
1.383 + for (int i = 0; i < _curr_length; ++i) {
1.384 e = _candidates[i];
1.385 - if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.386 + curr = _ns._state[e] * _ns._red_cost[e];
1.387 + if (curr < min) {
1.388 min = curr;
1.389 _ns._in_edge = e;
1.390 }
1.391 + if (curr >= 0) {
1.392 + _candidates[i--] = _candidates[--_curr_length];
1.393 + }
1.394 }
1.395 if (min < 0) return true;
1.396 }
1.397
1.398 - // Major iteration
1.399 - _candidates.clear();
1.400 - EdgeIt e = _next_edge;
1.401 + // Major iteration: building a new candidate list
1.402 + Edge e;
1.403 min = 0;
1.404 - for ( ; e != INVALID; ++e) {
1.405 + _curr_length = 0;
1.406 + int i;
1.407 + for (i = _next_edge; i < int(_edges.size()); ++i) {
1.408 + e = _edges[i];
1.409 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.410 - _candidates.push_back(e);
1.411 + _candidates[_curr_length++] = e;
1.412 if (curr < min) {
1.413 min = curr;
1.414 - _ns._in_edge = e;
1.415 + _min_edge = i;
1.416 }
1.417 - if (int(_candidates.size()) == _list_length) break;
1.418 + if (_curr_length == _list_length) break;
1.419 }
1.420 }
1.421 - if (int(_candidates.size()) < _list_length) {
1.422 - for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
1.423 + if (_curr_length < _list_length) {
1.424 + for (i = 0; i < _next_edge; ++i) {
1.425 + e = _edges[i];
1.426 if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.427 - _candidates.push_back(e);
1.428 + _candidates[_curr_length++] = e;
1.429 if (curr < min) {
1.430 min = curr;
1.431 - _ns._in_edge = e;
1.432 + _min_edge = i;
1.433 }
1.434 - if (int(_candidates.size()) == _list_length) break;
1.435 + if (_curr_length == _list_length) break;
1.436 }
1.437 }
1.438 }
1.439 - if (_candidates.size() == 0) return false;
1.440 + if (_curr_length == 0) return false;
1.441 _minor_count = 1;
1.442 - _next_edge = ++e;
1.443 + _ns._in_edge = _edges[_min_edge];
1.444 + _next_edge = i;
1.445 return true;
1.446 }
1.447 }; //class CandidateListPivotRule
1.448
1.449 + /// \brief Implementation of the "Altering Candidate List" pivot rule
1.450 + /// for the \ref NetworkSimplex "network simplex" algorithm.
1.451 + ///
1.452 + /// This class implements the "Altering Candidate List" pivot rule
1.453 + /// for the \ref NetworkSimplex "network simplex" algorithm.
1.454 + ///
1.455 + /// For more information see \ref NetworkSimplex::run().
1.456 + class AlteringListPivotRule
1.457 + {
1.458 + private:
1.459 +
1.460 + // References to the NetworkSimplex class
1.461 + NetworkSimplex &_ns;
1.462 + EdgeVector &_edges;
1.463 +
1.464 + EdgeVector _candidates;
1.465 + SCostMap _cand_cost;
1.466 + int _block_size, _head_length, _curr_length;
1.467 + int _next_edge;
1.468 +
1.469 + // Functor class to compare edges during sort of the candidate list
1.470 + class SortFunc
1.471 + {
1.472 + private:
1.473 + const SCostMap &_map;
1.474 + public:
1.475 + SortFunc(const SCostMap &map) : _map(map) {}
1.476 + bool operator()(const Edge &e1, const Edge &e2) {
1.477 + return _map[e1] < _map[e2];
1.478 + }
1.479 + };
1.480 +
1.481 + SortFunc _sort_func;
1.482 +
1.483 + public:
1.484 +
1.485 + /// Constructor
1.486 + AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.487 + _ns(ns), _edges(edges), _cand_cost(_ns._graph),
1.488 + _next_edge(0), _sort_func(_cand_cost)
1.489 + {
1.490 + // The main parameters of the pivot rule
1.491 + const double BLOCK_SIZE_FACTOR = 1.0;
1.492 + const int MIN_BLOCK_SIZE = 10;
1.493 + const double HEAD_LENGTH_FACTOR = 0.1;
1.494 + const int MIN_HEAD_LENGTH = 5;
1.495 +
1.496 + _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
1.497 + MIN_BLOCK_SIZE );
1.498 + _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
1.499 + MIN_HEAD_LENGTH );
1.500 + _candidates.resize(_head_length + _block_size);
1.501 + _curr_length = 0;
1.502 + }
1.503 +
1.504 + /// Find next entering edge
1.505 + inline bool findEnteringEdge() {
1.506 + // Checking the current candidate list
1.507 + Edge e;
1.508 + for (int idx = 0; idx < _curr_length; ++idx) {
1.509 + e = _candidates[idx];
1.510 + if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
1.511 + _candidates[idx--] = _candidates[--_curr_length];
1.512 + }
1.513 + }
1.514 +
1.515 + // Extending the list
1.516 + int cnt = _block_size;
1.517 + int last_edge = 0;
1.518 + int limit = _head_length;
1.519 + for (int i = _next_edge; i < int(_edges.size()); ++i) {
1.520 + e = _edges[i];
1.521 + if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.522 + _candidates[_curr_length++] = e;
1.523 + last_edge = i;
1.524 + }
1.525 + if (--cnt == 0) {
1.526 + if (_curr_length > limit) break;
1.527 + limit = 0;
1.528 + cnt = _block_size;
1.529 + }
1.530 + }
1.531 + if (_curr_length <= limit) {
1.532 + for (int i = 0; i < _next_edge; ++i) {
1.533 + e = _edges[i];
1.534 + if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.535 + _candidates[_curr_length++] = e;
1.536 + last_edge = i;
1.537 + }
1.538 + if (--cnt == 0) {
1.539 + if (_curr_length > limit) break;
1.540 + limit = 0;
1.541 + cnt = _block_size;
1.542 + }
1.543 + }
1.544 + }
1.545 + if (_curr_length == 0) return false;
1.546 + _next_edge = last_edge + 1;
1.547 +
1.548 + // Sorting the list partially
1.549 + EdgeVector::iterator sort_end = _candidates.begin();
1.550 + EdgeVector::iterator vector_end = _candidates.begin();
1.551 + for (int idx = 0; idx < _curr_length; ++idx) {
1.552 + ++vector_end;
1.553 + if (idx <= _head_length) ++sort_end;
1.554 + }
1.555 + partial_sort(_candidates.begin(), sort_end, vector_end, _sort_func);
1.556 +
1.557 + _ns._in_edge = _candidates[0];
1.558 + if (_curr_length > _head_length) {
1.559 + _candidates[0] = _candidates[_head_length - 1];
1.560 + _curr_length = _head_length - 1;
1.561 + } else {
1.562 + _candidates[0] = _candidates[_curr_length - 1];
1.563 + --_curr_length;
1.564 + }
1.565 +
1.566 + return true;
1.567 + }
1.568 + }; //class AlteringListPivotRule
1.569 +
1.570 private:
1.571
1.572 // State constants for edges
1.573 @@ -424,9 +527,6 @@
1.574 STATE_LOWER = 1
1.575 };
1.576
1.577 - // Constant for the combined pivot rule.
1.578 - static const int COMBINED_PIVOT_MAX_DEG = 5;
1.579 -
1.580 private:
1.581
1.582 // The directed graph the algorithm runs on
1.583 @@ -466,6 +566,9 @@
1.584 // The reduced cost map
1.585 ReducedCostMap _red_cost;
1.586
1.587 + // The non-artifical edges
1.588 + EdgeVector _edges;
1.589 +
1.590 // Members for handling the original graph
1.591 FlowMap *_flow_result;
1.592 PotentialMap *_potential_result;
1.593 @@ -672,9 +775,9 @@
1.594 if (_local_potential) delete _potential_result;
1.595 }
1.596
1.597 - /// \brief Sets the flow map.
1.598 + /// \brief Set the flow map.
1.599 ///
1.600 - /// Sets the flow map.
1.601 + /// Set the flow map.
1.602 ///
1.603 /// \return \c (*this)
1.604 NetworkSimplex& flowMap(FlowMap &map) {
1.605 @@ -686,9 +789,9 @@
1.606 return *this;
1.607 }
1.608
1.609 - /// \brief Sets the potential map.
1.610 + /// \brief Set the potential map.
1.611 ///
1.612 - /// Sets the potential map.
1.613 + /// Set the potential map.
1.614 ///
1.615 /// \return \c (*this)
1.616 NetworkSimplex& potentialMap(PotentialMap &map) {
1.617 @@ -701,8 +804,6 @@
1.618 }
1.619
1.620 /// \name Execution control
1.621 - /// The only way to execute the algorithm is to call the run()
1.622 - /// function.
1.623
1.624 /// @{
1.625
1.626 @@ -726,50 +827,49 @@
1.627 /// every iteration in a wraparound fashion and the best eligible
1.628 /// edge is selected from this block (\ref BlockSearchPivotRule).
1.629 ///
1.630 - /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
1.631 - /// examined in every iteration in a wraparound fashion and the best
1.632 - /// one is selected from them (\ref LimitedSearchPivotRule).
1.633 + /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
1.634 + /// built from eligible edges in a wraparound fashion and in the
1.635 + /// following minor iterations the best eligible edge is selected
1.636 + /// from this list (\ref CandidateListPivotRule).
1.637 ///
1.638 - /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
1.639 - /// built from eligible edges and it is used for edge selection in
1.640 - /// the following minor iterations (\ref CandidateListPivotRule).
1.641 + /// - ALTERING_LIST_PIVOT It is a modified version of the
1.642 + /// "Candidate List" pivot rule. It keeps only the several best
1.643 + /// eligible edges from the former candidate list and extends this
1.644 + /// list in every iteration (\ref AlteringListPivotRule).
1.645 ///
1.646 - /// - COMBINED_PIVOT This is a combined version of the two fastest
1.647 - /// pivot rules.
1.648 - /// For rather sparse graphs \ref LimitedSearchPivotRule
1.649 - /// "Limited Search" implementation is used, otherwise
1.650 - /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
1.651 - /// According to our benchmark tests this combined method is the
1.652 - /// most efficient.
1.653 + /// According to our comprehensive benchmark tests the "Block Search"
1.654 + /// pivot rule proved to be by far the fastest and the most robust
1.655 + /// on various test inputs. Thus it is the default option.
1.656 ///
1.657 /// \return \c true if a feasible flow can be found.
1.658 - bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
1.659 + bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
1.660 return init() && start(pivot_rule);
1.661 }
1.662
1.663 /// @}
1.664
1.665 /// \name Query Functions
1.666 - /// The result of the algorithm can be obtained using these
1.667 - /// functions.
1.668 - /// \n run() must be called before using them.
1.669 + /// The results of the algorithm can be obtained using these
1.670 + /// functions.\n
1.671 + /// \ref lemon::NetworkSimplex::run() "run()" must be called before
1.672 + /// using them.
1.673
1.674 /// @{
1.675
1.676 - /// \brief Returns a const reference to the edge map storing the
1.677 + /// \brief Return a const reference to the edge map storing the
1.678 /// found flow.
1.679 ///
1.680 - /// Returns a const reference to the edge map storing the found flow.
1.681 + /// Return a const reference to the edge map storing the found flow.
1.682 ///
1.683 /// \pre \ref run() must be called before using this function.
1.684 const FlowMap& flowMap() const {
1.685 return *_flow_result;
1.686 }
1.687
1.688 - /// \brief Returns a const reference to the node map storing the
1.689 + /// \brief Return a const reference to the node map storing the
1.690 /// found potentials (the dual solution).
1.691 ///
1.692 - /// Returns a const reference to the node map storing the found
1.693 + /// Return a const reference to the node map storing the found
1.694 /// potentials (the dual solution).
1.695 ///
1.696 /// \pre \ref run() must be called before using this function.
1.697 @@ -777,27 +877,27 @@
1.698 return *_potential_result;
1.699 }
1.700
1.701 - /// \brief Returns the flow on the given edge.
1.702 + /// \brief Return the flow on the given edge.
1.703 ///
1.704 - /// Returns the flow on the given edge.
1.705 + /// Return the flow on the given edge.
1.706 ///
1.707 /// \pre \ref run() must be called before using this function.
1.708 Capacity flow(const typename Graph::Edge& edge) const {
1.709 return (*_flow_result)[edge];
1.710 }
1.711
1.712 - /// \brief Returns the potential of the given node.
1.713 + /// \brief Return the potential of the given node.
1.714 ///
1.715 - /// Returns the potential of the given node.
1.716 + /// Return the potential of the given node.
1.717 ///
1.718 /// \pre \ref run() must be called before using this function.
1.719 Cost potential(const typename Graph::Node& node) const {
1.720 return (*_potential_result)[node];
1.721 }
1.722
1.723 - /// \brief Returns the total cost of the found flow.
1.724 + /// \brief Return the total cost of the found flow.
1.725 ///
1.726 - /// Returns the total cost of the found flow. The complexity of the
1.727 + /// Return the total cost of the found flow. The complexity of the
1.728 /// function is \f$ O(e) \f$.
1.729 ///
1.730 /// \pre \ref run() must be called before using this function.
1.731 @@ -812,7 +912,7 @@
1.732
1.733 private:
1.734
1.735 - /// \brief Extends the underlaying graph and initializes all the
1.736 + /// \brief Extend the underlying graph and initialize all the
1.737 /// node and edge maps.
1.738 bool init() {
1.739 if (!_valid_supply) return false;
1.740 @@ -827,8 +927,10 @@
1.741 _local_potential = true;
1.742 }
1.743
1.744 - // Initializing state and flow maps
1.745 + // Initializing the edge vector and the edge maps
1.746 + _edges.reserve(countEdges(_graph));
1.747 for (EdgeIt e(_graph); e != INVALID; ++e) {
1.748 + _edges.push_back(e);
1.749 _flow[e] = 0;
1.750 _state[e] = STATE_LOWER;
1.751 }
1.752 @@ -873,8 +975,8 @@
1.753 return true;
1.754 }
1.755
1.756 - /// Finds the join node.
1.757 - Node findJoinNode() {
1.758 + /// Find the join node.
1.759 + inline Node findJoinNode() {
1.760 Node u = _graph.source(_in_edge);
1.761 Node v = _graph.target(_in_edge);
1.762 while (u != v) {
1.763 @@ -888,9 +990,10 @@
1.764 return u;
1.765 }
1.766
1.767 - /// \brief Finds the leaving edge of the cycle. Returns \c true if
1.768 - /// the leaving edge is not the same as the entering edge.
1.769 - bool findLeavingEdge() {
1.770 + /// \brief Find the leaving edge of the cycle.
1.771 + /// \return \c true if the leaving edge is not the same as the
1.772 + /// entering edge.
1.773 + inline bool findLeavingEdge() {
1.774 // Initializing first and second nodes according to the direction
1.775 // of the cycle
1.776 if (_state[_in_edge] == STATE_LOWER) {
1.777 @@ -934,8 +1037,8 @@
1.778 return result;
1.779 }
1.780
1.781 - /// Changes \c flow and \c state edge maps.
1.782 - void changeFlows(bool change) {
1.783 + /// Change \c flow and \c state edge maps.
1.784 + inline void changeFlows(bool change) {
1.785 // Augmenting along the cycle
1.786 if (delta > 0) {
1.787 Capacity val = _state[_in_edge] * delta;
1.788 @@ -957,8 +1060,8 @@
1.789 }
1.790 }
1.791
1.792 - /// Updates \c thread and \c parent node maps.
1.793 - void updateThreadParent() {
1.794 + /// Update \c thread and \c parent node maps.
1.795 + inline void updateThreadParent() {
1.796 Node u;
1.797 v_out = _parent[u_out];
1.798
1.799 @@ -1016,8 +1119,8 @@
1.800 }
1.801 }
1.802
1.803 - /// Updates \c pred_edge and \c forward node maps.
1.804 - void updatePredEdge() {
1.805 + /// Update \c pred_edge and \c forward node maps.
1.806 + inline void updatePredEdge() {
1.807 Node u = u_out, v;
1.808 while (u != u_in) {
1.809 v = _parent[u];
1.810 @@ -1029,8 +1132,8 @@
1.811 _forward[u_in] = (u_in == _graph.source(_in_edge));
1.812 }
1.813
1.814 - /// Updates \c depth and \c potential node maps.
1.815 - void updateDepthPotential() {
1.816 + /// Update \c depth and \c potential node maps.
1.817 + inline void updateDepthPotential() {
1.818 _depth[u_in] = _depth[v_in] + 1;
1.819 _potential[u_in] = _forward[u_in] ?
1.820 _potential[v_in] - _cost[_pred_edge[u_in]] :
1.821 @@ -1049,8 +1152,9 @@
1.822 }
1.823 }
1.824
1.825 - /// Executes the algorithm.
1.826 + /// Execute the algorithm.
1.827 bool start(PivotRuleEnum pivot_rule) {
1.828 + // Selecting the pivot rule implementation
1.829 switch (pivot_rule) {
1.830 case FIRST_ELIGIBLE_PIVOT:
1.831 return start<FirstEligiblePivotRule>();
1.832 @@ -1058,23 +1162,17 @@
1.833 return start<BestEligiblePivotRule>();
1.834 case BLOCK_SEARCH_PIVOT:
1.835 return start<BlockSearchPivotRule>();
1.836 - case LIMITED_SEARCH_PIVOT:
1.837 - return start<LimitedSearchPivotRule>();
1.838 case CANDIDATE_LIST_PIVOT:
1.839 return start<CandidateListPivotRule>();
1.840 - case COMBINED_PIVOT:
1.841 - if ( countEdges(_graph) / countNodes(_graph) <=
1.842 - COMBINED_PIVOT_MAX_DEG )
1.843 - return start<LimitedSearchPivotRule>();
1.844 - else
1.845 - return start<BlockSearchPivotRule>();
1.846 + case ALTERING_LIST_PIVOT:
1.847 + return start<AlteringListPivotRule>();
1.848 }
1.849 return false;
1.850 }
1.851
1.852 template<class PivotRuleImplementation>
1.853 bool start() {
1.854 - PivotRuleImplementation pivot(*this);
1.855 + PivotRuleImplementation pivot(*this, _edges);
1.856
1.857 // Executing the network simplex algorithm
1.858 while (pivot.findEnteringEdge()) {