1.1 --- a/lemon/network_simplex.h Wed Feb 04 14:42:31 2009 +0000
1.2 +++ b/lemon/network_simplex.h Fri Feb 06 21:52:34 2009 +0000
1.3 @@ -28,9 +28,7 @@
1.4 #include <limits>
1.5 #include <algorithm>
1.6
1.7 -#include <lemon/graph_adaptor.h>
1.8 #include <lemon/graph_utils.h>
1.9 -#include <lemon/smart_graph.h>
1.10 #include <lemon/math.h>
1.11
1.12 namespace lemon {
1.13 @@ -72,29 +70,21 @@
1.14 typename SupplyMap = typename Graph::template NodeMap<int> >
1.15 class NetworkSimplex
1.16 {
1.17 + GRAPH_TYPEDEFS(typename Graph);
1.18 +
1.19 typedef typename CapacityMap::Value Capacity;
1.20 typedef typename CostMap::Value Cost;
1.21 typedef typename SupplyMap::Value Supply;
1.22
1.23 - typedef SmartGraph SGraph;
1.24 - GRAPH_TYPEDEFS(typename SGraph);
1.25 + typedef std::vector<Edge> EdgeVector;
1.26 + typedef std::vector<Node> NodeVector;
1.27 + typedef std::vector<int> IntVector;
1.28 + typedef std::vector<bool> BoolVector;
1.29 + typedef std::vector<Capacity> CapacityVector;
1.30 + typedef std::vector<Cost> CostVector;
1.31 + typedef std::vector<Supply> SupplyVector;
1.32
1.33 - typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
1.34 - typedef typename SGraph::template EdgeMap<Cost> SCostMap;
1.35 - typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
1.36 - typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
1.37 -
1.38 - typedef typename SGraph::template NodeMap<int> IntNodeMap;
1.39 - typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
1.40 - typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
1.41 - typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
1.42 - typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
1.43 - typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;
1.44 -
1.45 - typedef typename Graph::template NodeMap<Node> NodeRefMap;
1.46 - typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
1.47 -
1.48 - typedef std::vector<Edge> EdgeVector;
1.49 + typedef typename Graph::template NodeMap<int> IntNodeMap;
1.50
1.51 public:
1.52
1.53 @@ -116,35 +106,6 @@
1.54
1.55 private:
1.56
1.57 - /// \brief Map adaptor class for handling reduced edge costs.
1.58 - ///
1.59 - /// Map adaptor class for handling reduced edge costs.
1.60 - class ReducedCostMap : public MapBase<Edge, Cost>
1.61 - {
1.62 - private:
1.63 -
1.64 - const SGraph &_gr;
1.65 - const SCostMap &_cost_map;
1.66 - const SPotentialMap &_pot_map;
1.67 -
1.68 - public:
1.69 -
1.70 - ///\e
1.71 - ReducedCostMap( const SGraph &gr,
1.72 - const SCostMap &cost_map,
1.73 - const SPotentialMap &pot_map ) :
1.74 - _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
1.75 -
1.76 - ///\e
1.77 - Cost operator[](const Edge &e) const {
1.78 - return _cost_map[e] + _pot_map[_gr.source(e)]
1.79 - - _pot_map[_gr.target(e)];
1.80 - }
1.81 -
1.82 - }; //class ReducedCostMap
1.83 -
1.84 - private:
1.85 -
1.86 /// \brief Implementation of the "First Eligible" pivot rule for the
1.87 /// \ref NetworkSimplex "network simplex" algorithm.
1.88 ///
1.89 @@ -157,40 +118,52 @@
1.90 private:
1.91
1.92 // References to the NetworkSimplex class
1.93 - NetworkSimplex &_ns;
1.94 - EdgeVector &_edges;
1.95 + const EdgeVector &_edge;
1.96 + const IntVector &_source;
1.97 + const IntVector &_target;
1.98 + const CostVector &_cost;
1.99 + const IntVector &_state;
1.100 + const CostVector &_pi;
1.101 + int &_in_edge;
1.102 + int _edge_num;
1.103
1.104 + // Pivot rule data
1.105 int _next_edge;
1.106
1.107 public:
1.108
1.109 /// Constructor
1.110 - FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.111 - _ns(ns), _edges(edges), _next_edge(0) {}
1.112 + FirstEligiblePivotRule(NetworkSimplex &ns) :
1.113 + _edge(ns._edge), _source(ns._source), _target(ns._target),
1.114 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.115 + _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
1.116 + {}
1.117
1.118 /// Find next entering edge
1.119 bool findEnteringEdge() {
1.120 - Edge e;
1.121 - for (int i = _next_edge; i < int(_edges.size()); ++i) {
1.122 - e = _edges[i];
1.123 - if (_ns._state[e] * _ns._red_cost[e] < 0) {
1.124 - _ns._in_edge = e;
1.125 - _next_edge = i + 1;
1.126 + Cost c;
1.127 + for (int e = _next_edge; e < _edge_num; ++e) {
1.128 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.129 + if (c < 0) {
1.130 + _in_edge = e;
1.131 + _next_edge = e + 1;
1.132 return true;
1.133 }
1.134 }
1.135 - for (int i = 0; i < _next_edge; ++i) {
1.136 - e = _edges[i];
1.137 - if (_ns._state[e] * _ns._red_cost[e] < 0) {
1.138 - _ns._in_edge = e;
1.139 - _next_edge = i + 1;
1.140 + for (int e = 0; e < _next_edge; ++e) {
1.141 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.142 + if (c < 0) {
1.143 + _in_edge = e;
1.144 + _next_edge = e + 1;
1.145 return true;
1.146 }
1.147 }
1.148 return false;
1.149 }
1.150 +
1.151 }; //class FirstEligiblePivotRule
1.152
1.153 +
1.154 /// \brief Implementation of the "Best Eligible" pivot rule for the
1.155 /// \ref NetworkSimplex "network simplex" algorithm.
1.156 ///
1.157 @@ -203,30 +176,40 @@
1.158 private:
1.159
1.160 // References to the NetworkSimplex class
1.161 - NetworkSimplex &_ns;
1.162 - EdgeVector &_edges;
1.163 + const EdgeVector &_edge;
1.164 + const IntVector &_source;
1.165 + const IntVector &_target;
1.166 + const CostVector &_cost;
1.167 + const IntVector &_state;
1.168 + const CostVector &_pi;
1.169 + int &_in_edge;
1.170 + int _edge_num;
1.171
1.172 public:
1.173
1.174 /// Constructor
1.175 - BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.176 - _ns(ns), _edges(edges) {}
1.177 + BestEligiblePivotRule(NetworkSimplex &ns) :
1.178 + _edge(ns._edge), _source(ns._source), _target(ns._target),
1.179 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.180 + _in_edge(ns._in_edge), _edge_num(ns._edge_num)
1.181 + {}
1.182
1.183 /// Find next entering edge
1.184 bool findEnteringEdge() {
1.185 - Cost min = 0;
1.186 - Edge e;
1.187 - for (int i = 0; i < int(_edges.size()); ++i) {
1.188 - e = _edges[i];
1.189 - if (_ns._state[e] * _ns._red_cost[e] < min) {
1.190 - min = _ns._state[e] * _ns._red_cost[e];
1.191 - _ns._in_edge = e;
1.192 + Cost c, min = 0;
1.193 + for (int e = 0; e < _edge_num; ++e) {
1.194 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.195 + if (c < min) {
1.196 + min = c;
1.197 + _in_edge = e;
1.198 }
1.199 }
1.200 return min < 0;
1.201 }
1.202 +
1.203 }; //class BestEligiblePivotRule
1.204
1.205 +
1.206 /// \brief Implementation of the "Block Search" pivot rule for the
1.207 /// \ref NetworkSimplex "network simplex" algorithm.
1.208 ///
1.209 @@ -239,37 +222,45 @@
1.210 private:
1.211
1.212 // References to the NetworkSimplex class
1.213 - NetworkSimplex &_ns;
1.214 - EdgeVector &_edges;
1.215 + const EdgeVector &_edge;
1.216 + const IntVector &_source;
1.217 + const IntVector &_target;
1.218 + const CostVector &_cost;
1.219 + const IntVector &_state;
1.220 + const CostVector &_pi;
1.221 + int &_in_edge;
1.222 + int _edge_num;
1.223
1.224 + // Pivot rule data
1.225 int _block_size;
1.226 - int _next_edge, _min_edge;
1.227 + int _next_edge;
1.228
1.229 public:
1.230
1.231 /// Constructor
1.232 - BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.233 - _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
1.234 + BlockSearchPivotRule(NetworkSimplex &ns) :
1.235 + _edge(ns._edge), _source(ns._source), _target(ns._target),
1.236 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.237 + _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
1.238 {
1.239 // The main parameters of the pivot rule
1.240 const double BLOCK_SIZE_FACTOR = 2.0;
1.241 const int MIN_BLOCK_SIZE = 10;
1.242
1.243 - _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
1.244 + _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
1.245 MIN_BLOCK_SIZE );
1.246 }
1.247
1.248 /// Find next entering edge
1.249 bool findEnteringEdge() {
1.250 - Cost curr, min = 0;
1.251 - Edge e;
1.252 + Cost c, min = 0;
1.253 int cnt = _block_size;
1.254 - int i;
1.255 - for (i = _next_edge; i < int(_edges.size()); ++i) {
1.256 - e = _edges[i];
1.257 - if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.258 - min = curr;
1.259 - _min_edge = i;
1.260 + int e, min_edge = _next_edge;
1.261 + for (e = _next_edge; e < _edge_num; ++e) {
1.262 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.263 + if (c < min) {
1.264 + min = c;
1.265 + min_edge = e;
1.266 }
1.267 if (--cnt == 0) {
1.268 if (min < 0) break;
1.269 @@ -277,11 +268,11 @@
1.270 }
1.271 }
1.272 if (min == 0 || cnt > 0) {
1.273 - for (i = 0; i < _next_edge; ++i) {
1.274 - e = _edges[i];
1.275 - if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
1.276 - min = curr;
1.277 - _min_edge = i;
1.278 + for (e = 0; e < _next_edge; ++e) {
1.279 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.280 + if (c < min) {
1.281 + min = c;
1.282 + min_edge = e;
1.283 }
1.284 if (--cnt == 0) {
1.285 if (min < 0) break;
1.286 @@ -290,12 +281,14 @@
1.287 }
1.288 }
1.289 if (min >= 0) return false;
1.290 - _ns._in_edge = _edges[_min_edge];
1.291 - _next_edge = i;
1.292 + _in_edge = min_edge;
1.293 + _next_edge = e;
1.294 return true;
1.295 }
1.296 +
1.297 }; //class BlockSearchPivotRule
1.298
1.299 +
1.300 /// \brief Implementation of the "Candidate List" pivot rule for the
1.301 /// \ref NetworkSimplex "network simplex" algorithm.
1.302 ///
1.303 @@ -308,19 +301,28 @@
1.304 private:
1.305
1.306 // References to the NetworkSimplex class
1.307 - NetworkSimplex &_ns;
1.308 - EdgeVector &_edges;
1.309 + const EdgeVector &_edge;
1.310 + const IntVector &_source;
1.311 + const IntVector &_target;
1.312 + const CostVector &_cost;
1.313 + const IntVector &_state;
1.314 + const CostVector &_pi;
1.315 + int &_in_edge;
1.316 + int _edge_num;
1.317
1.318 - EdgeVector _candidates;
1.319 + // Pivot rule data
1.320 + IntVector _candidates;
1.321 int _list_length, _minor_limit;
1.322 int _curr_length, _minor_count;
1.323 - int _next_edge, _min_edge;
1.324 + int _next_edge;
1.325
1.326 public:
1.327
1.328 /// Constructor
1.329 - CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.330 - _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
1.331 + CandidateListPivotRule(NetworkSimplex &ns) :
1.332 + _edge(ns._edge), _source(ns._source), _target(ns._target),
1.333 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.334 + _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
1.335 {
1.336 // The main parameters of the pivot rule
1.337 const double LIST_LENGTH_FACTOR = 1.0;
1.338 @@ -328,7 +330,7 @@
1.339 const double MINOR_LIMIT_FACTOR = 0.1;
1.340 const int MIN_MINOR_LIMIT = 3;
1.341
1.342 - _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())),
1.343 + _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge_num)),
1.344 MIN_LIST_LENGTH );
1.345 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
1.346 MIN_MINOR_LIMIT );
1.347 @@ -338,51 +340,52 @@
1.348
1.349 /// Find next entering edge
1.350 bool findEnteringEdge() {
1.351 - Cost min, curr;
1.352 + Cost min, c;
1.353 + int e, min_edge = _next_edge;
1.354 if (_curr_length > 0 && _minor_count < _minor_limit) {
1.355 // Minor iteration: select the best eligible edge from the
1.356 // current candidate list
1.357 ++_minor_count;
1.358 - Edge e;
1.359 min = 0;
1.360 for (int i = 0; i < _curr_length; ++i) {
1.361 e = _candidates[i];
1.362 - curr = _ns._state[e] * _ns._red_cost[e];
1.363 - if (curr < min) {
1.364 - min = curr;
1.365 - _ns._in_edge = e;
1.366 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.367 + if (c < min) {
1.368 + min = c;
1.369 + min_edge = e;
1.370 }
1.371 - if (curr >= 0) {
1.372 + if (c >= 0) {
1.373 _candidates[i--] = _candidates[--_curr_length];
1.374 }
1.375 }
1.376 - if (min < 0) return true;
1.377 + if (min < 0) {
1.378 + _in_edge = min_edge;
1.379 + return true;
1.380 + }
1.381 }
1.382
1.383 // Major iteration: build a new candidate list
1.384 - Edge e;
1.385 min = 0;
1.386 _curr_length = 0;
1.387 - int i;
1.388 - for (i = _next_edge; i < int(_edges.size()); ++i) {
1.389 - e = _edges[i];
1.390 - if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.391 + for (e = _next_edge; e < _edge_num; ++e) {
1.392 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.393 + if (c < 0) {
1.394 _candidates[_curr_length++] = e;
1.395 - if (curr < min) {
1.396 - min = curr;
1.397 - _min_edge = i;
1.398 + if (c < min) {
1.399 + min = c;
1.400 + min_edge = e;
1.401 }
1.402 if (_curr_length == _list_length) break;
1.403 }
1.404 }
1.405 if (_curr_length < _list_length) {
1.406 - for (i = 0; i < _next_edge; ++i) {
1.407 - e = _edges[i];
1.408 - if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.409 + for (e = 0; e < _next_edge; ++e) {
1.410 + c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.411 + if (c < 0) {
1.412 _candidates[_curr_length++] = e;
1.413 - if (curr < min) {
1.414 - min = curr;
1.415 - _min_edge = i;
1.416 + if (c < min) {
1.417 + min = c;
1.418 + min_edge = e;
1.419 }
1.420 if (_curr_length == _list_length) break;
1.421 }
1.422 @@ -390,12 +393,14 @@
1.423 }
1.424 if (_curr_length == 0) return false;
1.425 _minor_count = 1;
1.426 - _ns._in_edge = _edges[_min_edge];
1.427 - _next_edge = i;
1.428 + _in_edge = min_edge;
1.429 + _next_edge = e;
1.430 return true;
1.431 }
1.432 +
1.433 }; //class CandidateListPivotRule
1.434
1.435 +
1.436 /// \brief Implementation of the "Altering Candidate List" pivot rule
1.437 /// for the \ref NetworkSimplex "network simplex" algorithm.
1.438 ///
1.439 @@ -408,23 +413,29 @@
1.440 private:
1.441
1.442 // References to the NetworkSimplex class
1.443 - NetworkSimplex &_ns;
1.444 - EdgeVector &_edges;
1.445 + const EdgeVector &_edge;
1.446 + const IntVector &_source;
1.447 + const IntVector &_target;
1.448 + const CostVector &_cost;
1.449 + const IntVector &_state;
1.450 + const CostVector &_pi;
1.451 + int &_in_edge;
1.452 + int _edge_num;
1.453
1.454 - EdgeVector _candidates;
1.455 - SCostMap _cand_cost;
1.456 int _block_size, _head_length, _curr_length;
1.457 int _next_edge;
1.458 + IntVector _candidates;
1.459 + CostVector _cand_cost;
1.460
1.461 // Functor class to compare edges during sort of the candidate list
1.462 class SortFunc
1.463 {
1.464 private:
1.465 - const SCostMap &_map;
1.466 + const CostVector &_map;
1.467 public:
1.468 - SortFunc(const SCostMap &map) : _map(map) {}
1.469 - bool operator()(const Edge &e1, const Edge &e2) {
1.470 - return _map[e1] > _map[e2];
1.471 + SortFunc(const CostVector &map) : _map(map) {}
1.472 + bool operator()(int left, int right) {
1.473 + return _map[left] > _map[right];
1.474 }
1.475 };
1.476
1.477 @@ -433,9 +444,11 @@
1.478 public:
1.479
1.480 /// Constructor
1.481 - AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
1.482 - _ns(ns), _edges(edges), _cand_cost(_ns._graph),
1.483 - _next_edge(0), _sort_func(_cand_cost)
1.484 + AlteringListPivotRule(NetworkSimplex &ns) :
1.485 + _edge(ns._edge), _source(ns._source), _target(ns._target),
1.486 + _cost(ns._cost), _state(ns._state), _pi(ns._pi),
1.487 + _in_edge(ns._in_edge), _edge_num(ns._edge_num),
1.488 + _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost)
1.489 {
1.490 // The main parameters of the pivot rule
1.491 const double BLOCK_SIZE_FACTOR = 1.5;
1.492 @@ -443,7 +456,7 @@
1.493 const double HEAD_LENGTH_FACTOR = 0.1;
1.494 const int MIN_HEAD_LENGTH = 3;
1.495
1.496 - _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
1.497 + _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
1.498 MIN_BLOCK_SIZE );
1.499 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
1.500 MIN_HEAD_LENGTH );
1.501 @@ -454,11 +467,13 @@
1.502 /// Find next entering edge
1.503 bool findEnteringEdge() {
1.504 // Check the current candidate list
1.505 - Edge e;
1.506 - for (int ix = 0; ix < _curr_length; ++ix) {
1.507 - e = _candidates[ix];
1.508 - if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
1.509 - _candidates[ix--] = _candidates[--_curr_length];
1.510 + int e;
1.511 + for (int i = 0; i < _curr_length; ++i) {
1.512 + e = _candidates[i];
1.513 + _cand_cost[e] = _state[e] *
1.514 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.515 + if (_cand_cost[e] >= 0) {
1.516 + _candidates[i--] = _candidates[--_curr_length];
1.517 }
1.518 }
1.519
1.520 @@ -466,11 +481,13 @@
1.521 int cnt = _block_size;
1.522 int last_edge = 0;
1.523 int limit = _head_length;
1.524 - for (int i = _next_edge; i < int(_edges.size()); ++i) {
1.525 - e = _edges[i];
1.526 - if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.527 +
1.528 + for (int e = _next_edge; e < _edge_num; ++e) {
1.529 + _cand_cost[e] = _state[e] *
1.530 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.531 + if (_cand_cost[e] < 0) {
1.532 _candidates[_curr_length++] = e;
1.533 - last_edge = i;
1.534 + last_edge = e;
1.535 }
1.536 if (--cnt == 0) {
1.537 if (_curr_length > limit) break;
1.538 @@ -479,11 +496,12 @@
1.539 }
1.540 }
1.541 if (_curr_length <= limit) {
1.542 - for (int i = 0; i < _next_edge; ++i) {
1.543 - e = _edges[i];
1.544 - if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
1.545 + for (int e = 0; e < _next_edge; ++e) {
1.546 + _cand_cost[e] = _state[e] *
1.547 + (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
1.548 + if (_cand_cost[e] < 0) {
1.549 _candidates[_curr_length++] = e;
1.550 - last_edge = i;
1.551 + last_edge = e;
1.552 }
1.553 if (--cnt == 0) {
1.554 if (_curr_length > limit) break;
1.555 @@ -500,12 +518,13 @@
1.556 _sort_func );
1.557
1.558 // Pop the first element of the heap
1.559 - _ns._in_edge = _candidates[0];
1.560 + _in_edge = _candidates[0];
1.561 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
1.562 _sort_func );
1.563 _curr_length = std::min(_head_length, _curr_length - 1);
1.564 return true;
1.565 }
1.566 +
1.567 }; //class AlteringListPivotRule
1.568
1.569 private:
1.570 @@ -519,66 +538,87 @@
1.571
1.572 private:
1.573
1.574 - // The directed graph the algorithm runs on
1.575 - SGraph _graph;
1.576 // The original graph
1.577 - const Graph &_graph_ref;
1.578 + const Graph &_orig_graph;
1.579 // The original lower bound map
1.580 - const LowerMap *_lower;
1.581 - // The capacity map
1.582 - SCapacityMap _capacity;
1.583 - // The cost map
1.584 - SCostMap _cost;
1.585 - // The supply map
1.586 - SSupplyMap _supply;
1.587 - bool _valid_supply;
1.588 + const LowerMap *_orig_lower;
1.589 + // The original capacity map
1.590 + const CapacityMap &_orig_cap;
1.591 + // The original cost map
1.592 + const CostMap &_orig_cost;
1.593 + // The original supply map
1.594 + const SupplyMap *_orig_supply;
1.595 + // The source node (if no supply map was given)
1.596 + Node _orig_source;
1.597 + // The target node (if no supply map was given)
1.598 + Node _orig_target;
1.599 + // The flow value (if no supply map was given)
1.600 + Capacity _orig_flow_value;
1.601
1.602 - // Edge map of the current flow
1.603 - SCapacityMap _flow;
1.604 - // Node map of the current potentials
1.605 - SPotentialMap _potential;
1.606 + // The flow result map
1.607 + FlowMap *_flow_result;
1.608 + // The potential result map
1.609 + PotentialMap *_potential_result;
1.610 + // Indicate if the flow result map is local
1.611 + bool _local_flow;
1.612 + // Indicate if the potential result map is local
1.613 + bool _local_potential;
1.614
1.615 - // The depth node map of the spanning tree structure
1.616 - IntNodeMap _depth;
1.617 - // The parent node map of the spanning tree structure
1.618 - NodeNodeMap _parent;
1.619 - // The pred_edge node map of the spanning tree structure
1.620 - EdgeNodeMap _pred_edge;
1.621 - // The thread node map of the spanning tree structure
1.622 - NodeNodeMap _thread;
1.623 - // The forward node map of the spanning tree structure
1.624 - BoolNodeMap _forward;
1.625 - // The state edge map
1.626 - IntEdgeMap _state;
1.627 - // The artificial root node of the spanning tree structure
1.628 - Node _root;
1.629 + // The edge references
1.630 + EdgeVector _edge;
1.631 + // The node references
1.632 + NodeVector _node;
1.633 + // The node ids
1.634 + IntNodeMap _node_id;
1.635 + // The source nodes of the edges
1.636 + IntVector _source;
1.637 + // The target nodess of the edges
1.638 + IntVector _target;
1.639
1.640 - // The reduced cost map
1.641 - ReducedCostMap _red_cost;
1.642 + // The (modified) capacity vector
1.643 + CapacityVector _cap;
1.644 + // The cost vector
1.645 + CostVector _cost;
1.646 + // The (modified) supply vector
1.647 + CostVector _supply;
1.648 + // The current flow vector
1.649 + CapacityVector _flow;
1.650 + // The current potential vector
1.651 + CostVector _pi;
1.652
1.653 - // The non-artifical edges
1.654 - EdgeVector _edges;
1.655 + // The number of nodes in the original graph
1.656 + int _node_num;
1.657 + // The number of edges in the original graph
1.658 + int _edge_num;
1.659
1.660 - // Members for handling the original graph
1.661 - FlowMap *_flow_result;
1.662 - PotentialMap *_potential_result;
1.663 - bool _local_flow;
1.664 - bool _local_potential;
1.665 - NodeRefMap _node_ref;
1.666 - EdgeRefMap _edge_ref;
1.667 + // The depth vector of the spanning tree structure
1.668 + IntVector _depth;
1.669 + // The parent vector of the spanning tree structure
1.670 + IntVector _parent;
1.671 + // The pred_edge vector of the spanning tree structure
1.672 + IntVector _pred;
1.673 + // The thread vector of the spanning tree structure
1.674 + IntVector _thread;
1.675 + // The forward vector of the spanning tree structure
1.676 + BoolVector _forward;
1.677 + // The state vector
1.678 + IntVector _state;
1.679 + // The root node
1.680 + int _root;
1.681
1.682 // The entering edge of the current pivot iteration
1.683 - Edge _in_edge;
1.684 + int _in_edge;
1.685
1.686 // Temporary nodes used in the current pivot iteration
1.687 - Node join, u_in, v_in, u_out, v_out;
1.688 - Node right, first, second, last;
1.689 - Node stem, par_stem, new_stem;
1.690 + int join, u_in, v_in, u_out, v_out;
1.691 + int right, first, second, last;
1.692 + int stem, par_stem, new_stem;
1.693 +
1.694 // The maximum augment amount along the found cycle in the current
1.695 // pivot iteration
1.696 Capacity delta;
1.697
1.698 - public :
1.699 + public:
1.700
1.701 /// \brief General constructor (with lower bounds).
1.702 ///
1.703 @@ -594,41 +634,12 @@
1.704 const CapacityMap &capacity,
1.705 const CostMap &cost,
1.706 const SupplyMap &supply ) :
1.707 - _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
1.708 - _cost(_graph), _supply(_graph), _flow(_graph),
1.709 - _potential(_graph), _depth(_graph), _parent(_graph),
1.710 - _pred_edge(_graph), _thread(_graph), _forward(_graph),
1.711 - _state(_graph), _red_cost(_graph, _cost, _potential),
1.712 + _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
1.713 + _orig_cost(cost), _orig_supply(&supply),
1.714 _flow_result(NULL), _potential_result(NULL),
1.715 _local_flow(false), _local_potential(false),
1.716 - _node_ref(graph), _edge_ref(graph)
1.717 - {
1.718 - // Check the sum of supply values
1.719 - Supply sum = 0;
1.720 - for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1.721 - sum += supply[n];
1.722 - if (!(_valid_supply = sum == 0)) return;
1.723 -
1.724 - // Copy _graph_ref to _graph
1.725 - _graph.reserveNode(countNodes(_graph_ref) + 1);
1.726 - _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
1.727 - copyGraph(_graph, _graph_ref)
1.728 - .edgeMap(_capacity, capacity)
1.729 - .edgeMap(_cost, cost)
1.730 - .nodeMap(_supply, supply)
1.731 - .nodeRef(_node_ref)
1.732 - .edgeRef(_edge_ref)
1.733 - .run();
1.734 -
1.735 - // Remove non-zero lower bounds
1.736 - for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
1.737 - if (lower[e] != 0) {
1.738 - _capacity[_edge_ref[e]] = capacity[e] - lower[e];
1.739 - _supply[_node_ref[_graph_ref.source(e)]] -= lower[e];
1.740 - _supply[_node_ref[_graph_ref.target(e)]] += lower[e];
1.741 - }
1.742 - }
1.743 - }
1.744 + _node_id(graph)
1.745 + {}
1.746
1.747 /// \brief General constructor (without lower bounds).
1.748 ///
1.749 @@ -642,32 +653,12 @@
1.750 const CapacityMap &capacity,
1.751 const CostMap &cost,
1.752 const SupplyMap &supply ) :
1.753 - _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
1.754 - _cost(_graph), _supply(_graph), _flow(_graph),
1.755 - _potential(_graph), _depth(_graph), _parent(_graph),
1.756 - _pred_edge(_graph), _thread(_graph), _forward(_graph),
1.757 - _state(_graph), _red_cost(_graph, _cost, _potential),
1.758 + _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
1.759 + _orig_cost(cost), _orig_supply(&supply),
1.760 _flow_result(NULL), _potential_result(NULL),
1.761 _local_flow(false), _local_potential(false),
1.762 - _node_ref(graph), _edge_ref(graph)
1.763 - {
1.764 - // Check the sum of supply values
1.765 - Supply sum = 0;
1.766 - for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1.767 - sum += supply[n];
1.768 - if (!(_valid_supply = sum == 0)) return;
1.769 -
1.770 - // Copy _graph_ref to _graph
1.771 - _graph.reserveNode(countNodes(_graph_ref) + 1);
1.772 - _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
1.773 - copyGraph(_graph, _graph_ref)
1.774 - .edgeMap(_capacity, capacity)
1.775 - .edgeMap(_cost, cost)
1.776 - .nodeMap(_supply, supply)
1.777 - .nodeRef(_node_ref)
1.778 - .edgeRef(_edge_ref)
1.779 - .run();
1.780 - }
1.781 + _node_id(graph)
1.782 + {}
1.783
1.784 /// \brief Simple constructor (with lower bounds).
1.785 ///
1.786 @@ -685,40 +676,15 @@
1.787 const LowerMap &lower,
1.788 const CapacityMap &capacity,
1.789 const CostMap &cost,
1.790 - typename Graph::Node s,
1.791 - typename Graph::Node t,
1.792 - typename SupplyMap::Value flow_value ) :
1.793 - _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
1.794 - _cost(_graph), _supply(_graph, 0), _flow(_graph),
1.795 - _potential(_graph), _depth(_graph), _parent(_graph),
1.796 - _pred_edge(_graph), _thread(_graph), _forward(_graph),
1.797 - _state(_graph), _red_cost(_graph, _cost, _potential),
1.798 + Node s, Node t,
1.799 + Capacity flow_value ) :
1.800 + _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
1.801 + _orig_cost(cost), _orig_supply(NULL),
1.802 + _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
1.803 _flow_result(NULL), _potential_result(NULL),
1.804 _local_flow(false), _local_potential(false),
1.805 - _node_ref(graph), _edge_ref(graph)
1.806 - {
1.807 - // Copy _graph_ref to graph
1.808 - _graph.reserveNode(countNodes(_graph_ref) + 1);
1.809 - _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
1.810 - copyGraph(_graph, _graph_ref)
1.811 - .edgeMap(_capacity, capacity)
1.812 - .edgeMap(_cost, cost)
1.813 - .nodeRef(_node_ref)
1.814 - .edgeRef(_edge_ref)
1.815 - .run();
1.816 -
1.817 - // Remove non-zero lower bounds
1.818 - _supply[_node_ref[s]] = flow_value;
1.819 - _supply[_node_ref[t]] = -flow_value;
1.820 - for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
1.821 - if (lower[e] != 0) {
1.822 - _capacity[_edge_ref[e]] = capacity[e] - lower[e];
1.823 - _supply[_node_ref[_graph_ref.source(e)]] -= lower[e];
1.824 - _supply[_node_ref[_graph_ref.target(e)]] += lower[e];
1.825 - }
1.826 - }
1.827 - _valid_supply = true;
1.828 - }
1.829 + _node_id(graph)
1.830 + {}
1.831
1.832 /// \brief Simple constructor (without lower bounds).
1.833 ///
1.834 @@ -734,31 +700,15 @@
1.835 NetworkSimplex( const Graph &graph,
1.836 const CapacityMap &capacity,
1.837 const CostMap &cost,
1.838 - typename Graph::Node s,
1.839 - typename Graph::Node t,
1.840 - typename SupplyMap::Value flow_value ) :
1.841 - _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
1.842 - _cost(_graph), _supply(_graph, 0), _flow(_graph),
1.843 - _potential(_graph), _depth(_graph), _parent(_graph),
1.844 - _pred_edge(_graph), _thread(_graph), _forward(_graph),
1.845 - _state(_graph), _red_cost(_graph, _cost, _potential),
1.846 + Node s, Node t,
1.847 + Capacity flow_value ) :
1.848 + _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
1.849 + _orig_cost(cost), _orig_supply(NULL),
1.850 + _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
1.851 _flow_result(NULL), _potential_result(NULL),
1.852 _local_flow(false), _local_potential(false),
1.853 - _node_ref(graph), _edge_ref(graph)
1.854 - {
1.855 - // Copy _graph_ref to graph
1.856 - _graph.reserveNode(countNodes(_graph_ref) + 1);
1.857 - _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
1.858 - copyGraph(_graph, _graph_ref)
1.859 - .edgeMap(_capacity, capacity)
1.860 - .edgeMap(_cost, cost)
1.861 - .nodeRef(_node_ref)
1.862 - .edgeRef(_edge_ref)
1.863 - .run();
1.864 - _supply[_node_ref[s]] = flow_value;
1.865 - _supply[_node_ref[t]] = -flow_value;
1.866 - _valid_supply = true;
1.867 - }
1.868 + _node_id(graph)
1.869 + {}
1.870
1.871 /// Destructor.
1.872 ~NetworkSimplex() {
1.873 @@ -894,8 +844,8 @@
1.874 /// \pre \ref run() must be called before using this function.
1.875 Cost totalCost() const {
1.876 Cost c = 0;
1.877 - for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1.878 - c += (*_flow_result)[e] * _cost[_edge_ref[e]];
1.879 + for (EdgeIt e(_orig_graph); e != INVALID; ++e)
1.880 + c += (*_flow_result)[e] * _orig_cost[e];
1.881 return c;
1.882 }
1.883
1.884 @@ -903,153 +853,216 @@
1.885
1.886 private:
1.887
1.888 - // Extend the underlying graph and initialize all the node and edge maps
1.889 + // Initialize internal data structures
1.890 bool init() {
1.891 - if (!_valid_supply) return false;
1.892 -
1.893 // Initialize result maps
1.894 if (!_flow_result) {
1.895 - _flow_result = new FlowMap(_graph_ref);
1.896 + _flow_result = new FlowMap(_orig_graph);
1.897 _local_flow = true;
1.898 }
1.899 if (!_potential_result) {
1.900 - _potential_result = new PotentialMap(_graph_ref);
1.901 + _potential_result = new PotentialMap(_orig_graph);
1.902 _local_potential = true;
1.903 }
1.904 -
1.905 - // Initialize the edge vector and the edge maps
1.906 - _edges.reserve(countEdges(_graph));
1.907 - for (EdgeIt e(_graph); e != INVALID; ++e) {
1.908 - _edges.push_back(e);
1.909 - _flow[e] = 0;
1.910 - _state[e] = STATE_LOWER;
1.911 +
1.912 + // Initialize vectors
1.913 + _node_num = countNodes(_orig_graph);
1.914 + _edge_num = countEdges(_orig_graph);
1.915 + int all_node_num = _node_num + 1;
1.916 + int all_edge_num = _edge_num + _node_num;
1.917 +
1.918 + _edge.resize(_edge_num);
1.919 + _node.reserve(_node_num);
1.920 + _source.resize(all_edge_num);
1.921 + _target.resize(all_edge_num);
1.922 +
1.923 + _cap.resize(all_edge_num);
1.924 + _cost.resize(all_edge_num);
1.925 + _supply.resize(all_node_num);
1.926 + _flow.resize(all_edge_num, 0);
1.927 + _pi.resize(all_node_num, 0);
1.928 +
1.929 + _depth.resize(all_node_num);
1.930 + _parent.resize(all_node_num);
1.931 + _pred.resize(all_node_num);
1.932 + _thread.resize(all_node_num);
1.933 + _forward.resize(all_node_num);
1.934 + _state.resize(all_edge_num, STATE_LOWER);
1.935 +
1.936 + // Initialize node related data
1.937 + bool valid_supply = true;
1.938 + if (_orig_supply) {
1.939 + Supply sum = 0;
1.940 + int i = 0;
1.941 + for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
1.942 + _node.push_back(n);
1.943 + _node_id[n] = i;
1.944 + _supply[i] = (*_orig_supply)[n];
1.945 + sum += _supply[i];
1.946 + }
1.947 + valid_supply = (sum == 0);
1.948 + } else {
1.949 + int i = 0;
1.950 + for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
1.951 + _node.push_back(n);
1.952 + _node_id[n] = i;
1.953 + _supply[i] = 0;
1.954 + }
1.955 + _supply[_node_id[_orig_source]] = _orig_flow_value;
1.956 + _supply[_node_id[_orig_target]] = -_orig_flow_value;
1.957 + }
1.958 + if (!valid_supply) return false;
1.959 +
1.960 + // Set data for the artificial root node
1.961 + _root = _node_num;
1.962 + _depth[_root] = 0;
1.963 + _parent[_root] = -1;
1.964 + _pred[_root] = -1;
1.965 + _thread[_root] = 0;
1.966 + _supply[_root] = 0;
1.967 + _pi[_root] = 0;
1.968 +
1.969 + // Store the edges in a mixed order
1.970 + int k = std::max(int(sqrt(_edge_num)), 10);
1.971 + int i = 0;
1.972 + for (EdgeIt e(_orig_graph); e != INVALID; ++e) {
1.973 + _edge[i] = e;
1.974 + if ((i += k) >= _edge_num) i = (i % k) + 1;
1.975 }
1.976
1.977 - // Add an artificial root node to the graph
1.978 - _root = _graph.addNode();
1.979 - _parent[_root] = INVALID;
1.980 - _pred_edge[_root] = INVALID;
1.981 - _depth[_root] = 0;
1.982 - _supply[_root] = 0;
1.983 - _potential[_root] = 0;
1.984 + // Initialize edge maps
1.985 + for (int i = 0; i != _edge_num; ++i) {
1.986 + Edge e = _edge[i];
1.987 + _source[i] = _node_id[_orig_graph.source(e)];
1.988 + _target[i] = _node_id[_orig_graph.target(e)];
1.989 + _cost[i] = _orig_cost[e];
1.990 + _cap[i] = _orig_cap[e];
1.991 + }
1.992
1.993 - // Add artificial edges to the graph and initialize the node maps
1.994 - // of the spanning tree data structure
1.995 - Node last = _root;
1.996 - Edge e;
1.997 + // Remove non-zero lower bounds
1.998 + if (_orig_lower) {
1.999 + for (int i = 0; i != _edge_num; ++i) {
1.1000 + Capacity c = (*_orig_lower)[_edge[i]];
1.1001 + if (c != 0) {
1.1002 + _cap[i] -= c;
1.1003 + _supply[_source[i]] -= c;
1.1004 + _supply[_target[i]] += c;
1.1005 + }
1.1006 + }
1.1007 + }
1.1008 +
1.1009 + // Add artificial edges and initialize the spanning tree data structure
1.1010 Cost max_cost = std::numeric_limits<Cost>::max() / 4;
1.1011 - for (NodeIt u(_graph); u != INVALID; ++u) {
1.1012 - if (u == _root) continue;
1.1013 - _thread[last] = u;
1.1014 - last = u;
1.1015 + Capacity max_cap = std::numeric_limits<Capacity>::max();
1.1016 + for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
1.1017 + _thread[u] = u + 1;
1.1018 + _depth[u] = 1;
1.1019 _parent[u] = _root;
1.1020 - _depth[u] = 1;
1.1021 + _pred[u] = e;
1.1022 if (_supply[u] >= 0) {
1.1023 - e = _graph.addEdge(u, _root);
1.1024 _flow[e] = _supply[u];
1.1025 _forward[u] = true;
1.1026 - _potential[u] = -max_cost;
1.1027 + _pi[u] = -max_cost;
1.1028 } else {
1.1029 - e = _graph.addEdge(_root, u);
1.1030 _flow[e] = -_supply[u];
1.1031 _forward[u] = false;
1.1032 - _potential[u] = max_cost;
1.1033 + _pi[u] = max_cost;
1.1034 }
1.1035 _cost[e] = max_cost;
1.1036 - _capacity[e] = std::numeric_limits<Capacity>::max();
1.1037 + _cap[e] = max_cap;
1.1038 _state[e] = STATE_TREE;
1.1039 - _pred_edge[u] = e;
1.1040 }
1.1041 - _thread[last] = _root;
1.1042
1.1043 return true;
1.1044 }
1.1045
1.1046 // Find the join node
1.1047 void findJoinNode() {
1.1048 - Node u = _graph.source(_in_edge);
1.1049 - Node v = _graph.target(_in_edge);
1.1050 + int u = _source[_in_edge];
1.1051 + int v = _target[_in_edge];
1.1052 + while (_depth[u] > _depth[v]) u = _parent[u];
1.1053 + while (_depth[v] > _depth[u]) v = _parent[v];
1.1054 while (u != v) {
1.1055 - if (_depth[u] == _depth[v]) {
1.1056 - u = _parent[u];
1.1057 - v = _parent[v];
1.1058 - }
1.1059 - else if (_depth[u] > _depth[v]) u = _parent[u];
1.1060 - else v = _parent[v];
1.1061 + u = _parent[u];
1.1062 + v = _parent[v];
1.1063 }
1.1064 join = u;
1.1065 }
1.1066
1.1067 - // Find the leaving edge of the cycle and returns true if the
1.1068 + // Find the leaving edge of the cycle and returns true if the
1.1069 // leaving edge is not the same as the entering edge
1.1070 bool findLeavingEdge() {
1.1071 // Initialize first and second nodes according to the direction
1.1072 // of the cycle
1.1073 if (_state[_in_edge] == STATE_LOWER) {
1.1074 - first = _graph.source(_in_edge);
1.1075 - second = _graph.target(_in_edge);
1.1076 + first = _source[_in_edge];
1.1077 + second = _target[_in_edge];
1.1078 } else {
1.1079 - first = _graph.target(_in_edge);
1.1080 - second = _graph.source(_in_edge);
1.1081 + first = _target[_in_edge];
1.1082 + second = _source[_in_edge];
1.1083 }
1.1084 - delta = _capacity[_in_edge];
1.1085 - bool result = false;
1.1086 + delta = _cap[_in_edge];
1.1087 + int result = 0;
1.1088 Capacity d;
1.1089 - Edge e;
1.1090 + int e;
1.1091
1.1092 // Search the cycle along the path form the first node to the root
1.1093 - for (Node u = first; u != join; u = _parent[u]) {
1.1094 - e = _pred_edge[u];
1.1095 - d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
1.1096 + for (int u = first; u != join; u = _parent[u]) {
1.1097 + e = _pred[u];
1.1098 + d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
1.1099 if (d < delta) {
1.1100 delta = d;
1.1101 u_out = u;
1.1102 - u_in = first;
1.1103 - v_in = second;
1.1104 - result = true;
1.1105 + result = 1;
1.1106 }
1.1107 }
1.1108 // Search the cycle along the path form the second node to the root
1.1109 - for (Node u = second; u != join; u = _parent[u]) {
1.1110 - e = _pred_edge[u];
1.1111 - d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
1.1112 + for (int u = second; u != join; u = _parent[u]) {
1.1113 + e = _pred[u];
1.1114 + d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1.1115 if (d <= delta) {
1.1116 delta = d;
1.1117 u_out = u;
1.1118 - u_in = second;
1.1119 - v_in = first;
1.1120 - result = true;
1.1121 + result = 2;
1.1122 }
1.1123 }
1.1124 - return result;
1.1125 +
1.1126 + if (result == 1) {
1.1127 + u_in = first;
1.1128 + v_in = second;
1.1129 + } else {
1.1130 + u_in = second;
1.1131 + v_in = first;
1.1132 + }
1.1133 + return result != 0;
1.1134 }
1.1135
1.1136 - // Change _flow and _state edge maps
1.1137 - void changeFlows(bool change) {
1.1138 + // Change _flow and _state vectors
1.1139 + void changeFlow(bool change) {
1.1140 // Augment along the cycle
1.1141 if (delta > 0) {
1.1142 Capacity val = _state[_in_edge] * delta;
1.1143 _flow[_in_edge] += val;
1.1144 - for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
1.1145 - _flow[_pred_edge[u]] += _forward[u] ? -val : val;
1.1146 + for (int u = _source[_in_edge]; u != join; u = _parent[u]) {
1.1147 + _flow[_pred[u]] += _forward[u] ? -val : val;
1.1148 }
1.1149 - for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
1.1150 - _flow[_pred_edge[u]] += _forward[u] ? val : -val;
1.1151 + for (int u = _target[_in_edge]; u != join; u = _parent[u]) {
1.1152 + _flow[_pred[u]] += _forward[u] ? val : -val;
1.1153 }
1.1154 }
1.1155 // Update the state of the entering and leaving edges
1.1156 if (change) {
1.1157 _state[_in_edge] = STATE_TREE;
1.1158 - _state[_pred_edge[u_out]] =
1.1159 - (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1.1160 + _state[_pred[u_out]] =
1.1161 + (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
1.1162 } else {
1.1163 _state[_in_edge] = -_state[_in_edge];
1.1164 }
1.1165 }
1.1166
1.1167 - // Update _thread and _parent node maps
1.1168 + // Update _thread and _parent vectors
1.1169 void updateThreadParent() {
1.1170 - Node u;
1.1171 + int u;
1.1172 v_out = _parent[u_out];
1.1173
1.1174 // Handle the case when join and v_out coincide
1.1175 @@ -1105,30 +1118,30 @@
1.1176 }
1.1177 }
1.1178
1.1179 - // Update _pred_edge and _forward node maps.
1.1180 + // Update _pred and _forward vectors
1.1181 void updatePredEdge() {
1.1182 - Node u = u_out, v;
1.1183 + int u = u_out, v;
1.1184 while (u != u_in) {
1.1185 v = _parent[u];
1.1186 - _pred_edge[u] = _pred_edge[v];
1.1187 + _pred[u] = _pred[v];
1.1188 _forward[u] = !_forward[v];
1.1189 u = v;
1.1190 }
1.1191 - _pred_edge[u_in] = _in_edge;
1.1192 - _forward[u_in] = (u_in == _graph.source(_in_edge));
1.1193 + _pred[u_in] = _in_edge;
1.1194 + _forward[u_in] = (u_in == _source[_in_edge]);
1.1195 }
1.1196
1.1197 - // Update _depth and _potential node maps
1.1198 + // Update _depth and _potential vectors
1.1199 void updateDepthPotential() {
1.1200 _depth[u_in] = _depth[v_in] + 1;
1.1201 Cost sigma = _forward[u_in] ?
1.1202 - _potential[v_in] - _potential[u_in] - _cost[_pred_edge[u_in]] :
1.1203 - _potential[v_in] - _potential[u_in] + _cost[_pred_edge[u_in]];
1.1204 - _potential[u_in] += sigma;
1.1205 - for(Node u = _thread[u_in]; _parent[u] != INVALID; u = _thread[u]) {
1.1206 + _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1.1207 + _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1.1208 + _pi[u_in] += sigma;
1.1209 + for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
1.1210 _depth[u] = _depth[_parent[u]] + 1;
1.1211 if (_depth[u] <= _depth[u_in]) break;
1.1212 - _potential[u] += sigma;
1.1213 + _pi[u] += sigma;
1.1214 }
1.1215 }
1.1216
1.1217 @@ -1152,13 +1165,13 @@
1.1218
1.1219 template<class PivotRuleImplementation>
1.1220 bool start() {
1.1221 - PivotRuleImplementation pivot(*this, _edges);
1.1222 + PivotRuleImplementation pivot(*this);
1.1223
1.1224 // Execute the network simplex algorithm
1.1225 while (pivot.findEnteringEdge()) {
1.1226 findJoinNode();
1.1227 bool change = findLeavingEdge();
1.1228 - changeFlows(change);
1.1229 + changeFlow(change);
1.1230 if (change) {
1.1231 updateThreadParent();
1.1232 updatePredEdge();
1.1233 @@ -1167,22 +1180,25 @@
1.1234 }
1.1235
1.1236 // Check if the flow amount equals zero on all the artificial edges
1.1237 - for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
1.1238 + for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
1.1239 if (_flow[e] > 0) return false;
1.1240 - for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
1.1241 - if (_flow[e] > 0) return false;
1.1242 + }
1.1243
1.1244 // Copy flow values to _flow_result
1.1245 - if (_lower) {
1.1246 - for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1.1247 - (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
1.1248 + if (_orig_lower) {
1.1249 + for (int i = 0; i != _edge_num; ++i) {
1.1250 + Edge e = _edge[i];
1.1251 + (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
1.1252 + }
1.1253 } else {
1.1254 - for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
1.1255 - (*_flow_result)[e] = _flow[_edge_ref[e]];
1.1256 + for (int i = 0; i != _edge_num; ++i) {
1.1257 + (*_flow_result)[_edge[i]] = _flow[i];
1.1258 + }
1.1259 }
1.1260 // Copy potential values to _potential_result
1.1261 - for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
1.1262 - (*_potential_result)[n] = _potential[_node_ref[n]];
1.1263 + for (int i = 0; i != _node_num; ++i) {
1.1264 + (*_potential_result)[_node[i]] = _pi[i];
1.1265 + }
1.1266
1.1267 return true;
1.1268 }