[2440] | 1 | /* -*- C++ -*- |
---|
| 2 | * |
---|
| 3 | * This file is a part of LEMON, a generic C++ optimization library |
---|
| 4 | * |
---|
[2553] | 5 | * Copyright (C) 2003-2008 |
---|
[2440] | 6 | * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
---|
| 7 | * (Egervary Research Group on Combinatorial Optimization, EGRES). |
---|
| 8 | * |
---|
| 9 | * Permission to use, modify and distribute this software is granted |
---|
| 10 | * provided that this copyright notice appears in all copies. For |
---|
| 11 | * precise terms see the accompanying LICENSE file. |
---|
| 12 | * |
---|
| 13 | * This software is provided "AS IS" with no warranty of any kind, |
---|
| 14 | * express or implied, and with no claim as to its suitability for any |
---|
| 15 | * purpose. |
---|
| 16 | * |
---|
| 17 | */ |
---|
| 18 | |
---|
| 19 | #ifndef LEMON_NETWORK_SIMPLEX_H |
---|
| 20 | #define LEMON_NETWORK_SIMPLEX_H |
---|
| 21 | |
---|
| 22 | /// \ingroup min_cost_flow |
---|
| 23 | /// |
---|
| 24 | /// \file |
---|
[2575] | 25 | /// \brief Network simplex algorithm for finding a minimum cost flow. |
---|
[2440] | 26 | |
---|
[2575] | 27 | #include <vector> |
---|
[2440] | 28 | #include <limits> |
---|
[2575] | 29 | |
---|
[2509] | 30 | #include <lemon/graph_adaptor.h> |
---|
| 31 | #include <lemon/graph_utils.h> |
---|
[2440] | 32 | #include <lemon/smart_graph.h> |
---|
[2575] | 33 | #include <lemon/math.h> |
---|
[2440] | 34 | |
---|
| 35 | namespace lemon { |
---|
| 36 | |
---|
| 37 | /// \addtogroup min_cost_flow |
---|
| 38 | /// @{ |
---|
| 39 | |
---|
| 40 | /// \brief Implementation of the network simplex algorithm for |
---|
| 41 | /// finding a minimum cost flow. |
---|
| 42 | /// |
---|
[2556] | 43 | /// \ref NetworkSimplex implements the network simplex algorithm for |
---|
| 44 | /// finding a minimum cost flow. |
---|
[2440] | 45 | /// |
---|
[2575] | 46 | /// \tparam Graph The directed graph type the algorithm runs on. |
---|
| 47 | /// \tparam LowerMap The type of the lower bound map. |
---|
| 48 | /// \tparam CapacityMap The type of the capacity (upper bound) map. |
---|
| 49 | /// \tparam CostMap The type of the cost (length) map. |
---|
| 50 | /// \tparam SupplyMap The type of the supply map. |
---|
[2440] | 51 | /// |
---|
| 52 | /// \warning |
---|
[2575] | 53 | /// - Edge capacities and costs should be \e non-negative \e integers. |
---|
| 54 | /// - Supply values should be \e signed \e integers. |
---|
| 55 | /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value. |
---|
| 56 | /// - \c CapacityMap::Value and \c SupplyMap::Value must be |
---|
| 57 | /// convertible to each other. |
---|
| 58 | /// - All value types must be convertible to \c CostMap::Value, which |
---|
| 59 | /// must be signed type. |
---|
| 60 | /// |
---|
| 61 | /// \note \ref NetworkSimplex provides six different pivot rule |
---|
| 62 | /// implementations that significantly affect the efficiency of the |
---|
| 63 | /// algorithm. |
---|
| 64 | /// By default a combined pivot rule is used, which is the fastest |
---|
| 65 | /// implementation according to our benchmark tests. |
---|
| 66 | /// Another pivot rule can be selected using \ref run() function |
---|
| 67 | /// with the proper parameter. |
---|
[2440] | 68 | /// |
---|
| 69 | /// \author Peter Kovacs |
---|
| 70 | |
---|
[2533] | 71 | template < typename Graph, |
---|
| 72 | typename LowerMap = typename Graph::template EdgeMap<int>, |
---|
[2575] | 73 | typename CapacityMap = typename Graph::template EdgeMap<int>, |
---|
[2533] | 74 | typename CostMap = typename Graph::template EdgeMap<int>, |
---|
[2575] | 75 | typename SupplyMap = typename Graph::template NodeMap<int> > |
---|
[2440] | 76 | class NetworkSimplex |
---|
| 77 | { |
---|
| 78 | typedef typename CapacityMap::Value Capacity; |
---|
| 79 | typedef typename CostMap::Value Cost; |
---|
| 80 | typedef typename SupplyMap::Value Supply; |
---|
| 81 | |
---|
| 82 | typedef SmartGraph SGraph; |
---|
[2556] | 83 | GRAPH_TYPEDEFS(typename SGraph); |
---|
[2440] | 84 | |
---|
| 85 | typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap; |
---|
| 86 | typedef typename SGraph::template EdgeMap<Cost> SCostMap; |
---|
| 87 | typedef typename SGraph::template NodeMap<Supply> SSupplyMap; |
---|
| 88 | typedef typename SGraph::template NodeMap<Cost> SPotentialMap; |
---|
| 89 | |
---|
| 90 | typedef typename SGraph::template NodeMap<int> IntNodeMap; |
---|
| 91 | typedef typename SGraph::template NodeMap<bool> BoolNodeMap; |
---|
| 92 | typedef typename SGraph::template NodeMap<Node> NodeNodeMap; |
---|
| 93 | typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap; |
---|
| 94 | typedef typename SGraph::template EdgeMap<int> IntEdgeMap; |
---|
| 95 | |
---|
| 96 | typedef typename Graph::template NodeMap<Node> NodeRefMap; |
---|
| 97 | typedef typename Graph::template EdgeMap<Edge> EdgeRefMap; |
---|
| 98 | |
---|
| 99 | public: |
---|
| 100 | |
---|
[2556] | 101 | /// The type of the flow map. |
---|
[2440] | 102 | typedef typename Graph::template EdgeMap<Capacity> FlowMap; |
---|
[2556] | 103 | /// The type of the potential map. |
---|
[2440] | 104 | typedef typename Graph::template NodeMap<Cost> PotentialMap; |
---|
| 105 | |
---|
[2575] | 106 | public: |
---|
[2440] | 107 | |
---|
[2575] | 108 | /// Enum type to select the pivot rule used by \ref run(). |
---|
| 109 | enum PivotRuleEnum { |
---|
| 110 | FIRST_ELIGIBLE_PIVOT, |
---|
| 111 | BEST_ELIGIBLE_PIVOT, |
---|
| 112 | BLOCK_SEARCH_PIVOT, |
---|
| 113 | LIMITED_SEARCH_PIVOT, |
---|
| 114 | CANDIDATE_LIST_PIVOT, |
---|
| 115 | COMBINED_PIVOT |
---|
| 116 | }; |
---|
| 117 | |
---|
| 118 | private: |
---|
| 119 | |
---|
| 120 | /// \brief Map adaptor class for handling reduced edge costs. |
---|
| 121 | /// |
---|
[2556] | 122 | /// Map adaptor class for handling reduced edge costs. |
---|
[2440] | 123 | class ReducedCostMap : public MapBase<Edge, Cost> |
---|
| 124 | { |
---|
| 125 | private: |
---|
| 126 | |
---|
[2575] | 127 | const SGraph &_gr; |
---|
| 128 | const SCostMap &_cost_map; |
---|
| 129 | const SPotentialMap &_pot_map; |
---|
[2440] | 130 | |
---|
| 131 | public: |
---|
| 132 | |
---|
[2575] | 133 | ///\e |
---|
| 134 | ReducedCostMap( const SGraph &gr, |
---|
| 135 | const SCostMap &cost_map, |
---|
| 136 | const SPotentialMap &pot_map ) : |
---|
| 137 | _gr(gr), _cost_map(cost_map), _pot_map(pm) {} |
---|
[2440] | 138 | |
---|
[2575] | 139 | ///\e |
---|
[2509] | 140 | Cost operator[](const Edge &e) const { |
---|
[2575] | 141 | return _cost_map[e] + _pot_map[_gr.source(e)] |
---|
| 142 | - _pot_map[_gr.target(e)]; |
---|
[2440] | 143 | } |
---|
| 144 | |
---|
| 145 | }; //class ReducedCostMap |
---|
| 146 | |
---|
[2575] | 147 | private: |
---|
[2440] | 148 | |
---|
[2575] | 149 | /// \brief Implementation of the "First Eligible" pivot rule for the |
---|
| 150 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 151 | /// |
---|
| 152 | /// This class implements the "First Eligible" pivot rule |
---|
| 153 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
| 154 | class FirstEligiblePivotRule |
---|
| 155 | { |
---|
| 156 | private: |
---|
[2440] | 157 | |
---|
[2575] | 158 | NetworkSimplex &_ns; |
---|
| 159 | EdgeIt _next_edge; |
---|
[2440] | 160 | |
---|
[2575] | 161 | public: |
---|
[2440] | 162 | |
---|
[2575] | 163 | /// Constructor. |
---|
| 164 | FirstEligiblePivotRule(NetworkSimplex &ns) : |
---|
| 165 | _ns(ns), _next_edge(ns._graph) {} |
---|
| 166 | |
---|
| 167 | /// Finds the next entering edge. |
---|
| 168 | bool findEnteringEdge() { |
---|
| 169 | for (EdgeIt e = _next_edge; e != INVALID; ++e) { |
---|
| 170 | if (_ns._state[e] * _ns._red_cost[e] < 0) { |
---|
| 171 | _ns._in_edge = e; |
---|
| 172 | _next_edge = ++e; |
---|
| 173 | return true; |
---|
| 174 | } |
---|
| 175 | } |
---|
| 176 | for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { |
---|
| 177 | if (_ns._state[e] * _ns._red_cost[e] < 0) { |
---|
| 178 | _ns._in_edge = e; |
---|
| 179 | _next_edge = ++e; |
---|
| 180 | return true; |
---|
| 181 | } |
---|
| 182 | } |
---|
| 183 | return false; |
---|
| 184 | } |
---|
| 185 | }; //class FirstEligiblePivotRule |
---|
| 186 | |
---|
| 187 | /// \brief Implementation of the "Best Eligible" pivot rule for the |
---|
| 188 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 189 | /// |
---|
| 190 | /// This class implements the "Best Eligible" pivot rule |
---|
| 191 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
| 192 | class BestEligiblePivotRule |
---|
| 193 | { |
---|
| 194 | private: |
---|
| 195 | |
---|
| 196 | NetworkSimplex &_ns; |
---|
| 197 | |
---|
| 198 | public: |
---|
| 199 | |
---|
| 200 | /// Constructor. |
---|
| 201 | BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {} |
---|
| 202 | |
---|
| 203 | /// Finds the next entering edge. |
---|
| 204 | bool findEnteringEdge() { |
---|
| 205 | Cost min = 0; |
---|
| 206 | for (EdgeIt e(_ns._graph); e != INVALID; ++e) { |
---|
| 207 | if (_ns._state[e] * _ns._red_cost[e] < min) { |
---|
| 208 | min = _ns._state[e] * _ns._red_cost[e]; |
---|
| 209 | _ns._in_edge = e; |
---|
| 210 | } |
---|
| 211 | } |
---|
| 212 | return min < 0; |
---|
| 213 | } |
---|
| 214 | }; //class BestEligiblePivotRule |
---|
| 215 | |
---|
| 216 | /// \brief Implementation of the "Block Search" pivot rule for the |
---|
| 217 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 218 | /// |
---|
| 219 | /// This class implements the "Block Search" pivot rule |
---|
| 220 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
| 221 | class BlockSearchPivotRule |
---|
| 222 | { |
---|
| 223 | private: |
---|
| 224 | |
---|
| 225 | NetworkSimplex &_ns; |
---|
| 226 | EdgeIt _next_edge, _min_edge; |
---|
| 227 | int _block_size; |
---|
| 228 | |
---|
| 229 | static const int MIN_BLOCK_SIZE = 10; |
---|
| 230 | |
---|
| 231 | public: |
---|
| 232 | |
---|
| 233 | /// Constructor. |
---|
| 234 | BlockSearchPivotRule(NetworkSimplex &ns) : |
---|
| 235 | _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) |
---|
| 236 | { |
---|
| 237 | _block_size = 2 * int(sqrt(countEdges(_ns._graph))); |
---|
| 238 | if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE; |
---|
| 239 | } |
---|
| 240 | |
---|
| 241 | /// Finds the next entering edge. |
---|
| 242 | bool findEnteringEdge() { |
---|
| 243 | Cost curr, min = 0; |
---|
| 244 | int cnt = 0; |
---|
| 245 | for (EdgeIt e = _next_edge; e != INVALID; ++e) { |
---|
| 246 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
---|
| 247 | min = curr; |
---|
| 248 | _min_edge = e; |
---|
| 249 | } |
---|
| 250 | if (++cnt == _block_size) { |
---|
| 251 | if (min < 0) break; |
---|
| 252 | cnt = 0; |
---|
| 253 | } |
---|
| 254 | } |
---|
| 255 | if (min == 0) { |
---|
| 256 | for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { |
---|
| 257 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
---|
| 258 | min = curr; |
---|
| 259 | _min_edge = e; |
---|
| 260 | } |
---|
| 261 | if (++cnt == _block_size) { |
---|
| 262 | if (min < 0) break; |
---|
| 263 | cnt = 0; |
---|
| 264 | } |
---|
| 265 | } |
---|
| 266 | } |
---|
| 267 | _ns._in_edge = _min_edge; |
---|
| 268 | _next_edge = ++_min_edge; |
---|
| 269 | return min < 0; |
---|
| 270 | } |
---|
| 271 | }; //class BlockSearchPivotRule |
---|
| 272 | |
---|
| 273 | /// \brief Implementation of the "Limited Search" pivot rule for the |
---|
| 274 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 275 | /// |
---|
| 276 | /// This class implements the "Limited Search" pivot rule |
---|
| 277 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
| 278 | class LimitedSearchPivotRule |
---|
| 279 | { |
---|
| 280 | private: |
---|
| 281 | |
---|
| 282 | NetworkSimplex &_ns; |
---|
| 283 | EdgeIt _next_edge, _min_edge; |
---|
| 284 | int _sample_size; |
---|
| 285 | |
---|
| 286 | static const int MIN_SAMPLE_SIZE = 10; |
---|
| 287 | static const double SAMPLE_SIZE_FACTOR = 0.0015; |
---|
| 288 | |
---|
| 289 | public: |
---|
| 290 | |
---|
| 291 | /// Constructor. |
---|
| 292 | LimitedSearchPivotRule(NetworkSimplex &ns) : |
---|
| 293 | _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph) |
---|
| 294 | { |
---|
| 295 | _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph)); |
---|
| 296 | if (_sample_size < MIN_SAMPLE_SIZE) |
---|
| 297 | _sample_size = MIN_SAMPLE_SIZE; |
---|
| 298 | } |
---|
| 299 | |
---|
| 300 | /// Finds the next entering edge. |
---|
| 301 | bool findEnteringEdge() { |
---|
| 302 | Cost curr, min = 0; |
---|
| 303 | int cnt = 0; |
---|
| 304 | for (EdgeIt e = _next_edge; e != INVALID; ++e) { |
---|
| 305 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
---|
| 306 | min = curr; |
---|
| 307 | _min_edge = e; |
---|
| 308 | } |
---|
| 309 | if (curr < 0 && ++cnt == _sample_size) break; |
---|
| 310 | } |
---|
| 311 | if (min == 0) { |
---|
| 312 | for (EdgeIt e(_ns._graph); e != _next_edge; ++e) { |
---|
| 313 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
---|
| 314 | min = curr; |
---|
| 315 | _min_edge = e; |
---|
| 316 | } |
---|
| 317 | if (curr < 0 && ++cnt == _sample_size) break; |
---|
| 318 | } |
---|
| 319 | } |
---|
| 320 | _ns._in_edge = _min_edge; |
---|
| 321 | _next_edge = ++_min_edge; |
---|
| 322 | return min < 0; |
---|
| 323 | } |
---|
| 324 | }; //class LimitedSearchPivotRule |
---|
| 325 | |
---|
| 326 | /// \brief Implementation of the "Candidate List" pivot rule for the |
---|
| 327 | /// \ref NetworkSimplex "network simplex" algorithm. |
---|
| 328 | /// |
---|
| 329 | /// This class implements the "Candidate List" pivot rule |
---|
| 330 | /// for the \ref NetworkSimplex "network simplex" algorithm. |
---|
| 331 | class CandidateListPivotRule |
---|
| 332 | { |
---|
| 333 | private: |
---|
| 334 | |
---|
| 335 | NetworkSimplex &_ns; |
---|
| 336 | |
---|
| 337 | // The list of candidate edges. |
---|
| 338 | std::vector<Edge> _candidates; |
---|
| 339 | // The maximum length of the edge list. |
---|
| 340 | int _list_length; |
---|
| 341 | // The maximum number of minor iterations between two major |
---|
| 342 | // itarations. |
---|
| 343 | int _minor_limit; |
---|
| 344 | |
---|
| 345 | int _minor_count; |
---|
| 346 | EdgeIt _next_edge; |
---|
| 347 | |
---|
| 348 | static const double LIST_LENGTH_FACTOR = 0.002; |
---|
| 349 | static const double MINOR_LIMIT_FACTOR = 0.1; |
---|
| 350 | static const int MIN_LIST_LENGTH = 10; |
---|
| 351 | static const int MIN_MINOR_LIMIT = 2; |
---|
| 352 | |
---|
| 353 | public: |
---|
| 354 | |
---|
| 355 | /// Constructor. |
---|
| 356 | CandidateListPivotRule(NetworkSimplex &ns) : |
---|
| 357 | _ns(ns), _next_edge(ns._graph) |
---|
| 358 | { |
---|
| 359 | int edge_num = countEdges(_ns._graph); |
---|
| 360 | _minor_count = 0; |
---|
| 361 | _list_length = int(edge_num * LIST_LENGTH_FACTOR); |
---|
| 362 | if (_list_length < MIN_LIST_LENGTH) |
---|
| 363 | _list_length = MIN_LIST_LENGTH; |
---|
| 364 | _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR); |
---|
| 365 | if (_minor_limit < MIN_MINOR_LIMIT) |
---|
| 366 | _minor_limit = MIN_MINOR_LIMIT; |
---|
| 367 | } |
---|
| 368 | |
---|
| 369 | /// Finds the next entering edge. |
---|
| 370 | bool findEnteringEdge() { |
---|
| 371 | Cost min, curr; |
---|
| 372 | if (_minor_count < _minor_limit && _candidates.size() > 0) { |
---|
| 373 | // Minor iteration |
---|
| 374 | ++_minor_count; |
---|
| 375 | Edge e; |
---|
| 376 | min = 0; |
---|
| 377 | for (int i = 0; i < int(_candidates.size()); ++i) { |
---|
| 378 | e = _candidates[i]; |
---|
| 379 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) { |
---|
| 380 | min = curr; |
---|
| 381 | _ns._in_edge = e; |
---|
| 382 | } |
---|
| 383 | } |
---|
| 384 | if (min < 0) return true; |
---|
| 385 | } |
---|
| 386 | |
---|
| 387 | // Major iteration |
---|
| 388 | _candidates.clear(); |
---|
| 389 | EdgeIt e = _next_edge; |
---|
| 390 | min = 0; |
---|
| 391 | for ( ; e != INVALID; ++e) { |
---|
| 392 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { |
---|
| 393 | _candidates.push_back(e); |
---|
| 394 | if (curr < min) { |
---|
| 395 | min = curr; |
---|
| 396 | _ns._in_edge = e; |
---|
| 397 | } |
---|
| 398 | if (int(_candidates.size()) == _list_length) break; |
---|
| 399 | } |
---|
| 400 | } |
---|
| 401 | if (int(_candidates.size()) < _list_length) { |
---|
| 402 | for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) { |
---|
| 403 | if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) { |
---|
| 404 | _candidates.push_back(e); |
---|
| 405 | if (curr < min) { |
---|
| 406 | min = curr; |
---|
| 407 | _ns._in_edge = e; |
---|
| 408 | } |
---|
| 409 | if (int(_candidates.size()) == _list_length) break; |
---|
| 410 | } |
---|
| 411 | } |
---|
| 412 | } |
---|
| 413 | if (_candidates.size() == 0) return false; |
---|
| 414 | _minor_count = 1; |
---|
| 415 | _next_edge = ++e; |
---|
| 416 | return true; |
---|
| 417 | } |
---|
| 418 | }; //class CandidateListPivotRule |
---|
| 419 | |
---|
| 420 | private: |
---|
| 421 | |
---|
| 422 | // State constant for edges at their lower bounds. |
---|
| 423 | static const int STATE_LOWER = 1; |
---|
| 424 | // State constant for edges in the spanning tree. |
---|
| 425 | static const int STATE_TREE = 0; |
---|
| 426 | // State constant for edges at their upper bounds. |
---|
| 427 | static const int STATE_UPPER = -1; |
---|
| 428 | |
---|
| 429 | // Constant for the combined pivot rule. |
---|
| 430 | static const int COMBINED_PIVOT_MAX_DEG = 5; |
---|
| 431 | |
---|
| 432 | private: |
---|
| 433 | |
---|
| 434 | // The directed graph the algorithm runs on |
---|
| 435 | SGraph _graph; |
---|
| 436 | // The original graph |
---|
| 437 | const Graph &_graph_ref; |
---|
| 438 | // The original lower bound map |
---|
| 439 | const LowerMap *_lower; |
---|
| 440 | // The capacity map |
---|
| 441 | SCapacityMap _capacity; |
---|
| 442 | // The cost map |
---|
| 443 | SCostMap _cost; |
---|
| 444 | // The supply map |
---|
| 445 | SSupplyMap _supply; |
---|
| 446 | bool _valid_supply; |
---|
| 447 | |
---|
| 448 | // Edge map of the current flow |
---|
| 449 | SCapacityMap _flow; |
---|
| 450 | // Node map of the current potentials |
---|
| 451 | SPotentialMap _potential; |
---|
| 452 | |
---|
| 453 | // The depth node map of the spanning tree structure |
---|
| 454 | IntNodeMap _depth; |
---|
| 455 | // The parent node map of the spanning tree structure |
---|
| 456 | NodeNodeMap _parent; |
---|
| 457 | // The pred_edge node map of the spanning tree structure |
---|
| 458 | EdgeNodeMap _pred_edge; |
---|
| 459 | // The thread node map of the spanning tree structure |
---|
| 460 | NodeNodeMap _thread; |
---|
| 461 | // The forward node map of the spanning tree structure |
---|
| 462 | BoolNodeMap _forward; |
---|
| 463 | // The state edge map |
---|
| 464 | IntEdgeMap _state; |
---|
| 465 | // The root node of the starting spanning tree |
---|
| 466 | Node _root; |
---|
| 467 | |
---|
| 468 | // The reduced cost map |
---|
| 469 | ReducedCostMap _red_cost; |
---|
| 470 | |
---|
| 471 | // Members for handling the original graph |
---|
| 472 | FlowMap _flow_result; |
---|
| 473 | PotentialMap _potential_result; |
---|
| 474 | NodeRefMap _node_ref; |
---|
| 475 | EdgeRefMap _edge_ref; |
---|
[2440] | 476 | |
---|
[2556] | 477 | // The entering edge of the current pivot iteration. |
---|
[2575] | 478 | Edge _in_edge; |
---|
| 479 | |
---|
[2556] | 480 | // Temporary nodes used in the current pivot iteration. |
---|
| 481 | Node join, u_in, v_in, u_out, v_out; |
---|
| 482 | Node right, first, second, last; |
---|
| 483 | Node stem, par_stem, new_stem; |
---|
| 484 | // The maximum augment amount along the found cycle in the current |
---|
| 485 | // pivot iteration. |
---|
| 486 | Capacity delta; |
---|
[2440] | 487 | |
---|
| 488 | public : |
---|
| 489 | |
---|
| 490 | /// \brief General constructor of the class (with lower bounds). |
---|
| 491 | /// |
---|
| 492 | /// General constructor of the class (with lower bounds). |
---|
| 493 | /// |
---|
[2575] | 494 | /// \param graph The directed graph the algorithm runs on. |
---|
| 495 | /// \param lower The lower bounds of the edges. |
---|
| 496 | /// \param capacity The capacities (upper bounds) of the edges. |
---|
| 497 | /// \param cost The cost (length) values of the edges. |
---|
| 498 | /// \param supply The supply values of the nodes (signed). |
---|
| 499 | NetworkSimplex( const Graph &graph, |
---|
| 500 | const LowerMap &lower, |
---|
| 501 | const CapacityMap &capacity, |
---|
| 502 | const CostMap &cost, |
---|
| 503 | const SupplyMap &supply ) : |
---|
| 504 | _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), |
---|
| 505 | _cost(_graph), _supply(_graph), _flow(_graph), |
---|
| 506 | _potential(_graph), _depth(_graph), _parent(_graph), |
---|
| 507 | _pred_edge(_graph), _thread(_graph), _forward(_graph), |
---|
| 508 | _state(_graph), _red_cost(_graph, _cost, _potential), |
---|
| 509 | _flow_result(graph), _potential_result(graph), |
---|
| 510 | _node_ref(graph), _edge_ref(graph) |
---|
[2440] | 511 | { |
---|
| 512 | // Checking the sum of supply values |
---|
| 513 | Supply sum = 0; |
---|
[2575] | 514 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
---|
| 515 | sum += supply[n]; |
---|
| 516 | if (!(_valid_supply = sum == 0)) return; |
---|
[2440] | 517 | |
---|
[2575] | 518 | // Copying _graph_ref to _graph |
---|
| 519 | _graph.reserveNode(countNodes(_graph_ref) + 1); |
---|
| 520 | _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref)); |
---|
| 521 | copyGraph(_graph, _graph_ref) |
---|
| 522 | .edgeMap(_cost, cost) |
---|
| 523 | .nodeRef(_node_ref) |
---|
| 524 | .edgeRef(_edge_ref) |
---|
[2556] | 525 | .run(); |
---|
[2440] | 526 | |
---|
[2556] | 527 | // Removing non-zero lower bounds |
---|
[2575] | 528 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { |
---|
| 529 | _capacity[_edge_ref[e]] = capacity[e] - lower[e]; |
---|
[2440] | 530 | } |
---|
[2575] | 531 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { |
---|
| 532 | Supply s = supply[n]; |
---|
| 533 | for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
---|
| 534 | s += lower[e]; |
---|
| 535 | for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
---|
| 536 | s -= lower[e]; |
---|
| 537 | _supply[_node_ref[n]] = s; |
---|
[2440] | 538 | } |
---|
| 539 | } |
---|
| 540 | |
---|
| 541 | /// \brief General constructor of the class (without lower bounds). |
---|
| 542 | /// |
---|
| 543 | /// General constructor of the class (without lower bounds). |
---|
| 544 | /// |
---|
[2575] | 545 | /// \param graph The directed graph the algorithm runs on. |
---|
| 546 | /// \param capacity The capacities (upper bounds) of the edges. |
---|
| 547 | /// \param cost The cost (length) values of the edges. |
---|
| 548 | /// \param supply The supply values of the nodes (signed). |
---|
| 549 | NetworkSimplex( const Graph &graph, |
---|
| 550 | const CapacityMap &capacity, |
---|
| 551 | const CostMap &cost, |
---|
| 552 | const SupplyMap &supply ) : |
---|
| 553 | _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), |
---|
| 554 | _cost(_graph), _supply(_graph), _flow(_graph), |
---|
| 555 | _potential(_graph), _depth(_graph), _parent(_graph), |
---|
| 556 | _pred_edge(_graph), _thread(_graph), _forward(_graph), |
---|
| 557 | _state(_graph), _red_cost(_graph, _cost, _potential), |
---|
| 558 | _flow_result(graph), _potential_result(graph), |
---|
| 559 | _node_ref(graph), _edge_ref(graph) |
---|
[2440] | 560 | { |
---|
| 561 | // Checking the sum of supply values |
---|
| 562 | Supply sum = 0; |
---|
[2575] | 563 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
---|
| 564 | sum += supply[n]; |
---|
| 565 | if (!(_valid_supply = sum == 0)) return; |
---|
[2440] | 566 | |
---|
[2575] | 567 | // Copying _graph_ref to graph |
---|
| 568 | copyGraph(_graph, _graph_ref) |
---|
| 569 | .edgeMap(_capacity, capacity) |
---|
| 570 | .edgeMap(_cost, cost) |
---|
| 571 | .nodeMap(_supply, supply) |
---|
| 572 | .nodeRef(_node_ref) |
---|
| 573 | .edgeRef(_edge_ref) |
---|
[2556] | 574 | .run(); |
---|
[2440] | 575 | } |
---|
| 576 | |
---|
| 577 | /// \brief Simple constructor of the class (with lower bounds). |
---|
| 578 | /// |
---|
| 579 | /// Simple constructor of the class (with lower bounds). |
---|
| 580 | /// |
---|
[2575] | 581 | /// \param graph The directed graph the algorithm runs on. |
---|
| 582 | /// \param lower The lower bounds of the edges. |
---|
| 583 | /// \param capacity The capacities (upper bounds) of the edges. |
---|
| 584 | /// \param cost The cost (length) values of the edges. |
---|
| 585 | /// \param s The source node. |
---|
| 586 | /// \param t The target node. |
---|
| 587 | /// \param flow_value The required amount of flow from node \c s |
---|
| 588 | /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
---|
| 589 | NetworkSimplex( const Graph &graph, |
---|
| 590 | const LowerMap &lower, |
---|
| 591 | const CapacityMap &capacity, |
---|
| 592 | const CostMap &cost, |
---|
| 593 | typename Graph::Node s, |
---|
| 594 | typename Graph::Node t, |
---|
| 595 | typename SupplyMap::Value flow_value ) : |
---|
| 596 | _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph), |
---|
| 597 | _cost(_graph), _supply(_graph), _flow(_graph), |
---|
| 598 | _potential(_graph), _depth(_graph), _parent(_graph), |
---|
| 599 | _pred_edge(_graph), _thread(_graph), _forward(_graph), |
---|
| 600 | _state(_graph), _red_cost(_graph, _cost, _potential), |
---|
| 601 | _flow_result(graph), _potential_result(graph), |
---|
| 602 | _node_ref(graph), _edge_ref(graph) |
---|
[2440] | 603 | { |
---|
[2575] | 604 | // Copying _graph_ref to graph |
---|
| 605 | copyGraph(_graph, _graph_ref) |
---|
| 606 | .edgeMap(_cost, cost) |
---|
| 607 | .nodeRef(_node_ref) |
---|
| 608 | .edgeRef(_edge_ref) |
---|
[2556] | 609 | .run(); |
---|
[2440] | 610 | |
---|
[2556] | 611 | // Removing non-zero lower bounds |
---|
[2575] | 612 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) { |
---|
| 613 | _capacity[_edge_ref[e]] = capacity[e] - lower[e]; |
---|
[2440] | 614 | } |
---|
[2575] | 615 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) { |
---|
| 616 | Supply sum = 0; |
---|
| 617 | if (n == s) sum = flow_value; |
---|
| 618 | if (n == t) sum = -flow_value; |
---|
| 619 | for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
---|
| 620 | sum += lower[e]; |
---|
| 621 | for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e) |
---|
| 622 | sum -= lower[e]; |
---|
| 623 | _supply[_node_ref[n]] = sum; |
---|
[2440] | 624 | } |
---|
[2575] | 625 | _valid_supply = true; |
---|
[2440] | 626 | } |
---|
| 627 | |
---|
| 628 | /// \brief Simple constructor of the class (without lower bounds). |
---|
| 629 | /// |
---|
| 630 | /// Simple constructor of the class (without lower bounds). |
---|
| 631 | /// |
---|
[2575] | 632 | /// \param graph The directed graph the algorithm runs on. |
---|
| 633 | /// \param capacity The capacities (upper bounds) of the edges. |
---|
| 634 | /// \param cost The cost (length) values of the edges. |
---|
| 635 | /// \param s The source node. |
---|
| 636 | /// \param t The target node. |
---|
| 637 | /// \param flow_value The required amount of flow from node \c s |
---|
| 638 | /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
---|
| 639 | NetworkSimplex( const Graph &graph, |
---|
| 640 | const CapacityMap &capacity, |
---|
| 641 | const CostMap &cost, |
---|
| 642 | typename Graph::Node s, |
---|
| 643 | typename Graph::Node t, |
---|
| 644 | typename SupplyMap::Value flow_value ) : |
---|
| 645 | _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph), |
---|
| 646 | _cost(_graph), _supply(_graph, 0), _flow(_graph), |
---|
| 647 | _potential(_graph), _depth(_graph), _parent(_graph), |
---|
| 648 | _pred_edge(_graph), _thread(_graph), _forward(_graph), |
---|
| 649 | _state(_graph), _red_cost(_graph, _cost, _potential), |
---|
| 650 | _flow_result(graph), _potential_result(graph), |
---|
| 651 | _node_ref(graph), _edge_ref(graph) |
---|
[2440] | 652 | { |
---|
[2575] | 653 | // Copying _graph_ref to graph |
---|
| 654 | copyGraph(_graph, _graph_ref) |
---|
| 655 | .edgeMap(_capacity, capacity) |
---|
| 656 | .edgeMap(_cost, cost) |
---|
| 657 | .nodeRef(_node_ref) |
---|
| 658 | .edgeRef(_edge_ref) |
---|
[2556] | 659 | .run(); |
---|
[2575] | 660 | _supply[_node_ref[s]] = flow_value; |
---|
| 661 | _supply[_node_ref[t]] = -flow_value; |
---|
| 662 | _valid_supply = true; |
---|
[2440] | 663 | } |
---|
| 664 | |
---|
[2556] | 665 | /// \brief Runs the algorithm. |
---|
| 666 | /// |
---|
| 667 | /// Runs the algorithm. |
---|
| 668 | /// |
---|
[2575] | 669 | /// \param pivot_rule The pivot rule that is used during the |
---|
| 670 | /// algorithm. |
---|
| 671 | /// |
---|
| 672 | /// The available pivot rules: |
---|
| 673 | /// |
---|
| 674 | /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in |
---|
| 675 | /// a wraparound fashion in every iteration |
---|
| 676 | /// (\ref FirstEligiblePivotRule). |
---|
| 677 | /// |
---|
| 678 | /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in |
---|
| 679 | /// every iteration (\ref BestEligiblePivotRule). |
---|
| 680 | /// |
---|
| 681 | /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in |
---|
| 682 | /// every iteration in a wraparound fashion and the best eligible |
---|
| 683 | /// edge is selected from this block (\ref BlockSearchPivotRule). |
---|
| 684 | /// |
---|
| 685 | /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are |
---|
| 686 | /// examined in every iteration in a wraparound fashion and the best |
---|
| 687 | /// one is selected from them (\ref LimitedSearchPivotRule). |
---|
| 688 | /// |
---|
| 689 | /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is |
---|
| 690 | /// built from eligible edges and it is used for edge selection in |
---|
| 691 | /// the following minor iterations (\ref CandidateListPivotRule). |
---|
| 692 | /// |
---|
| 693 | /// - COMBINED_PIVOT This is a combined version of the two fastest |
---|
| 694 | /// pivot rules. |
---|
| 695 | /// For rather sparse graphs \ref LimitedSearchPivotRule |
---|
| 696 | /// "Limited Search" implementation is used, otherwise |
---|
| 697 | /// \ref BlockSearchPivotRule "Block Search" pivot rule is used. |
---|
| 698 | /// According to our benchmark tests this combined method is the |
---|
| 699 | /// most efficient. |
---|
| 700 | /// |
---|
[2556] | 701 | /// \return \c true if a feasible flow can be found. |
---|
[2575] | 702 | bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) { |
---|
| 703 | return init() && start(pivot_rule); |
---|
[2556] | 704 | } |
---|
| 705 | |
---|
[2575] | 706 | /// \brief Returns a const reference to the edge map storing the |
---|
| 707 | /// found flow. |
---|
[2440] | 708 | /// |
---|
[2575] | 709 | /// Returns a const reference to the edge map storing the found flow. |
---|
[2440] | 710 | /// |
---|
| 711 | /// \pre \ref run() must be called before using this function. |
---|
| 712 | const FlowMap& flowMap() const { |
---|
[2575] | 713 | return _flow_result; |
---|
[2440] | 714 | } |
---|
| 715 | |
---|
[2575] | 716 | /// \brief Returns a const reference to the node map storing the |
---|
| 717 | /// found potentials (the dual solution). |
---|
[2440] | 718 | /// |
---|
[2575] | 719 | /// Returns a const reference to the node map storing the found |
---|
| 720 | /// potentials (the dual solution). |
---|
[2440] | 721 | /// |
---|
| 722 | /// \pre \ref run() must be called before using this function. |
---|
| 723 | const PotentialMap& potentialMap() const { |
---|
[2575] | 724 | return _potential_result; |
---|
[2440] | 725 | } |
---|
| 726 | |
---|
| 727 | /// \brief Returns the total cost of the found flow. |
---|
| 728 | /// |
---|
| 729 | /// Returns the total cost of the found flow. The complexity of the |
---|
| 730 | /// function is \f$ O(e) \f$. |
---|
| 731 | /// |
---|
| 732 | /// \pre \ref run() must be called before using this function. |
---|
| 733 | Cost totalCost() const { |
---|
| 734 | Cost c = 0; |
---|
[2575] | 735 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
---|
| 736 | c += _flow_result[e] * _cost[_edge_ref[e]]; |
---|
[2440] | 737 | return c; |
---|
| 738 | } |
---|
| 739 | |
---|
[2575] | 740 | private: |
---|
[2440] | 741 | |
---|
| 742 | /// \brief Extends the underlaying graph and initializes all the |
---|
| 743 | /// node and edge maps. |
---|
| 744 | bool init() { |
---|
[2575] | 745 | if (!_valid_supply) return false; |
---|
[2440] | 746 | |
---|
| 747 | // Initializing state and flow maps |
---|
[2575] | 748 | for (EdgeIt e(_graph); e != INVALID; ++e) { |
---|
| 749 | _flow[e] = 0; |
---|
| 750 | _state[e] = STATE_LOWER; |
---|
[2440] | 751 | } |
---|
| 752 | |
---|
| 753 | // Adding an artificial root node to the graph |
---|
[2575] | 754 | _root = _graph.addNode(); |
---|
| 755 | _parent[_root] = INVALID; |
---|
| 756 | _pred_edge[_root] = INVALID; |
---|
| 757 | _depth[_root] = 0; |
---|
| 758 | _supply[_root] = 0; |
---|
| 759 | _potential[_root] = 0; |
---|
[2440] | 760 | |
---|
| 761 | // Adding artificial edges to the graph and initializing the node |
---|
| 762 | // maps of the spanning tree data structure |
---|
[2575] | 763 | Node last = _root; |
---|
[2440] | 764 | Edge e; |
---|
| 765 | Cost max_cost = std::numeric_limits<Cost>::max() / 4; |
---|
[2575] | 766 | for (NodeIt u(_graph); u != INVALID; ++u) { |
---|
| 767 | if (u == _root) continue; |
---|
| 768 | _thread[last] = u; |
---|
[2556] | 769 | last = u; |
---|
[2575] | 770 | _parent[u] = _root; |
---|
| 771 | _depth[u] = 1; |
---|
| 772 | if (_supply[u] >= 0) { |
---|
| 773 | e = _graph.addEdge(u, _root); |
---|
| 774 | _flow[e] = _supply[u]; |
---|
| 775 | _forward[u] = true; |
---|
| 776 | _potential[u] = -max_cost; |
---|
[2556] | 777 | } else { |
---|
[2575] | 778 | e = _graph.addEdge(_root, u); |
---|
| 779 | _flow[e] = -_supply[u]; |
---|
| 780 | _forward[u] = false; |
---|
| 781 | _potential[u] = max_cost; |
---|
[2556] | 782 | } |
---|
[2575] | 783 | _cost[e] = max_cost; |
---|
| 784 | _capacity[e] = std::numeric_limits<Capacity>::max(); |
---|
| 785 | _state[e] = STATE_TREE; |
---|
| 786 | _pred_edge[u] = e; |
---|
[2440] | 787 | } |
---|
[2575] | 788 | _thread[last] = _root; |
---|
[2440] | 789 | |
---|
[2575] | 790 | return true; |
---|
[2440] | 791 | } |
---|
| 792 | |
---|
[2575] | 793 | /// Finds the join node. |
---|
| 794 | Node findJoinNode() { |
---|
| 795 | Node u = _graph.source(_in_edge); |
---|
| 796 | Node v = _graph.target(_in_edge); |
---|
| 797 | while (u != v) { |
---|
| 798 | if (_depth[u] == _depth[v]) { |
---|
| 799 | u = _parent[u]; |
---|
| 800 | v = _parent[v]; |
---|
[2556] | 801 | } |
---|
[2575] | 802 | else if (_depth[u] > _depth[v]) u = _parent[u]; |
---|
| 803 | else v = _parent[v]; |
---|
[2440] | 804 | } |
---|
| 805 | return u; |
---|
| 806 | } |
---|
| 807 | |
---|
| 808 | /// \brief Finds the leaving edge of the cycle. Returns \c true if |
---|
| 809 | /// the leaving edge is not the same as the entering edge. |
---|
| 810 | bool findLeavingEdge() { |
---|
| 811 | // Initializing first and second nodes according to the direction |
---|
| 812 | // of the cycle |
---|
[2575] | 813 | if (_state[_in_edge] == STATE_LOWER) { |
---|
| 814 | first = _graph.source(_in_edge); |
---|
| 815 | second = _graph.target(_in_edge); |
---|
[2440] | 816 | } else { |
---|
[2575] | 817 | first = _graph.target(_in_edge); |
---|
| 818 | second = _graph.source(_in_edge); |
---|
[2440] | 819 | } |
---|
[2575] | 820 | delta = _capacity[_in_edge]; |
---|
[2440] | 821 | bool result = false; |
---|
| 822 | Capacity d; |
---|
| 823 | Edge e; |
---|
| 824 | |
---|
| 825 | // Searching the cycle along the path form the first node to the |
---|
| 826 | // root node |
---|
[2575] | 827 | for (Node u = first; u != join; u = _parent[u]) { |
---|
| 828 | e = _pred_edge[u]; |
---|
| 829 | d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e]; |
---|
[2556] | 830 | if (d < delta) { |
---|
| 831 | delta = d; |
---|
| 832 | u_out = u; |
---|
| 833 | u_in = first; |
---|
| 834 | v_in = second; |
---|
| 835 | result = true; |
---|
| 836 | } |
---|
[2440] | 837 | } |
---|
| 838 | // Searching the cycle along the path form the second node to the |
---|
| 839 | // root node |
---|
[2575] | 840 | for (Node u = second; u != join; u = _parent[u]) { |
---|
| 841 | e = _pred_edge[u]; |
---|
| 842 | d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e]; |
---|
[2556] | 843 | if (d <= delta) { |
---|
| 844 | delta = d; |
---|
| 845 | u_out = u; |
---|
| 846 | u_in = second; |
---|
| 847 | v_in = first; |
---|
| 848 | result = true; |
---|
| 849 | } |
---|
[2440] | 850 | } |
---|
| 851 | return result; |
---|
| 852 | } |
---|
| 853 | |
---|
[2575] | 854 | /// Changes \c flow and \c state edge maps. |
---|
[2440] | 855 | void changeFlows(bool change) { |
---|
| 856 | // Augmenting along the cycle |
---|
| 857 | if (delta > 0) { |
---|
[2575] | 858 | Capacity val = _state[_in_edge] * delta; |
---|
| 859 | _flow[_in_edge] += val; |
---|
| 860 | for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) { |
---|
| 861 | _flow[_pred_edge[u]] += _forward[u] ? -val : val; |
---|
[2556] | 862 | } |
---|
[2575] | 863 | for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) { |
---|
| 864 | _flow[_pred_edge[u]] += _forward[u] ? val : -val; |
---|
[2556] | 865 | } |
---|
[2440] | 866 | } |
---|
| 867 | // Updating the state of the entering and leaving edges |
---|
| 868 | if (change) { |
---|
[2575] | 869 | _state[_in_edge] = STATE_TREE; |
---|
| 870 | _state[_pred_edge[u_out]] = |
---|
| 871 | (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; |
---|
[2440] | 872 | } else { |
---|
[2575] | 873 | _state[_in_edge] = -_state[_in_edge]; |
---|
[2440] | 874 | } |
---|
| 875 | } |
---|
| 876 | |
---|
[2575] | 877 | /// Updates \c thread and \c parent node maps. |
---|
[2440] | 878 | void updateThreadParent() { |
---|
| 879 | Node u; |
---|
[2575] | 880 | v_out = _parent[u_out]; |
---|
[2440] | 881 | |
---|
| 882 | // Handling the case when join and v_out coincide |
---|
| 883 | bool par_first = false; |
---|
| 884 | if (join == v_out) { |
---|
[2575] | 885 | for (u = join; u != u_in && u != v_in; u = _thread[u]) ; |
---|
[2556] | 886 | if (u == v_in) { |
---|
| 887 | par_first = true; |
---|
[2575] | 888 | while (_thread[u] != u_out) u = _thread[u]; |
---|
[2556] | 889 | first = u; |
---|
| 890 | } |
---|
[2440] | 891 | } |
---|
| 892 | |
---|
| 893 | // Finding the last successor of u_in (u) and the node after it |
---|
| 894 | // (right) according to the thread index |
---|
[2575] | 895 | for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ; |
---|
| 896 | right = _thread[u]; |
---|
| 897 | if (_thread[v_in] == u_out) { |
---|
| 898 | for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ; |
---|
| 899 | if (last == u_out) last = _thread[last]; |
---|
[2440] | 900 | } |
---|
[2575] | 901 | else last = _thread[v_in]; |
---|
[2440] | 902 | |
---|
| 903 | // Updating stem nodes |
---|
[2575] | 904 | _thread[v_in] = stem = u_in; |
---|
[2440] | 905 | par_stem = v_in; |
---|
| 906 | while (stem != u_out) { |
---|
[2575] | 907 | _thread[u] = new_stem = _parent[stem]; |
---|
[2440] | 908 | |
---|
[2556] | 909 | // Finding the node just before the stem node (u) according to |
---|
| 910 | // the original thread index |
---|
[2575] | 911 | for (u = new_stem; _thread[u] != stem; u = _thread[u]) ; |
---|
| 912 | _thread[u] = right; |
---|
[2440] | 913 | |
---|
[2556] | 914 | // Changing the parent node of stem and shifting stem and |
---|
| 915 | // par_stem nodes |
---|
[2575] | 916 | _parent[stem] = par_stem; |
---|
[2556] | 917 | par_stem = stem; |
---|
| 918 | stem = new_stem; |
---|
[2440] | 919 | |
---|
[2556] | 920 | // Finding the last successor of stem (u) and the node after it |
---|
| 921 | // (right) according to the thread index |
---|
[2575] | 922 | for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ; |
---|
| 923 | right = _thread[u]; |
---|
[2440] | 924 | } |
---|
[2575] | 925 | _parent[u_out] = par_stem; |
---|
| 926 | _thread[u] = last; |
---|
[2440] | 927 | |
---|
| 928 | if (join == v_out && par_first) { |
---|
[2575] | 929 | if (first != v_in) _thread[first] = right; |
---|
[2440] | 930 | } else { |
---|
[2575] | 931 | for (u = v_out; _thread[u] != u_out; u = _thread[u]) ; |
---|
| 932 | _thread[u] = right; |
---|
[2440] | 933 | } |
---|
| 934 | } |
---|
| 935 | |
---|
[2575] | 936 | /// Updates \c pred_edge and \c forward node maps. |
---|
[2440] | 937 | void updatePredEdge() { |
---|
| 938 | Node u = u_out, v; |
---|
| 939 | while (u != u_in) { |
---|
[2575] | 940 | v = _parent[u]; |
---|
| 941 | _pred_edge[u] = _pred_edge[v]; |
---|
| 942 | _forward[u] = !_forward[v]; |
---|
[2556] | 943 | u = v; |
---|
[2440] | 944 | } |
---|
[2575] | 945 | _pred_edge[u_in] = _in_edge; |
---|
| 946 | _forward[u_in] = (u_in == _graph.source(_in_edge)); |
---|
[2440] | 947 | } |
---|
| 948 | |
---|
[2575] | 949 | /// Updates \c depth and \c potential node maps. |
---|
[2440] | 950 | void updateDepthPotential() { |
---|
[2575] | 951 | _depth[u_in] = _depth[v_in] + 1; |
---|
| 952 | _potential[u_in] = _forward[u_in] ? |
---|
| 953 | _potential[v_in] - _cost[_pred_edge[u_in]] : |
---|
| 954 | _potential[v_in] + _cost[_pred_edge[u_in]]; |
---|
[2440] | 955 | |
---|
[2575] | 956 | Node u = _thread[u_in], v; |
---|
[2440] | 957 | while (true) { |
---|
[2575] | 958 | v = _parent[u]; |
---|
[2556] | 959 | if (v == INVALID) break; |
---|
[2575] | 960 | _depth[u] = _depth[v] + 1; |
---|
| 961 | _potential[u] = _forward[u] ? |
---|
| 962 | _potential[v] - _cost[_pred_edge[u]] : |
---|
| 963 | _potential[v] + _cost[_pred_edge[u]]; |
---|
| 964 | if (_depth[u] <= _depth[v_in]) break; |
---|
| 965 | u = _thread[u]; |
---|
[2440] | 966 | } |
---|
| 967 | } |
---|
| 968 | |
---|
[2575] | 969 | /// Executes the algorithm. |
---|
| 970 | bool start(PivotRuleEnum pivot_rule) { |
---|
| 971 | switch (pivot_rule) { |
---|
| 972 | case FIRST_ELIGIBLE_PIVOT: |
---|
| 973 | return start<FirstEligiblePivotRule>(); |
---|
| 974 | case BEST_ELIGIBLE_PIVOT: |
---|
| 975 | return start<BestEligiblePivotRule>(); |
---|
| 976 | case BLOCK_SEARCH_PIVOT: |
---|
| 977 | return start<BlockSearchPivotRule>(); |
---|
| 978 | case LIMITED_SEARCH_PIVOT: |
---|
| 979 | return start<LimitedSearchPivotRule>(); |
---|
| 980 | case CANDIDATE_LIST_PIVOT: |
---|
| 981 | return start<CandidateListPivotRule>(); |
---|
| 982 | case COMBINED_PIVOT: |
---|
| 983 | if ( countEdges(_graph) / countNodes(_graph) <= |
---|
| 984 | COMBINED_PIVOT_MAX_DEG ) |
---|
| 985 | return start<LimitedSearchPivotRule>(); |
---|
| 986 | else |
---|
| 987 | return start<BlockSearchPivotRule>(); |
---|
| 988 | } |
---|
| 989 | return false; |
---|
| 990 | } |
---|
| 991 | |
---|
| 992 | template<class PivotRuleImplementation> |
---|
[2440] | 993 | bool start() { |
---|
[2575] | 994 | PivotRuleImplementation pivot(*this); |
---|
| 995 | |
---|
| 996 | // Executing the network simplex algorithm |
---|
| 997 | while (pivot.findEnteringEdge()) { |
---|
[2556] | 998 | join = findJoinNode(); |
---|
| 999 | bool change = findLeavingEdge(); |
---|
| 1000 | changeFlows(change); |
---|
| 1001 | if (change) { |
---|
| 1002 | updateThreadParent(); |
---|
| 1003 | updatePredEdge(); |
---|
| 1004 | updateDepthPotential(); |
---|
| 1005 | } |
---|
[2440] | 1006 | } |
---|
| 1007 | |
---|
[2575] | 1008 | // Checking if the flow amount equals zero on all the artificial |
---|
| 1009 | // edges |
---|
| 1010 | for (InEdgeIt e(_graph, _root); e != INVALID; ++e) |
---|
| 1011 | if (_flow[e] > 0) return false; |
---|
| 1012 | for (OutEdgeIt e(_graph, _root); e != INVALID; ++e) |
---|
| 1013 | if (_flow[e] > 0) return false; |
---|
[2440] | 1014 | |
---|
[2575] | 1015 | // Copying flow values to _flow_result |
---|
| 1016 | if (_lower) { |
---|
| 1017 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
---|
| 1018 | _flow_result[e] = (*_lower)[e] + _flow[_edge_ref[e]]; |
---|
[2440] | 1019 | } else { |
---|
[2575] | 1020 | for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) |
---|
| 1021 | _flow_result[e] = _flow[_edge_ref[e]]; |
---|
[2440] | 1022 | } |
---|
[2575] | 1023 | // Copying potential values to _potential_result |
---|
| 1024 | for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) |
---|
| 1025 | _potential_result[n] = _potential[_node_ref[n]]; |
---|
[2440] | 1026 | |
---|
| 1027 | return true; |
---|
| 1028 | } |
---|
| 1029 | |
---|
| 1030 | }; //class NetworkSimplex |
---|
| 1031 | |
---|
| 1032 | ///@} |
---|
| 1033 | |
---|
| 1034 | } //namespace lemon |
---|
| 1035 | |
---|
| 1036 | #endif //LEMON_NETWORK_SIMPLEX_H |
---|