Major improvements in NetworkSimplex.
Main changes:
- Use -potenital[] instead of potential[] to conform to the usual
terminology.
- Use function parameter instead of #define commands to select pivot rule.
- Use much faster implementation for the candidate list pivot rule.
It is about 5-20 times faster now.
- Add a new pivot rule called "Limited Search" that is a modified
version of "Block Search". It is about 25 percent faster on rather
sparse graphs.
- By default "Limited Search" is used for sparse graphs and
"Block Search" is used otherwise. This combined method is the most
efficient on every input class.
- Change the name of private members to start with "_".
- Change the name of function parameters not to start with "_".
- Remove unnecessary documentation for private members.
- Many doc improvements.
1.1 --- a/lemon/network_simplex.h Mon Feb 18 03:30:53 2008 +0000
1.2 +++ b/lemon/network_simplex.h Mon Feb 18 03:32:06 2008 +0000
1.3 @@ -22,47 +22,15 @@
1.4 /// \ingroup min_cost_flow
1.5 ///
1.6 /// \file
1.7 -/// \brief The 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 +
1.13 #include <lemon/graph_adaptor.h>
1.14 #include <lemon/graph_utils.h>
1.15 #include <lemon/smart_graph.h>
1.16 -
1.17 -/// \brief The pivot rule used in the algorithm.
1.18 -//#define FIRST_ELIGIBLE_PIVOT
1.19 -//#define BEST_ELIGIBLE_PIVOT
1.20 -#define EDGE_BLOCK_PIVOT
1.21 -//#define CANDIDATE_LIST_PIVOT
1.22 -//#define SORTED_LIST_PIVOT
1.23 -
1.24 -//#define _DEBUG_ITER_
1.25 -
1.26 -
1.27 -// State constant for edges at their lower bounds.
1.28 -#define LOWER 1
1.29 -// State constant for edges in the spanning tree.
1.30 -#define TREE 0
1.31 -// State constant for edges at their upper bounds.
1.32 -#define UPPER -1
1.33 -
1.34 -#ifdef EDGE_BLOCK_PIVOT
1.35 - #include <lemon/math.h>
1.36 - #define MIN_BLOCK_SIZE 10
1.37 -#endif
1.38 -
1.39 -#ifdef CANDIDATE_LIST_PIVOT
1.40 - #include <vector>
1.41 - #define LIST_LENGTH_DIV 50
1.42 - #define MINOR_LIMIT_DIV 200
1.43 -#endif
1.44 -
1.45 -#ifdef SORTED_LIST_PIVOT
1.46 - #include <vector>
1.47 - #include <algorithm>
1.48 - #define LIST_LENGTH_DIV 100
1.49 - #define LOWER_DIV 4
1.50 -#endif
1.51 +#include <lemon/math.h>
1.52
1.53 namespace lemon {
1.54
1.55 @@ -75,31 +43,38 @@
1.56 /// \ref NetworkSimplex implements the network simplex algorithm for
1.57 /// finding a minimum cost flow.
1.58 ///
1.59 - /// \param Graph The directed graph type the algorithm runs on.
1.60 - /// \param LowerMap The type of the lower bound map.
1.61 - /// \param CapacityMap The type of the capacity (upper bound) map.
1.62 - /// \param CostMap The type of the cost (length) map.
1.63 - /// \param SupplyMap The type of the supply map.
1.64 + /// \tparam Graph The directed graph type the algorithm runs on.
1.65 + /// \tparam LowerMap The type of the lower bound map.
1.66 + /// \tparam CapacityMap The type of the capacity (upper bound) map.
1.67 + /// \tparam CostMap The type of the cost (length) map.
1.68 + /// \tparam SupplyMap The type of the supply map.
1.69 ///
1.70 /// \warning
1.71 - /// - Edge capacities and costs should be non-negative integers.
1.72 - /// However \c CostMap::Value should be signed type.
1.73 - /// - Supply values should be signed integers.
1.74 - /// - \c LowerMap::Value must be convertible to
1.75 - /// \c CapacityMap::Value and \c CapacityMap::Value must be
1.76 - /// convertible to \c SupplyMap::Value.
1.77 + /// - Edge capacities and costs should be \e non-negative \e integers.
1.78 + /// - Supply values should be \e signed \e integers.
1.79 + /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
1.80 + /// - \c CapacityMap::Value and \c SupplyMap::Value must be
1.81 + /// convertible to each other.
1.82 + /// - All value types must be convertible to \c CostMap::Value, which
1.83 + /// must be signed type.
1.84 + ///
1.85 + /// \note \ref NetworkSimplex provides six different pivot rule
1.86 + /// implementations that significantly affect the efficiency of the
1.87 + /// algorithm.
1.88 + /// By default a combined pivot rule is used, which is the fastest
1.89 + /// implementation according to our benchmark tests.
1.90 + /// Another pivot rule can be selected using \ref run() function
1.91 + /// with the proper parameter.
1.92 ///
1.93 /// \author Peter Kovacs
1.94
1.95 template < typename Graph,
1.96 typename LowerMap = typename Graph::template EdgeMap<int>,
1.97 - typename CapacityMap = LowerMap,
1.98 + typename CapacityMap = typename Graph::template EdgeMap<int>,
1.99 typename CostMap = typename Graph::template EdgeMap<int>,
1.100 - typename SupplyMap = typename Graph::template NodeMap
1.101 - <typename CapacityMap::Value> >
1.102 + typename SupplyMap = typename Graph::template NodeMap<int> >
1.103 class NetworkSimplex
1.104 {
1.105 - typedef typename LowerMap::Value Lower;
1.106 typedef typename CapacityMap::Value Capacity;
1.107 typedef typename CostMap::Value Cost;
1.108 typedef typename SupplyMap::Value Supply;
1.109 @@ -107,7 +82,6 @@
1.110 typedef SmartGraph SGraph;
1.111 GRAPH_TYPEDEFS(typename SGraph);
1.112
1.113 - typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
1.114 typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
1.115 typedef typename SGraph::template EdgeMap<Cost> SCostMap;
1.116 typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
1.117 @@ -129,77 +103,380 @@
1.118 /// The type of the potential map.
1.119 typedef typename Graph::template NodeMap<Cost> PotentialMap;
1.120
1.121 - protected:
1.122 + public:
1.123
1.124 + /// Enum type to select the pivot rule used by \ref run().
1.125 + enum PivotRuleEnum {
1.126 + FIRST_ELIGIBLE_PIVOT,
1.127 + BEST_ELIGIBLE_PIVOT,
1.128 + BLOCK_SEARCH_PIVOT,
1.129 + LIMITED_SEARCH_PIVOT,
1.130 + CANDIDATE_LIST_PIVOT,
1.131 + COMBINED_PIVOT
1.132 + };
1.133 +
1.134 + private:
1.135 +
1.136 + /// \brief Map adaptor class for handling reduced edge costs.
1.137 + ///
1.138 /// Map adaptor class for handling reduced edge costs.
1.139 class ReducedCostMap : public MapBase<Edge, Cost>
1.140 {
1.141 private:
1.142
1.143 - const SGraph &gr;
1.144 - const SCostMap &cost_map;
1.145 - const SPotentialMap &pot_map;
1.146 + const SGraph &_gr;
1.147 + const SCostMap &_cost_map;
1.148 + const SPotentialMap &_pot_map;
1.149
1.150 public:
1.151
1.152 - ReducedCostMap( const SGraph &_gr,
1.153 - const SCostMap &_cm,
1.154 - const SPotentialMap &_pm ) :
1.155 - gr(_gr), cost_map(_cm), pot_map(_pm) {}
1.156 + ///\e
1.157 + ReducedCostMap( const SGraph &gr,
1.158 + const SCostMap &cost_map,
1.159 + const SPotentialMap &pot_map ) :
1.160 + _gr(gr), _cost_map(cost_map), _pot_map(pm) {}
1.161
1.162 + ///\e
1.163 Cost operator[](const Edge &e) const {
1.164 - return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
1.165 + return _cost_map[e] + _pot_map[_gr.source(e)]
1.166 + - _pot_map[_gr.target(e)];
1.167 }
1.168
1.169 }; //class ReducedCostMap
1.170
1.171 - protected:
1.172 + private:
1.173
1.174 - /// The directed graph the algorithm runs on.
1.175 - SGraph graph;
1.176 - /// The original graph.
1.177 - const Graph &graph_ref;
1.178 - /// The original lower bound map.
1.179 - const LowerMap *lower;
1.180 - /// The capacity map.
1.181 - SCapacityMap capacity;
1.182 - /// The cost map.
1.183 - SCostMap cost;
1.184 - /// The supply map.
1.185 - SSupplyMap supply;
1.186 - /// The reduced cost map.
1.187 - ReducedCostMap red_cost;
1.188 - bool valid_supply;
1.189 + /// \brief Implementation of the "First Eligible" pivot rule for the
1.190 + /// \ref NetworkSimplex "network simplex" algorithm.
1.191 + ///
1.192 + /// This class implements the "First Eligible" pivot rule
1.193 + /// for the \ref NetworkSimplex "network simplex" algorithm.
1.194 + class FirstEligiblePivotRule
1.195 + {
1.196 + private:
1.197
1.198 - /// The edge map of the current flow.
1.199 - SCapacityMap flow;
1.200 - /// The potential node map.
1.201 - SPotentialMap potential;
1.202 - FlowMap flow_result;
1.203 - PotentialMap potential_result;
1.204 + NetworkSimplex &_ns;
1.205 + EdgeIt _next_edge;
1.206
1.207 - /// Node reference for the original graph.
1.208 - NodeRefMap node_ref;
1.209 - /// Edge reference for the original graph.
1.210 - EdgeRefMap edge_ref;
1.211 + public:
1.212
1.213 - /// The \c depth node map of the spanning tree structure.
1.214 - IntNodeMap depth;
1.215 - /// The \c parent node map of the spanning tree structure.
1.216 - NodeNodeMap parent;
1.217 - /// The \c pred_edge node map of the spanning tree structure.
1.218 - EdgeNodeMap pred_edge;
1.219 - /// The \c thread node map of the spanning tree structure.
1.220 - NodeNodeMap thread;
1.221 - /// The \c forward node map of the spanning tree structure.
1.222 - BoolNodeMap forward;
1.223 - /// The state edge map.
1.224 - IntEdgeMap state;
1.225 - /// The root node of the starting spanning tree.
1.226 - Node root;
1.227 + /// Constructor.
1.228 + FirstEligiblePivotRule(NetworkSimplex &ns) :
1.229 + _ns(ns), _next_edge(ns._graph) {}
1.230 +
1.231 + /// Finds the next entering edge.
1.232 + bool findEnteringEdge() {
1.233 + for (EdgeIt e = _next_edge; e != INVALID; ++e) {
1.234 + if (_ns._state[e] * _ns._red_cost[e] < 0) {
1.235 + _ns._in_edge = e;
1.236 + _next_edge = ++e;
1.237 + return true;
1.238 + }
1.239 + }
1.240 + for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
1.241 + if (_ns._state[e] * _ns._red_cost[e] < 0) {
1.242 + _ns._in_edge = e;
1.243 + _next_edge = ++e;
1.244 + return true;
1.245 + }
1.246 + }
1.247 + return false;
1.248 + }
1.249 + }; //class FirstEligiblePivotRule
1.250 +
1.251 + /// \brief Implementation of the "Best Eligible" pivot rule for the
1.252 + /// \ref NetworkSimplex "network simplex" algorithm.
1.253 + ///
1.254 + /// This class implements the "Best Eligible" pivot rule
1.255 + /// for the \ref NetworkSimplex "network simplex" algorithm.
1.256 + class BestEligiblePivotRule
1.257 + {
1.258 + private:
1.259 +
1.260 + NetworkSimplex &_ns;
1.261 +
1.262 + public:
1.263 +
1.264 + /// Constructor.
1.265 + BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
1.266 +
1.267 + /// Finds the next entering edge.
1.268 + bool findEnteringEdge() {
1.269 + Cost min = 0;
1.270 + for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
1.271 + if (_ns._state[e] * _ns._red_cost[e] < min) {
1.272 + min = _ns._state[e] * _ns._red_cost[e];
1.273 + _ns._in_edge = e;
1.274 + }
1.275 + }
1.276 + return min < 0;
1.277 + }
1.278 + }; //class BestEligiblePivotRule
1.279 +
1.280 + /// \brief Implementation of the "Block Search" pivot rule for the
1.281 + /// \ref NetworkSimplex "network simplex" algorithm.
1.282 + ///
1.283 + /// This class implements the "Block Search" pivot rule
1.284 + /// for the \ref NetworkSimplex "network simplex" algorithm.
1.285 + class BlockSearchPivotRule
1.286 + {
1.287 + private:
1.288 +
1.289 + NetworkSimplex &_ns;
1.290 + EdgeIt _next_edge, _min_edge;
1.291 + int _block_size;
1.292 +
1.293 + static const int MIN_BLOCK_SIZE = 10;
1.294 +
1.295 + public:
1.296 +
1.297 + /// Constructor.
1.298 + BlockSearchPivotRule(NetworkSimplex &ns) :
1.299 + _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
1.300 + {
1.301 + _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
1.302 + if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
1.303 + }
1.304 +
1.305 + /// Finds the next entering edge.
1.306 + bool findEnteringEdge() {
1.307 + Cost curr, min = 0;
1.308 + int cnt = 0;
1.309 + for (EdgeIt e = _next_edge; e != INVALID; ++e) {
1.310 + if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.311 + min = curr;
1.312 + _min_edge = e;
1.313 + }
1.314 + if (++cnt == _block_size) {
1.315 + if (min < 0) break;
1.316 + cnt = 0;
1.317 + }
1.318 + }
1.319 + if (min == 0) {
1.320 + for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
1.321 + if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.322 + min = curr;
1.323 + _min_edge = e;
1.324 + }
1.325 + if (++cnt == _block_size) {
1.326 + if (min < 0) break;
1.327 + cnt = 0;
1.328 + }
1.329 + }
1.330 + }
1.331 + _ns._in_edge = _min_edge;
1.332 + _next_edge = ++_min_edge;
1.333 + return min < 0;
1.334 + }
1.335 + }; //class BlockSearchPivotRule
1.336 +
1.337 + /// \brief Implementation of the "Limited Search" pivot rule for the
1.338 + /// \ref NetworkSimplex "network simplex" algorithm.
1.339 + ///
1.340 + /// This class implements the "Limited Search" pivot rule
1.341 + /// for the \ref NetworkSimplex "network simplex" algorithm.
1.342 + class LimitedSearchPivotRule
1.343 + {
1.344 + private:
1.345 +
1.346 + NetworkSimplex &_ns;
1.347 + EdgeIt _next_edge, _min_edge;
1.348 + int _sample_size;
1.349 +
1.350 + static const int MIN_SAMPLE_SIZE = 10;
1.351 + static const double SAMPLE_SIZE_FACTOR = 0.0015;
1.352 +
1.353 + public:
1.354 +
1.355 + /// Constructor.
1.356 + LimitedSearchPivotRule(NetworkSimplex &ns) :
1.357 + _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
1.358 + {
1.359 + _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph));
1.360 + if (_sample_size < MIN_SAMPLE_SIZE)
1.361 + _sample_size = MIN_SAMPLE_SIZE;
1.362 + }
1.363 +
1.364 + /// Finds the next entering edge.
1.365 + bool findEnteringEdge() {
1.366 + Cost curr, min = 0;
1.367 + int cnt = 0;
1.368 + for (EdgeIt e = _next_edge; e != INVALID; ++e) {
1.369 + if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.370 + min = curr;
1.371 + _min_edge = e;
1.372 + }
1.373 + if (curr < 0 && ++cnt == _sample_size) break;
1.374 + }
1.375 + if (min == 0) {
1.376 + for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
1.377 + if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.378 + min = curr;
1.379 + _min_edge = e;
1.380 + }
1.381 + if (curr < 0 && ++cnt == _sample_size) break;
1.382 + }
1.383 + }
1.384 + _ns._in_edge = _min_edge;
1.385 + _next_edge = ++_min_edge;
1.386 + return min < 0;
1.387 + }
1.388 + }; //class LimitedSearchPivotRule
1.389 +
1.390 + /// \brief Implementation of the "Candidate List" pivot rule for the
1.391 + /// \ref NetworkSimplex "network simplex" algorithm.
1.392 + ///
1.393 + /// This class implements the "Candidate List" pivot rule
1.394 + /// for the \ref NetworkSimplex "network simplex" algorithm.
1.395 + class CandidateListPivotRule
1.396 + {
1.397 + private:
1.398 +
1.399 + NetworkSimplex &_ns;
1.400 +
1.401 + // The list of candidate edges.
1.402 + std::vector<Edge> _candidates;
1.403 + // The maximum length of the edge list.
1.404 + int _list_length;
1.405 + // The maximum number of minor iterations between two major
1.406 + // itarations.
1.407 + int _minor_limit;
1.408 +
1.409 + int _minor_count;
1.410 + EdgeIt _next_edge;
1.411 +
1.412 + static const double LIST_LENGTH_FACTOR = 0.002;
1.413 + static const double MINOR_LIMIT_FACTOR = 0.1;
1.414 + static const int MIN_LIST_LENGTH = 10;
1.415 + static const int MIN_MINOR_LIMIT = 2;
1.416 +
1.417 + public:
1.418 +
1.419 + /// Constructor.
1.420 + CandidateListPivotRule(NetworkSimplex &ns) :
1.421 + _ns(ns), _next_edge(ns._graph)
1.422 + {
1.423 + int edge_num = countEdges(_ns._graph);
1.424 + _minor_count = 0;
1.425 + _list_length = int(edge_num * LIST_LENGTH_FACTOR);
1.426 + if (_list_length < MIN_LIST_LENGTH)
1.427 + _list_length = MIN_LIST_LENGTH;
1.428 + _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR);
1.429 + if (_minor_limit < MIN_MINOR_LIMIT)
1.430 + _minor_limit = MIN_MINOR_LIMIT;
1.431 + }
1.432 +
1.433 + /// Finds the next entering edge.
1.434 + bool findEnteringEdge() {
1.435 + Cost min, curr;
1.436 + if (_minor_count < _minor_limit && _candidates.size() > 0) {
1.437 + // Minor iteration
1.438 + ++_minor_count;
1.439 + Edge e;
1.440 + min = 0;
1.441 + for (int i = 0; i < int(_candidates.size()); ++i) {
1.442 + e = _candidates[i];
1.443 + if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.444 + min = curr;
1.445 + _ns._in_edge = e;
1.446 + }
1.447 + }
1.448 + if (min < 0) return true;
1.449 + }
1.450 +
1.451 + // Major iteration
1.452 + _candidates.clear();
1.453 + EdgeIt e = _next_edge;
1.454 + min = 0;
1.455 + for ( ; e != INVALID; ++e) {
1.456 + if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.457 + _candidates.push_back(e);
1.458 + if (curr < min) {
1.459 + min = curr;
1.460 + _ns._in_edge = e;
1.461 + }
1.462 + if (int(_candidates.size()) == _list_length) break;
1.463 + }
1.464 + }
1.465 + if (int(_candidates.size()) < _list_length) {
1.466 + for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
1.467 + if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.468 + _candidates.push_back(e);
1.469 + if (curr < min) {
1.470 + min = curr;
1.471 + _ns._in_edge = e;
1.472 + }
1.473 + if (int(_candidates.size()) == _list_length) break;
1.474 + }
1.475 + }
1.476 + }
1.477 + if (_candidates.size() == 0) return false;
1.478 + _minor_count = 1;
1.479 + _next_edge = ++e;
1.480 + return true;
1.481 + }
1.482 + }; //class CandidateListPivotRule
1.483 +
1.484 + private:
1.485 +
1.486 + // State constant for edges at their lower bounds.
1.487 + static const int STATE_LOWER = 1;
1.488 + // State constant for edges in the spanning tree.
1.489 + static const int STATE_TREE = 0;
1.490 + // State constant for edges at their upper bounds.
1.491 + static const int STATE_UPPER = -1;
1.492 +
1.493 + // Constant for the combined pivot rule.
1.494 + static const int COMBINED_PIVOT_MAX_DEG = 5;
1.495 +
1.496 + private:
1.497 +
1.498 + // The directed graph the algorithm runs on
1.499 + SGraph _graph;
1.500 + // The original graph
1.501 + const Graph &_graph_ref;
1.502 + // The original lower bound map
1.503 + const LowerMap *_lower;
1.504 + // The capacity map
1.505 + SCapacityMap _capacity;
1.506 + // The cost map
1.507 + SCostMap _cost;
1.508 + // The supply map
1.509 + SSupplyMap _supply;
1.510 + bool _valid_supply;
1.511 +
1.512 + // Edge map of the current flow
1.513 + SCapacityMap _flow;
1.514 + // Node map of the current potentials
1.515 + SPotentialMap _potential;
1.516 +
1.517 + // The depth node map of the spanning tree structure
1.518 + IntNodeMap _depth;
1.519 + // The parent node map of the spanning tree structure
1.520 + NodeNodeMap _parent;
1.521 + // The pred_edge node map of the spanning tree structure
1.522 + EdgeNodeMap _pred_edge;
1.523 + // The thread node map of the spanning tree structure
1.524 + NodeNodeMap _thread;
1.525 + // The forward node map of the spanning tree structure
1.526 + BoolNodeMap _forward;
1.527 + // The state edge map
1.528 + IntEdgeMap _state;
1.529 + // The root node of the starting spanning tree
1.530 + Node _root;
1.531 +
1.532 + // The reduced cost map
1.533 + ReducedCostMap _red_cost;
1.534 +
1.535 + // Members for handling the original graph
1.536 + FlowMap _flow_result;
1.537 + PotentialMap _potential_result;
1.538 + NodeRefMap _node_ref;
1.539 + EdgeRefMap _edge_ref;
1.540
1.541 // The entering edge of the current pivot iteration.
1.542 - Edge in_edge;
1.543 + Edge _in_edge;
1.544 +
1.545 // Temporary nodes used in the current pivot iteration.
1.546 Node join, u_in, v_in, u_out, v_out;
1.547 Node right, first, second, last;
1.548 @@ -208,82 +485,56 @@
1.549 // pivot iteration.
1.550 Capacity delta;
1.551
1.552 -#ifdef EDGE_BLOCK_PIVOT
1.553 - /// The size of blocks for the "Edge Block" pivot rule.
1.554 - int block_size;
1.555 -#endif
1.556 -#ifdef CANDIDATE_LIST_PIVOT
1.557 - /// \brief The list of candidate edges for the "Candidate List"
1.558 - /// pivot rule.
1.559 - std::vector<Edge> candidates;
1.560 - /// \brief The maximum length of the edge list for the
1.561 - /// "Candidate List" pivot rule.
1.562 - int list_length;
1.563 - /// \brief The maximum number of minor iterations between two major
1.564 - /// itarations.
1.565 - int minor_limit;
1.566 - /// \brief The number of minor iterations.
1.567 - int minor_count;
1.568 -#endif
1.569 -#ifdef SORTED_LIST_PIVOT
1.570 - /// \brief The list of candidate edges for the "Sorted List"
1.571 - /// pivot rule.
1.572 - std::vector<Edge> candidates;
1.573 - /// \brief The maximum length of the edge list for the
1.574 - /// "Sorted List" pivot rule.
1.575 - int list_length;
1.576 - int list_index;
1.577 -#endif
1.578 -
1.579 public :
1.580
1.581 /// \brief General constructor of the class (with lower bounds).
1.582 ///
1.583 /// General constructor of the class (with lower bounds).
1.584 ///
1.585 - /// \param _graph The directed graph the algorithm runs on.
1.586 - /// \param _lower The lower bounds of the edges.
1.587 - /// \param _capacity The capacities (upper bounds) of the edges.
1.588 - /// \param _cost The cost (length) values of the edges.
1.589 - /// \param _supply The supply values of the nodes (signed).
1.590 - NetworkSimplex( const Graph &_graph,
1.591 - const LowerMap &_lower,
1.592 - const CapacityMap &_capacity,
1.593 - const CostMap &_cost,
1.594 - const SupplyMap &_supply ) :
1.595 - graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
1.596 - supply(graph), flow(graph), flow_result(_graph), potential(graph),
1.597 - potential_result(_graph), depth(graph), parent(graph),
1.598 - pred_edge(graph), thread(graph), forward(graph), state(graph),
1.599 - node_ref(graph_ref), edge_ref(graph_ref),
1.600 - red_cost(graph, cost, potential)
1.601 + /// \param graph The directed graph the algorithm runs on.
1.602 + /// \param lower The lower bounds of the edges.
1.603 + /// \param capacity The capacities (upper bounds) of the edges.
1.604 + /// \param cost The cost (length) values of the edges.
1.605 + /// \param supply The supply values of the nodes (signed).
1.606 + NetworkSimplex( const Graph &graph,
1.607 + const LowerMap &lower,
1.608 + const CapacityMap &capacity,
1.609 + const CostMap &cost,
1.610 + const SupplyMap &supply ) :
1.611 + _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
1.612 + _cost(_graph), _supply(_graph), _flow(_graph),
1.613 + _potential(_graph), _depth(_graph), _parent(_graph),
1.614 + _pred_edge(_graph), _thread(_graph), _forward(_graph),
1.615 + _state(_graph), _red_cost(_graph, _cost, _potential),
1.616 + _flow_result(graph), _potential_result(graph),
1.617 + _node_ref(graph), _edge_ref(graph)
1.618 {
1.619 // Checking the sum of supply values
1.620 Supply sum = 0;
1.621 - for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
1.622 - sum += _supply[n];
1.623 - if (!(valid_supply = sum == 0)) return;
1.624 + for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1.625 + sum += supply[n];
1.626 + if (!(_valid_supply = sum == 0)) return;
1.627
1.628 - // Copying graph_ref to graph
1.629 - graph.reserveNode(countNodes(graph_ref) + 1);
1.630 - graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));
1.631 - copyGraph(graph, graph_ref)
1.632 - .edgeMap(cost, _cost)
1.633 - .nodeRef(node_ref)
1.634 - .edgeRef(edge_ref)
1.635 + // Copying _graph_ref to _graph
1.636 + _graph.reserveNode(countNodes(_graph_ref) + 1);
1.637 + _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
1.638 + copyGraph(_graph, _graph_ref)
1.639 + .edgeMap(_cost, cost)
1.640 + .nodeRef(_node_ref)
1.641 + .edgeRef(_edge_ref)
1.642 .run();
1.643
1.644 // Removing non-zero lower bounds
1.645 - for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
1.646 - capacity[edge_ref[e]] = _capacity[e] - _lower[e];
1.647 + for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
1.648 + _capacity[_edge_ref[e]] = capacity[e] - lower[e];
1.649 }
1.650 - for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
1.651 - Supply s = _supply[n];
1.652 - for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.653 - s += _lower[e];
1.654 - for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.655 - s -= _lower[e];
1.656 - supply[node_ref[n]] = s;
1.657 + for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
1.658 + Supply s = supply[n];
1.659 + for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
1.660 + s += lower[e];
1.661 + for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
1.662 + s -= lower[e];
1.663 + _supply[_node_ref[n]] = s;
1.664 }
1.665 }
1.666
1.667 @@ -291,34 +542,35 @@
1.668 ///
1.669 /// General constructor of the class (without lower bounds).
1.670 ///
1.671 - /// \param _graph The directed graph the algorithm runs on.
1.672 - /// \param _capacity The capacities (upper bounds) of the edges.
1.673 - /// \param _cost The cost (length) values of the edges.
1.674 - /// \param _supply The supply values of the nodes (signed).
1.675 - NetworkSimplex( const Graph &_graph,
1.676 - const CapacityMap &_capacity,
1.677 - const CostMap &_cost,
1.678 - const SupplyMap &_supply ) :
1.679 - graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
1.680 - supply(graph), flow(graph), flow_result(_graph), potential(graph),
1.681 - potential_result(_graph), depth(graph), parent(graph),
1.682 - pred_edge(graph), thread(graph), forward(graph), state(graph),
1.683 - node_ref(graph_ref), edge_ref(graph_ref),
1.684 - red_cost(graph, cost, potential)
1.685 + /// \param graph The directed graph the algorithm runs on.
1.686 + /// \param capacity The capacities (upper bounds) of the edges.
1.687 + /// \param cost The cost (length) values of the edges.
1.688 + /// \param supply The supply values of the nodes (signed).
1.689 + NetworkSimplex( const Graph &graph,
1.690 + const CapacityMap &capacity,
1.691 + const CostMap &cost,
1.692 + const SupplyMap &supply ) :
1.693 + _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
1.694 + _cost(_graph), _supply(_graph), _flow(_graph),
1.695 + _potential(_graph), _depth(_graph), _parent(_graph),
1.696 + _pred_edge(_graph), _thread(_graph), _forward(_graph),
1.697 + _state(_graph), _red_cost(_graph, _cost, _potential),
1.698 + _flow_result(graph), _potential_result(graph),
1.699 + _node_ref(graph), _edge_ref(graph)
1.700 {
1.701 // Checking the sum of supply values
1.702 Supply sum = 0;
1.703 - for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
1.704 - sum += _supply[n];
1.705 - if (!(valid_supply = sum == 0)) return;
1.706 + for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1.707 + sum += supply[n];
1.708 + if (!(_valid_supply = sum == 0)) return;
1.709
1.710 - // Copying graph_ref to graph
1.711 - copyGraph(graph, graph_ref)
1.712 - .edgeMap(capacity, _capacity)
1.713 - .edgeMap(cost, _cost)
1.714 - .nodeMap(supply, _supply)
1.715 - .nodeRef(node_ref)
1.716 - .edgeRef(edge_ref)
1.717 + // Copying _graph_ref to graph
1.718 + copyGraph(_graph, _graph_ref)
1.719 + .edgeMap(_capacity, capacity)
1.720 + .edgeMap(_cost, cost)
1.721 + .nodeMap(_supply, supply)
1.722 + .nodeRef(_node_ref)
1.723 + .edgeRef(_edge_ref)
1.724 .run();
1.725 }
1.726
1.727 @@ -326,115 +578,150 @@
1.728 ///
1.729 /// Simple constructor of the class (with lower bounds).
1.730 ///
1.731 - /// \param _graph The directed graph the algorithm runs on.
1.732 - /// \param _lower The lower bounds of the edges.
1.733 - /// \param _capacity The capacities (upper bounds) of the edges.
1.734 - /// \param _cost The cost (length) values of the edges.
1.735 - /// \param _s The source node.
1.736 - /// \param _t The target node.
1.737 - /// \param _flow_value The required amount of flow from node \c _s
1.738 - /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
1.739 - NetworkSimplex( const Graph &_graph,
1.740 - const LowerMap &_lower,
1.741 - const CapacityMap &_capacity,
1.742 - const CostMap &_cost,
1.743 - typename Graph::Node _s,
1.744 - typename Graph::Node _t,
1.745 - typename SupplyMap::Value _flow_value ) :
1.746 - graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
1.747 - supply(graph), flow(graph), flow_result(_graph), potential(graph),
1.748 - potential_result(_graph), depth(graph), parent(graph),
1.749 - pred_edge(graph), thread(graph), forward(graph), state(graph),
1.750 - node_ref(graph_ref), edge_ref(graph_ref),
1.751 - red_cost(graph, cost, potential)
1.752 + /// \param graph The directed graph the algorithm runs on.
1.753 + /// \param lower The lower bounds of the edges.
1.754 + /// \param capacity The capacities (upper bounds) of the edges.
1.755 + /// \param cost The cost (length) values of the edges.
1.756 + /// \param s The source node.
1.757 + /// \param t The target node.
1.758 + /// \param flow_value The required amount of flow from node \c s
1.759 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.760 + NetworkSimplex( const Graph &graph,
1.761 + const LowerMap &lower,
1.762 + const CapacityMap &capacity,
1.763 + const CostMap &cost,
1.764 + typename Graph::Node s,
1.765 + typename Graph::Node t,
1.766 + typename SupplyMap::Value flow_value ) :
1.767 + _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
1.768 + _cost(_graph), _supply(_graph), _flow(_graph),
1.769 + _potential(_graph), _depth(_graph), _parent(_graph),
1.770 + _pred_edge(_graph), _thread(_graph), _forward(_graph),
1.771 + _state(_graph), _red_cost(_graph, _cost, _potential),
1.772 + _flow_result(graph), _potential_result(graph),
1.773 + _node_ref(graph), _edge_ref(graph)
1.774 {
1.775 - // Copying graph_ref to graph
1.776 - copyGraph(graph, graph_ref)
1.777 - .edgeMap(cost, _cost)
1.778 - .nodeRef(node_ref)
1.779 - .edgeRef(edge_ref)
1.780 + // Copying _graph_ref to graph
1.781 + copyGraph(_graph, _graph_ref)
1.782 + .edgeMap(_cost, cost)
1.783 + .nodeRef(_node_ref)
1.784 + .edgeRef(_edge_ref)
1.785 .run();
1.786
1.787 // Removing non-zero lower bounds
1.788 - for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
1.789 - capacity[edge_ref[e]] = _capacity[e] - _lower[e];
1.790 + for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
1.791 + _capacity[_edge_ref[e]] = capacity[e] - lower[e];
1.792 }
1.793 - for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
1.794 - Supply s = 0;
1.795 - if (n == _s) s = _flow_value;
1.796 - if (n == _t) s = -_flow_value;
1.797 - for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.798 - s += _lower[e];
1.799 - for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.800 - s -= _lower[e];
1.801 - supply[node_ref[n]] = s;
1.802 + for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
1.803 + Supply sum = 0;
1.804 + if (n == s) sum = flow_value;
1.805 + if (n == t) sum = -flow_value;
1.806 + for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
1.807 + sum += lower[e];
1.808 + for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
1.809 + sum -= lower[e];
1.810 + _supply[_node_ref[n]] = sum;
1.811 }
1.812 - valid_supply = true;
1.813 + _valid_supply = true;
1.814 }
1.815
1.816 /// \brief Simple constructor of the class (without lower bounds).
1.817 ///
1.818 /// Simple constructor of the class (without lower bounds).
1.819 ///
1.820 - /// \param _graph The directed graph the algorithm runs on.
1.821 - /// \param _capacity The capacities (upper bounds) of the edges.
1.822 - /// \param _cost The cost (length) values of the edges.
1.823 - /// \param _s The source node.
1.824 - /// \param _t The target node.
1.825 - /// \param _flow_value The required amount of flow from node \c _s
1.826 - /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
1.827 - NetworkSimplex( const Graph &_graph,
1.828 - const CapacityMap &_capacity,
1.829 - const CostMap &_cost,
1.830 - typename Graph::Node _s,
1.831 - typename Graph::Node _t,
1.832 - typename SupplyMap::Value _flow_value ) :
1.833 - graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
1.834 - supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),
1.835 - potential_result(_graph), depth(graph), parent(graph),
1.836 - pred_edge(graph), thread(graph), forward(graph), state(graph),
1.837 - node_ref(graph_ref), edge_ref(graph_ref),
1.838 - red_cost(graph, cost, potential)
1.839 + /// \param graph The directed graph the algorithm runs on.
1.840 + /// \param capacity The capacities (upper bounds) of the edges.
1.841 + /// \param cost The cost (length) values of the edges.
1.842 + /// \param s The source node.
1.843 + /// \param t The target node.
1.844 + /// \param flow_value The required amount of flow from node \c s
1.845 + /// to node \c t (i.e. the supply of \c s and the demand of \c t).
1.846 + NetworkSimplex( const Graph &graph,
1.847 + const CapacityMap &capacity,
1.848 + const CostMap &cost,
1.849 + typename Graph::Node s,
1.850 + typename Graph::Node t,
1.851 + typename SupplyMap::Value flow_value ) :
1.852 + _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
1.853 + _cost(_graph), _supply(_graph, 0), _flow(_graph),
1.854 + _potential(_graph), _depth(_graph), _parent(_graph),
1.855 + _pred_edge(_graph), _thread(_graph), _forward(_graph),
1.856 + _state(_graph), _red_cost(_graph, _cost, _potential),
1.857 + _flow_result(graph), _potential_result(graph),
1.858 + _node_ref(graph), _edge_ref(graph)
1.859 {
1.860 - // Copying graph_ref to graph
1.861 - copyGraph(graph, graph_ref)
1.862 - .edgeMap(capacity, _capacity)
1.863 - .edgeMap(cost, _cost)
1.864 - .nodeRef(node_ref)
1.865 - .edgeRef(edge_ref)
1.866 + // Copying _graph_ref to graph
1.867 + copyGraph(_graph, _graph_ref)
1.868 + .edgeMap(_capacity, capacity)
1.869 + .edgeMap(_cost, cost)
1.870 + .nodeRef(_node_ref)
1.871 + .edgeRef(_edge_ref)
1.872 .run();
1.873 - supply[node_ref[_s]] = _flow_value;
1.874 - supply[node_ref[_t]] = -_flow_value;
1.875 - valid_supply = true;
1.876 + _supply[_node_ref[s]] = flow_value;
1.877 + _supply[_node_ref[t]] = -flow_value;
1.878 + _valid_supply = true;
1.879 }
1.880
1.881 /// \brief Runs the algorithm.
1.882 ///
1.883 /// Runs the algorithm.
1.884 ///
1.885 + /// \param pivot_rule The pivot rule that is used during the
1.886 + /// algorithm.
1.887 + ///
1.888 + /// The available pivot rules:
1.889 + ///
1.890 + /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
1.891 + /// a wraparound fashion in every iteration
1.892 + /// (\ref FirstEligiblePivotRule).
1.893 + ///
1.894 + /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
1.895 + /// every iteration (\ref BestEligiblePivotRule).
1.896 + ///
1.897 + /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
1.898 + /// every iteration in a wraparound fashion and the best eligible
1.899 + /// edge is selected from this block (\ref BlockSearchPivotRule).
1.900 + ///
1.901 + /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
1.902 + /// examined in every iteration in a wraparound fashion and the best
1.903 + /// one is selected from them (\ref LimitedSearchPivotRule).
1.904 + ///
1.905 + /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
1.906 + /// built from eligible edges and it is used for edge selection in
1.907 + /// the following minor iterations (\ref CandidateListPivotRule).
1.908 + ///
1.909 + /// - COMBINED_PIVOT This is a combined version of the two fastest
1.910 + /// pivot rules.
1.911 + /// For rather sparse graphs \ref LimitedSearchPivotRule
1.912 + /// "Limited Search" implementation is used, otherwise
1.913 + /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
1.914 + /// According to our benchmark tests this combined method is the
1.915 + /// most efficient.
1.916 + ///
1.917 /// \return \c true if a feasible flow can be found.
1.918 - bool run() {
1.919 - return init() && start();
1.920 + bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
1.921 + return init() && start(pivot_rule);
1.922 }
1.923
1.924 - /// \brief Returns a const reference to the flow map.
1.925 + /// \brief Returns a const reference to the edge map storing the
1.926 + /// found flow.
1.927 ///
1.928 - /// Returns a const reference to the flow map.
1.929 + /// Returns a const reference to the edge map storing the found flow.
1.930 ///
1.931 /// \pre \ref run() must be called before using this function.
1.932 const FlowMap& flowMap() const {
1.933 - return flow_result;
1.934 + return _flow_result;
1.935 }
1.936
1.937 - /// \brief Returns a const reference to the potential map (the dual
1.938 - /// solution).
1.939 + /// \brief Returns a const reference to the node map storing the
1.940 + /// found potentials (the dual solution).
1.941 ///
1.942 - /// Returns a const reference to the potential map (the dual
1.943 - /// solution).
1.944 + /// Returns a const reference to the node map storing the found
1.945 + /// potentials (the dual solution).
1.946 ///
1.947 /// \pre \ref run() must be called before using this function.
1.948 const PotentialMap& potentialMap() const {
1.949 - return potential_result;
1.950 + return _potential_result;
1.951 }
1.952
1.953 /// \brief Returns the total cost of the found flow.
1.954 @@ -445,248 +732,75 @@
1.955 /// \pre \ref run() must be called before using this function.
1.956 Cost totalCost() const {
1.957 Cost c = 0;
1.958 - for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
1.959 - c += flow_result[e] * cost[edge_ref[e]];
1.960 + for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1.961 + c += _flow_result[e] * _cost[_edge_ref[e]];
1.962 return c;
1.963 }
1.964
1.965 - protected:
1.966 + private:
1.967
1.968 /// \brief Extends the underlaying graph and initializes all the
1.969 /// node and edge maps.
1.970 bool init() {
1.971 - if (!valid_supply) return false;
1.972 + if (!_valid_supply) return false;
1.973
1.974 // Initializing state and flow maps
1.975 - for (EdgeIt e(graph); e != INVALID; ++e) {
1.976 - flow[e] = 0;
1.977 - state[e] = LOWER;
1.978 + for (EdgeIt e(_graph); e != INVALID; ++e) {
1.979 + _flow[e] = 0;
1.980 + _state[e] = STATE_LOWER;
1.981 }
1.982
1.983 // Adding an artificial root node to the graph
1.984 - root = graph.addNode();
1.985 - parent[root] = INVALID;
1.986 - pred_edge[root] = INVALID;
1.987 - depth[root] = 0;
1.988 - supply[root] = 0;
1.989 - potential[root] = 0;
1.990 + _root = _graph.addNode();
1.991 + _parent[_root] = INVALID;
1.992 + _pred_edge[_root] = INVALID;
1.993 + _depth[_root] = 0;
1.994 + _supply[_root] = 0;
1.995 + _potential[_root] = 0;
1.996
1.997 // Adding artificial edges to the graph and initializing the node
1.998 // maps of the spanning tree data structure
1.999 - Supply sum = 0;
1.1000 - Node last = root;
1.1001 + Node last = _root;
1.1002 Edge e;
1.1003 Cost max_cost = std::numeric_limits<Cost>::max() / 4;
1.1004 - for (NodeIt u(graph); u != INVALID; ++u) {
1.1005 - if (u == root) continue;
1.1006 - thread[last] = u;
1.1007 + for (NodeIt u(_graph); u != INVALID; ++u) {
1.1008 + if (u == _root) continue;
1.1009 + _thread[last] = u;
1.1010 last = u;
1.1011 - parent[u] = root;
1.1012 - depth[u] = 1;
1.1013 - sum += supply[u];
1.1014 - if (supply[u] >= 0) {
1.1015 - e = graph.addEdge(u, root);
1.1016 - flow[e] = supply[u];
1.1017 - forward[u] = true;
1.1018 - potential[u] = max_cost;
1.1019 + _parent[u] = _root;
1.1020 + _depth[u] = 1;
1.1021 + if (_supply[u] >= 0) {
1.1022 + e = _graph.addEdge(u, _root);
1.1023 + _flow[e] = _supply[u];
1.1024 + _forward[u] = true;
1.1025 + _potential[u] = -max_cost;
1.1026 } else {
1.1027 - e = graph.addEdge(root, u);
1.1028 - flow[e] = -supply[u];
1.1029 - forward[u] = false;
1.1030 - potential[u] = -max_cost;
1.1031 + e = _graph.addEdge(_root, u);
1.1032 + _flow[e] = -_supply[u];
1.1033 + _forward[u] = false;
1.1034 + _potential[u] = max_cost;
1.1035 }
1.1036 - cost[e] = max_cost;
1.1037 - capacity[e] = std::numeric_limits<Capacity>::max();
1.1038 - state[e] = TREE;
1.1039 - pred_edge[u] = e;
1.1040 + _cost[e] = max_cost;
1.1041 + _capacity[e] = std::numeric_limits<Capacity>::max();
1.1042 + _state[e] = STATE_TREE;
1.1043 + _pred_edge[u] = e;
1.1044 }
1.1045 - thread[last] = root;
1.1046 + _thread[last] = _root;
1.1047
1.1048 -#ifdef EDGE_BLOCK_PIVOT
1.1049 - // Initializing block_size for the edge block pivot rule
1.1050 - int edge_num = countEdges(graph);
1.1051 - block_size = 2 * int(sqrt(countEdges(graph)));
1.1052 - if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE;
1.1053 -#endif
1.1054 -#ifdef CANDIDATE_LIST_PIVOT
1.1055 - int edge_num = countEdges(graph);
1.1056 - minor_count = 0;
1.1057 - list_length = edge_num / LIST_LENGTH_DIV;
1.1058 - minor_limit = edge_num / MINOR_LIMIT_DIV;
1.1059 -#endif
1.1060 -#ifdef SORTED_LIST_PIVOT
1.1061 - int edge_num = countEdges(graph);
1.1062 - list_index = 0;
1.1063 - list_length = edge_num / LIST_LENGTH_DIV;
1.1064 -#endif
1.1065 -
1.1066 - return sum == 0;
1.1067 + return true;
1.1068 }
1.1069
1.1070 -#ifdef FIRST_ELIGIBLE_PIVOT
1.1071 - /// \brief Finds entering edge according to the "First Eligible"
1.1072 - /// pivot rule.
1.1073 - bool findEnteringEdge(EdgeIt &next_edge) {
1.1074 - for (EdgeIt e = next_edge; e != INVALID; ++e) {
1.1075 - if (state[e] * red_cost[e] < 0) {
1.1076 - in_edge = e;
1.1077 - next_edge = ++e;
1.1078 - return true;
1.1079 + /// Finds the join node.
1.1080 + Node findJoinNode() {
1.1081 + Node u = _graph.source(_in_edge);
1.1082 + Node v = _graph.target(_in_edge);
1.1083 + while (u != v) {
1.1084 + if (_depth[u] == _depth[v]) {
1.1085 + u = _parent[u];
1.1086 + v = _parent[v];
1.1087 }
1.1088 - }
1.1089 - for (EdgeIt e(graph); e != next_edge; ++e) {
1.1090 - if (state[e] * red_cost[e] < 0) {
1.1091 - in_edge = e;
1.1092 - next_edge = ++e;
1.1093 - return true;
1.1094 - }
1.1095 - }
1.1096 - return false;
1.1097 - }
1.1098 -#endif
1.1099 -
1.1100 -#ifdef BEST_ELIGIBLE_PIVOT
1.1101 - /// \brief Finds entering edge according to the "Best Eligible"
1.1102 - /// pivot rule.
1.1103 - bool findEnteringEdge() {
1.1104 - Cost min = 0;
1.1105 - for (EdgeIt e(graph); e != INVALID; ++e) {
1.1106 - if (state[e] * red_cost[e] < min) {
1.1107 - min = state[e] * red_cost[e];
1.1108 - in_edge = e;
1.1109 - }
1.1110 - }
1.1111 - return min < 0;
1.1112 - }
1.1113 -#endif
1.1114 -
1.1115 -#ifdef EDGE_BLOCK_PIVOT
1.1116 - /// \brief Finds entering edge according to the "Edge Block"
1.1117 - /// pivot rule.
1.1118 - bool findEnteringEdge(EdgeIt &next_edge) {
1.1119 - // Performing edge block selection
1.1120 - Cost curr, min = 0;
1.1121 - EdgeIt min_edge(graph);
1.1122 - int cnt = 0;
1.1123 - for (EdgeIt e = next_edge; e != INVALID; ++e) {
1.1124 - if ((curr = state[e] * red_cost[e]) < min) {
1.1125 - min = curr;
1.1126 - min_edge = e;
1.1127 - }
1.1128 - if (++cnt == block_size) {
1.1129 - if (min < 0) break;
1.1130 - cnt = 0;
1.1131 - }
1.1132 - }
1.1133 - if (!(min < 0)) {
1.1134 - for (EdgeIt e(graph); e != next_edge; ++e) {
1.1135 - if ((curr = state[e] * red_cost[e]) < min) {
1.1136 - min = curr;
1.1137 - min_edge = e;
1.1138 - }
1.1139 - if (++cnt == block_size) {
1.1140 - if (min < 0) break;
1.1141 - cnt = 0;
1.1142 - }
1.1143 - }
1.1144 - }
1.1145 - in_edge = min_edge;
1.1146 - if ((next_edge = ++min_edge) == INVALID)
1.1147 - next_edge = EdgeIt(graph);
1.1148 - return min < 0;
1.1149 - }
1.1150 -#endif
1.1151 -
1.1152 -#ifdef CANDIDATE_LIST_PIVOT
1.1153 - /// \brief Finds entering edge according to the "Candidate List"
1.1154 - /// pivot rule.
1.1155 - bool findEnteringEdge() {
1.1156 - typedef typename std::vector<Edge>::iterator ListIt;
1.1157 -
1.1158 - if (minor_count >= minor_limit || candidates.size() == 0) {
1.1159 - // Major iteration
1.1160 - candidates.clear();
1.1161 - for (EdgeIt e(graph); e != INVALID; ++e) {
1.1162 - if (state[e] * red_cost[e] < 0) {
1.1163 - candidates.push_back(e);
1.1164 - if (candidates.size() == list_length) break;
1.1165 - }
1.1166 - }
1.1167 - if (candidates.size() == 0) return false;
1.1168 - }
1.1169 -
1.1170 - // Minor iteration
1.1171 - ++minor_count;
1.1172 - Cost min = 0;
1.1173 - Edge e;
1.1174 - for (int i = 0; i < candidates.size(); ++i) {
1.1175 - e = candidates[i];
1.1176 - if (state[e] * red_cost[e] < min) {
1.1177 - min = state[e] * red_cost[e];
1.1178 - in_edge = e;
1.1179 - }
1.1180 - }
1.1181 - return true;
1.1182 - }
1.1183 -#endif
1.1184 -
1.1185 -#ifdef SORTED_LIST_PIVOT
1.1186 - /// \brief Functor class to compare edges during sort of the
1.1187 - /// candidate list.
1.1188 - class SortFunc
1.1189 - {
1.1190 - private:
1.1191 - const IntEdgeMap &st;
1.1192 - const ReducedCostMap &rc;
1.1193 - public:
1.1194 - SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
1.1195 - st(_st), rc(_rc) {}
1.1196 - bool operator()(const Edge &e1, const Edge &e2) {
1.1197 - return st[e1] * rc[e1] < st[e2] * rc[e2];
1.1198 - }
1.1199 - };
1.1200 -
1.1201 - /// \brief Finds entering edge according to the "Sorted List"
1.1202 - /// pivot rule.
1.1203 - bool findEnteringEdge() {
1.1204 - static SortFunc sort_func(state, red_cost);
1.1205 -
1.1206 - // Minor iteration
1.1207 - while (list_index < candidates.size()) {
1.1208 - in_edge = candidates[list_index++];
1.1209 - if (state[in_edge] * red_cost[in_edge] < 0) return true;
1.1210 - }
1.1211 -
1.1212 - // Major iteration
1.1213 - candidates.clear();
1.1214 - Cost curr, min = 0;
1.1215 - for (EdgeIt e(graph); e != INVALID; ++e) {
1.1216 - if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
1.1217 - candidates.push_back(e);
1.1218 - if (curr < min) min = curr;
1.1219 - if (candidates.size() == list_length) break;
1.1220 - }
1.1221 - }
1.1222 - if (candidates.size() == 0) return false;
1.1223 - sort(candidates.begin(), candidates.end(), sort_func);
1.1224 - in_edge = candidates[0];
1.1225 - list_index = 1;
1.1226 - return true;
1.1227 - }
1.1228 -#endif
1.1229 -
1.1230 - /// \brief Finds the join node.
1.1231 - Node findJoinNode() {
1.1232 - // Finding the join node
1.1233 - Node u = graph.source(in_edge);
1.1234 - Node v = graph.target(in_edge);
1.1235 - while (u != v) {
1.1236 - if (depth[u] == depth[v]) {
1.1237 - u = parent[u];
1.1238 - v = parent[v];
1.1239 - }
1.1240 - else if (depth[u] > depth[v]) u = parent[u];
1.1241 - else v = parent[v];
1.1242 + else if (_depth[u] > _depth[v]) u = _parent[u];
1.1243 + else v = _parent[v];
1.1244 }
1.1245 return u;
1.1246 }
1.1247 @@ -696,23 +810,23 @@
1.1248 bool findLeavingEdge() {
1.1249 // Initializing first and second nodes according to the direction
1.1250 // of the cycle
1.1251 - if (state[in_edge] == LOWER) {
1.1252 - first = graph.source(in_edge);
1.1253 - second = graph.target(in_edge);
1.1254 + if (_state[_in_edge] == STATE_LOWER) {
1.1255 + first = _graph.source(_in_edge);
1.1256 + second = _graph.target(_in_edge);
1.1257 } else {
1.1258 - first = graph.target(in_edge);
1.1259 - second = graph.source(in_edge);
1.1260 + first = _graph.target(_in_edge);
1.1261 + second = _graph.source(_in_edge);
1.1262 }
1.1263 - delta = capacity[in_edge];
1.1264 + delta = _capacity[_in_edge];
1.1265 bool result = false;
1.1266 Capacity d;
1.1267 Edge e;
1.1268
1.1269 // Searching the cycle along the path form the first node to the
1.1270 // root node
1.1271 - for (Node u = first; u != join; u = parent[u]) {
1.1272 - e = pred_edge[u];
1.1273 - d = forward[u] ? flow[e] : capacity[e] - flow[e];
1.1274 + for (Node u = first; u != join; u = _parent[u]) {
1.1275 + e = _pred_edge[u];
1.1276 + d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
1.1277 if (d < delta) {
1.1278 delta = d;
1.1279 u_out = u;
1.1280 @@ -723,9 +837,9 @@
1.1281 }
1.1282 // Searching the cycle along the path form the second node to the
1.1283 // root node
1.1284 - for (Node u = second; u != join; u = parent[u]) {
1.1285 - e = pred_edge[u];
1.1286 - d = forward[u] ? capacity[e] - flow[e] : flow[e];
1.1287 + for (Node u = second; u != join; u = _parent[u]) {
1.1288 + e = _pred_edge[u];
1.1289 + d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
1.1290 if (d <= delta) {
1.1291 delta = d;
1.1292 u_out = u;
1.1293 @@ -737,134 +851,150 @@
1.1294 return result;
1.1295 }
1.1296
1.1297 - /// \brief Changes \c flow and \c state edge maps.
1.1298 + /// Changes \c flow and \c state edge maps.
1.1299 void changeFlows(bool change) {
1.1300 // Augmenting along the cycle
1.1301 if (delta > 0) {
1.1302 - Capacity val = state[in_edge] * delta;
1.1303 - flow[in_edge] += val;
1.1304 - for (Node u = graph.source(in_edge); u != join; u = parent[u]) {
1.1305 - flow[pred_edge[u]] += forward[u] ? -val : val;
1.1306 + Capacity val = _state[_in_edge] * delta;
1.1307 + _flow[_in_edge] += val;
1.1308 + for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
1.1309 + _flow[_pred_edge[u]] += _forward[u] ? -val : val;
1.1310 }
1.1311 - for (Node u = graph.target(in_edge); u != join; u = parent[u]) {
1.1312 - flow[pred_edge[u]] += forward[u] ? val : -val;
1.1313 + for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
1.1314 + _flow[_pred_edge[u]] += _forward[u] ? val : -val;
1.1315 }
1.1316 }
1.1317 // Updating the state of the entering and leaving edges
1.1318 if (change) {
1.1319 - state[in_edge] = TREE;
1.1320 - state[pred_edge[u_out]] =
1.1321 - (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
1.1322 + _state[_in_edge] = STATE_TREE;
1.1323 + _state[_pred_edge[u_out]] =
1.1324 + (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1.1325 } else {
1.1326 - state[in_edge] = -state[in_edge];
1.1327 + _state[_in_edge] = -_state[_in_edge];
1.1328 }
1.1329 }
1.1330
1.1331 - /// \brief Updates \c thread and \c parent node maps.
1.1332 + /// Updates \c thread and \c parent node maps.
1.1333 void updateThreadParent() {
1.1334 Node u;
1.1335 - v_out = parent[u_out];
1.1336 + v_out = _parent[u_out];
1.1337
1.1338 // Handling the case when join and v_out coincide
1.1339 bool par_first = false;
1.1340 if (join == v_out) {
1.1341 - for (u = join; u != u_in && u != v_in; u = thread[u]) ;
1.1342 + for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
1.1343 if (u == v_in) {
1.1344 par_first = true;
1.1345 - while (thread[u] != u_out) u = thread[u];
1.1346 + while (_thread[u] != u_out) u = _thread[u];
1.1347 first = u;
1.1348 }
1.1349 }
1.1350
1.1351 // Finding the last successor of u_in (u) and the node after it
1.1352 // (right) according to the thread index
1.1353 - for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;
1.1354 - right = thread[u];
1.1355 - if (thread[v_in] == u_out) {
1.1356 - for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
1.1357 - if (last == u_out) last = thread[last];
1.1358 + for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
1.1359 + right = _thread[u];
1.1360 + if (_thread[v_in] == u_out) {
1.1361 + for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
1.1362 + if (last == u_out) last = _thread[last];
1.1363 }
1.1364 - else last = thread[v_in];
1.1365 + else last = _thread[v_in];
1.1366
1.1367 // Updating stem nodes
1.1368 - thread[v_in] = stem = u_in;
1.1369 + _thread[v_in] = stem = u_in;
1.1370 par_stem = v_in;
1.1371 while (stem != u_out) {
1.1372 - thread[u] = new_stem = parent[stem];
1.1373 + _thread[u] = new_stem = _parent[stem];
1.1374
1.1375 // Finding the node just before the stem node (u) according to
1.1376 // the original thread index
1.1377 - for (u = new_stem; thread[u] != stem; u = thread[u]) ;
1.1378 - thread[u] = right;
1.1379 + for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
1.1380 + _thread[u] = right;
1.1381
1.1382 // Changing the parent node of stem and shifting stem and
1.1383 // par_stem nodes
1.1384 - parent[stem] = par_stem;
1.1385 + _parent[stem] = par_stem;
1.1386 par_stem = stem;
1.1387 stem = new_stem;
1.1388
1.1389 // Finding the last successor of stem (u) and the node after it
1.1390 // (right) according to the thread index
1.1391 - for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
1.1392 - right = thread[u];
1.1393 + for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
1.1394 + right = _thread[u];
1.1395 }
1.1396 - parent[u_out] = par_stem;
1.1397 - thread[u] = last;
1.1398 + _parent[u_out] = par_stem;
1.1399 + _thread[u] = last;
1.1400
1.1401 if (join == v_out && par_first) {
1.1402 - if (first != v_in) thread[first] = right;
1.1403 + if (first != v_in) _thread[first] = right;
1.1404 } else {
1.1405 - for (u = v_out; thread[u] != u_out; u = thread[u]) ;
1.1406 - thread[u] = right;
1.1407 + for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
1.1408 + _thread[u] = right;
1.1409 }
1.1410 }
1.1411
1.1412 - /// \brief Updates \c pred_edge and \c forward node maps.
1.1413 + /// Updates \c pred_edge and \c forward node maps.
1.1414 void updatePredEdge() {
1.1415 Node u = u_out, v;
1.1416 while (u != u_in) {
1.1417 - v = parent[u];
1.1418 - pred_edge[u] = pred_edge[v];
1.1419 - forward[u] = !forward[v];
1.1420 + v = _parent[u];
1.1421 + _pred_edge[u] = _pred_edge[v];
1.1422 + _forward[u] = !_forward[v];
1.1423 u = v;
1.1424 }
1.1425 - pred_edge[u_in] = in_edge;
1.1426 - forward[u_in] = (u_in == graph.source(in_edge));
1.1427 + _pred_edge[u_in] = _in_edge;
1.1428 + _forward[u_in] = (u_in == _graph.source(_in_edge));
1.1429 }
1.1430
1.1431 - /// \brief Updates \c depth and \c potential node maps.
1.1432 + /// Updates \c depth and \c potential node maps.
1.1433 void updateDepthPotential() {
1.1434 - depth[u_in] = depth[v_in] + 1;
1.1435 - potential[u_in] = forward[u_in] ?
1.1436 - potential[v_in] + cost[pred_edge[u_in]] :
1.1437 - potential[v_in] - cost[pred_edge[u_in]];
1.1438 + _depth[u_in] = _depth[v_in] + 1;
1.1439 + _potential[u_in] = _forward[u_in] ?
1.1440 + _potential[v_in] - _cost[_pred_edge[u_in]] :
1.1441 + _potential[v_in] + _cost[_pred_edge[u_in]];
1.1442
1.1443 - Node u = thread[u_in], v;
1.1444 + Node u = _thread[u_in], v;
1.1445 while (true) {
1.1446 - v = parent[u];
1.1447 + v = _parent[u];
1.1448 if (v == INVALID) break;
1.1449 - depth[u] = depth[v] + 1;
1.1450 - potential[u] = forward[u] ?
1.1451 - potential[v] + cost[pred_edge[u]] :
1.1452 - potential[v] - cost[pred_edge[u]];
1.1453 - if (depth[u] <= depth[v_in]) break;
1.1454 - u = thread[u];
1.1455 + _depth[u] = _depth[v] + 1;
1.1456 + _potential[u] = _forward[u] ?
1.1457 + _potential[v] - _cost[_pred_edge[u]] :
1.1458 + _potential[v] + _cost[_pred_edge[u]];
1.1459 + if (_depth[u] <= _depth[v_in]) break;
1.1460 + u = _thread[u];
1.1461 }
1.1462 }
1.1463
1.1464 - /// \brief Executes the algorithm.
1.1465 + /// Executes the algorithm.
1.1466 + bool start(PivotRuleEnum pivot_rule) {
1.1467 + switch (pivot_rule) {
1.1468 + case FIRST_ELIGIBLE_PIVOT:
1.1469 + return start<FirstEligiblePivotRule>();
1.1470 + case BEST_ELIGIBLE_PIVOT:
1.1471 + return start<BestEligiblePivotRule>();
1.1472 + case BLOCK_SEARCH_PIVOT:
1.1473 + return start<BlockSearchPivotRule>();
1.1474 + case LIMITED_SEARCH_PIVOT:
1.1475 + return start<LimitedSearchPivotRule>();
1.1476 + case CANDIDATE_LIST_PIVOT:
1.1477 + return start<CandidateListPivotRule>();
1.1478 + case COMBINED_PIVOT:
1.1479 + if ( countEdges(_graph) / countNodes(_graph) <=
1.1480 + COMBINED_PIVOT_MAX_DEG )
1.1481 + return start<LimitedSearchPivotRule>();
1.1482 + else
1.1483 + return start<BlockSearchPivotRule>();
1.1484 + }
1.1485 + return false;
1.1486 + }
1.1487 +
1.1488 + template<class PivotRuleImplementation>
1.1489 bool start() {
1.1490 - // Processing pivots
1.1491 -#ifdef _DEBUG_ITER_
1.1492 - int iter_num = 0;
1.1493 -#endif
1.1494 -#if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT)
1.1495 - EdgeIt next_edge(graph);
1.1496 - while (findEnteringEdge(next_edge))
1.1497 -#else
1.1498 - while (findEnteringEdge())
1.1499 -#endif
1.1500 - {
1.1501 + PivotRuleImplementation pivot(*this);
1.1502 +
1.1503 + // Executing the network simplex algorithm
1.1504 + while (pivot.findEnteringEdge()) {
1.1505 join = findJoinNode();
1.1506 bool change = findLeavingEdge();
1.1507 changeFlows(change);
1.1508 @@ -873,34 +1003,26 @@
1.1509 updatePredEdge();
1.1510 updateDepthPotential();
1.1511 }
1.1512 -#ifdef _DEBUG_ITER_
1.1513 - ++iter_num;
1.1514 -#endif
1.1515 }
1.1516
1.1517 -#ifdef _DEBUG_ITER_
1.1518 - std::cout << "Network Simplex algorithm finished. " << iter_num
1.1519 - << " pivot iterations performed." << std::endl;
1.1520 -#endif
1.1521 + // Checking if the flow amount equals zero on all the artificial
1.1522 + // edges
1.1523 + for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
1.1524 + if (_flow[e] > 0) return false;
1.1525 + for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
1.1526 + if (_flow[e] > 0) return false;
1.1527
1.1528 - // Checking if the flow amount equals zero on all the
1.1529 - // artificial edges
1.1530 - for (InEdgeIt e(graph, root); e != INVALID; ++e)
1.1531 - if (flow[e] > 0) return false;
1.1532 - for (OutEdgeIt e(graph, root); e != INVALID; ++e)
1.1533 - if (flow[e] > 0) return false;
1.1534 -
1.1535 - // Copying flow values to flow_result
1.1536 - if (lower) {
1.1537 - for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
1.1538 - flow_result[e] = (*lower)[e] + flow[edge_ref[e]];
1.1539 + // Copying flow values to _flow_result
1.1540 + if (_lower) {
1.1541 + for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1.1542 + _flow_result[e] = (*_lower)[e] + _flow[_edge_ref[e]];
1.1543 } else {
1.1544 - for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
1.1545 - flow_result[e] = flow[edge_ref[e]];
1.1546 + for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1.1547 + _flow_result[e] = _flow[_edge_ref[e]];
1.1548 }
1.1549 - // Copying potential values to potential_result
1.1550 - for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
1.1551 - potential_result[n] = potential[node_ref[n]];
1.1552 + // Copying potential values to _potential_result
1.1553 + for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1.1554 + _potential_result[n] = _potential[_node_ref[n]];
1.1555
1.1556 return true;
1.1557 }