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