1 /* -*- mode: C++; indent-tabs-mode: nil; -*- 2 * 3 * This file is a part of LEMON, a generic C++ optimization library. 4 * 5 * Copyright (C) 2003-2013 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_algs 23 /// 24 /// \file 25 /// \brief Network Simplex algorithm for finding a minimum cost flow. 26 27 #include <vector> 28 #include <limits> 29 #include <algorithm> 30 31 #include <lemon/core.h> 32 #include <lemon/math.h> 33 34 namespace lemon { 35 36 /// \addtogroup min_cost_flow_algs 37 /// @{ 38 39 /// \brief Implementation of the primal Network Simplex algorithm 40 /// for finding a \ref min_cost_flow "minimum cost flow". 41 /// 42 /// \ref NetworkSimplex implements the primal Network Simplex algorithm 43 /// for finding a \ref min_cost_flow "minimum cost flow" 44 /// \cite amo93networkflows, \cite dantzig63linearprog, 45 /// \cite kellyoneill91netsimplex. 46 /// This algorithm is a highly efficient specialized version of the 47 /// linear programming simplex method directly for the minimum cost 48 /// flow problem. 49 /// 50 /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest 51 /// implementations available in LEMON for solving this problem. 52 /// (For more information, see \ref min_cost_flow_algs "the module page".) 53 /// Furthermore, this class supports both directions of the supply/demand 54 /// inequality constraints. For more information, see \ref SupplyType. 55 /// 56 /// Most of the parameters of the problem (except for the digraph) 57 /// can be given using separate functions, and the algorithm can be 58 /// executed using the \ref run() function. If some parameters are not 59 /// specified, then default values will be used. 60 /// 61 /// \tparam GR The digraph type the algorithm runs on. 62 /// \tparam V The number type used for flow amounts, capacity bounds 63 /// and supply values in the algorithm. By default, it is \c int. 64 /// \tparam C The number type used for costs and potentials in the 65 /// algorithm. By default, it is the same as \c V. 66 /// 67 /// \warning Both \c V and \c C must be signed number types. 68 /// \warning All input data (capacities, supply values, and costs) must 69 /// be integer. 70 /// 71 /// \note %NetworkSimplex provides five different pivot rule 72 /// implementations, from which the most efficient one is used 73 /// by default. For more information, see \ref PivotRule. 74 template <typename GR, typename V = int, typename C = V> 75 class NetworkSimplex 76 { 77 public: 78 79 /// The type of the flow amounts, capacity bounds and supply values 80 typedef V Value; 81 /// The type of the arc costs 82 typedef C Cost; 83 84 public: 85 86 /// \brief Problem type constants for the \c run() function. 87 /// 88 /// Enum type containing the problem type constants that can be 89 /// returned by the \ref run() function of the algorithm. 90 enum ProblemType { 91 /// The problem has no feasible solution (flow). 92 INFEASIBLE, 93 /// The problem has optimal solution (i.e. it is feasible and 94 /// bounded), and the algorithm has found optimal flow and node 95 /// potentials (primal and dual solutions). 96 OPTIMAL, 97 /// The objective function of the problem is unbounded, i.e. 98 /// there is a directed cycle having negative total cost and 99 /// infinite upper bound. 100 UNBOUNDED 101 }; 102 103 /// \brief Constants for selecting the type of the supply constraints. 104 /// 105 /// Enum type containing constants for selecting the supply type, 106 /// i.e. the direction of the inequalities in the supply/demand 107 /// constraints of the \ref min_cost_flow "minimum cost flow problem". 108 /// 109 /// The default supply type is \c GEQ, the \c LEQ type can be 110 /// selected using \ref supplyType(). 111 /// The equality form is a special case of both supply types. 112 enum SupplyType { 113 /// This option means that there are <em>"greater or equal"</em> 114 /// supply/demand constraints in the definition of the problem. 115 GEQ, 116 /// This option means that there are <em>"less or equal"</em> 117 /// supply/demand constraints in the definition of the problem. 118 LEQ 119 }; 120 121 /// \brief Constants for selecting the pivot rule. 122 /// 123 /// Enum type containing constants for selecting the pivot rule for 124 /// the \ref run() function. 125 /// 126 /// \ref NetworkSimplex provides five different implementations for 127 /// the pivot strategy that significantly affects the running time 128 /// of the algorithm. 129 /// According to experimental tests conducted on various problem 130 /// instances, \ref BLOCK_SEARCH "Block Search" and 131 /// \ref ALTERING_LIST "Altering Candidate List" rules turned out 132 /// to be the most efficient. 133 /// Since \ref BLOCK_SEARCH "Block Search" is a simpler strategy that 134 /// seemed to be slightly more robust, it is used by default. 135 /// However, another pivot rule can easily be selected using the 136 /// \ref run() function with the proper parameter. 137 enum PivotRule { 138 139 /// The \e First \e Eligible pivot rule. 140 /// The next eligible arc is selected in a wraparound fashion 141 /// in every iteration. 142 FIRST_ELIGIBLE, 143 144 /// The \e Best \e Eligible pivot rule. 145 /// The best eligible arc is selected in every iteration. 146 BEST_ELIGIBLE, 147 148 /// The \e Block \e Search pivot rule. 149 /// A specified number of arcs are examined in every iteration 150 /// in a wraparound fashion and the best eligible arc is selected 151 /// from this block. 152 BLOCK_SEARCH, 153 154 /// The \e Candidate \e List pivot rule. 155 /// In a major iteration a candidate list is built from eligible arcs 156 /// in a wraparound fashion and in the following minor iterations 157 /// the best eligible arc is selected from this list. 158 CANDIDATE_LIST, 159 160 /// The \e Altering \e Candidate \e List pivot rule. 161 /// It is a modified version of the Candidate List method. 162 /// It keeps only a few of the best eligible arcs from the former 163 /// candidate list and extends this list in every iteration. 164 ALTERING_LIST 165 }; 166 167 private: 168 169 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 170 171 typedef std::vector<int> IntVector; 172 typedef std::vector<Value> ValueVector; 173 typedef std::vector<Cost> CostVector; 174 typedef std::vector<signed char> CharVector; 175 // Note: vector<signed char> is used instead of vector<ArcState> and 176 // vector<ArcDirection> for efficiency reasons 177 178 // State constants for arcs 179 enum ArcState { 180 STATE_UPPER = -1, 181 STATE_TREE = 0, 182 STATE_LOWER = 1 183 }; 184 185 // Direction constants for tree arcs 186 enum ArcDirection { 187 DIR_DOWN = -1, 188 DIR_UP = 1 189 }; 190 191 private: 192 193 // Data related to the underlying digraph 194 const GR &_graph; 195 int _node_num; 196 int _arc_num; 197 int _all_arc_num; 198 int _search_arc_num; 199 200 // Parameters of the problem 201 bool _has_lower; 202 SupplyType _stype; 203 Value _sum_supply; 204 205 // Data structures for storing the digraph 206 IntNodeMap _node_id; 207 IntArcMap _arc_id; 208 IntVector _source; 209 IntVector _target; 210 bool _arc_mixing; 211 212 // Node and arc data 213 ValueVector _lower; 214 ValueVector _upper; 215 ValueVector _cap; 216 CostVector _cost; 217 ValueVector _supply; 218 ValueVector _flow; 219 CostVector _pi; 220 221 // Data for storing the spanning tree structure 222 IntVector _parent; 223 IntVector _pred; 224 IntVector _thread; 225 IntVector _rev_thread; 226 IntVector _succ_num; 227 IntVector _last_succ; 228 CharVector _pred_dir; 229 CharVector _state; 230 IntVector _dirty_revs; 231 int _root; 232 233 // Temporary data used in the current pivot iteration 234 int in_arc, join, u_in, v_in, u_out, v_out; 235 Value delta; 236 237 const Value MAX_VALUE; 238 239 public: 240 241 /// \brief Constant for infinite upper bounds (capacities). 242 /// 243 /// Constant for infinite upper bounds (capacities). 244 /// It is \c std::numeric_limits<Value>::infinity() if available, 245 /// \c std::numeric_limits<Value>::max() otherwise. 246 const Value INF; 247 248 private: 249 250 // Implementation of the First Eligible pivot rule 251 class FirstEligiblePivotRule 252 { 253 private: 254 255 // References to the NetworkSimplex class 256 const IntVector &_source; 257 const IntVector &_target; 258 const CostVector &_cost; 259 const CharVector &_state; 260 const CostVector &_pi; 261 int &_in_arc; 262 int _search_arc_num; 263 264 // Pivot rule data 265 int _next_arc; 266 267 public: 268 269 // Constructor FirstEligiblePivotRule(NetworkSimplex & ns)270 FirstEligiblePivotRule(NetworkSimplex &ns) : 271 _source(ns._source), _target(ns._target), 272 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 273 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 274 _next_arc(0) 275 {} 276 277 // Find next entering arc findEnteringArc()278 bool findEnteringArc() { 279 Cost c; 280 for (int e = _next_arc; e != _search_arc_num; ++e) { 281 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 282 if (c < 0) { 283 _in_arc = e; 284 _next_arc = e + 1; 285 return true; 286 } 287 } 288 for (int e = 0; e != _next_arc; ++e) { 289 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 290 if (c < 0) { 291 _in_arc = e; 292 _next_arc = e + 1; 293 return true; 294 } 295 } 296 return false; 297 } 298 299 }; //class FirstEligiblePivotRule 300 301 302 // Implementation of the Best Eligible pivot rule 303 class BestEligiblePivotRule 304 { 305 private: 306 307 // References to the NetworkSimplex class 308 const IntVector &_source; 309 const IntVector &_target; 310 const CostVector &_cost; 311 const CharVector &_state; 312 const CostVector &_pi; 313 int &_in_arc; 314 int _search_arc_num; 315 316 public: 317 318 // Constructor BestEligiblePivotRule(NetworkSimplex & ns)319 BestEligiblePivotRule(NetworkSimplex &ns) : 320 _source(ns._source), _target(ns._target), 321 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 322 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) 323 {} 324 325 // Find next entering arc findEnteringArc()326 bool findEnteringArc() { 327 Cost c, min = 0; 328 for (int e = 0; e != _search_arc_num; ++e) { 329 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 330 if (c < min) { 331 min = c; 332 _in_arc = e; 333 } 334 } 335 return min < 0; 336 } 337 338 }; //class BestEligiblePivotRule 339 340 341 // Implementation of the Block Search pivot rule 342 class BlockSearchPivotRule 343 { 344 private: 345 346 // References to the NetworkSimplex class 347 const IntVector &_source; 348 const IntVector &_target; 349 const CostVector &_cost; 350 const CharVector &_state; 351 const CostVector &_pi; 352 int &_in_arc; 353 int _search_arc_num; 354 355 // Pivot rule data 356 int _block_size; 357 int _next_arc; 358 359 public: 360 361 // Constructor BlockSearchPivotRule(NetworkSimplex & ns)362 BlockSearchPivotRule(NetworkSimplex &ns) : 363 _source(ns._source), _target(ns._target), 364 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 365 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 366 _next_arc(0) 367 { 368 // The main parameters of the pivot rule 369 const double BLOCK_SIZE_FACTOR = 1.0; 370 const int MIN_BLOCK_SIZE = 10; 371 372 _block_size = std::max( int(BLOCK_SIZE_FACTOR * 373 std::sqrt(double(_search_arc_num))), 374 MIN_BLOCK_SIZE ); 375 } 376 377 // Find next entering arc findEnteringArc()378 bool findEnteringArc() { 379 Cost c, min = 0; 380 int cnt = _block_size; 381 int e; 382 for (e = _next_arc; e != _search_arc_num; ++e) { 383 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 384 if (c < min) { 385 min = c; 386 _in_arc = e; 387 } 388 if (--cnt == 0) { 389 if (min < 0) goto search_end; 390 cnt = _block_size; 391 } 392 } 393 for (e = 0; e != _next_arc; ++e) { 394 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 395 if (c < min) { 396 min = c; 397 _in_arc = e; 398 } 399 if (--cnt == 0) { 400 if (min < 0) goto search_end; 401 cnt = _block_size; 402 } 403 } 404 if (min >= 0) return false; 405 406 search_end: 407 _next_arc = e; 408 return true; 409 } 410 411 }; //class BlockSearchPivotRule 412 413 414 // Implementation of the Candidate List pivot rule 415 class CandidateListPivotRule 416 { 417 private: 418 419 // References to the NetworkSimplex class 420 const IntVector &_source; 421 const IntVector &_target; 422 const CostVector &_cost; 423 const CharVector &_state; 424 const CostVector &_pi; 425 int &_in_arc; 426 int _search_arc_num; 427 428 // Pivot rule data 429 IntVector _candidates; 430 int _list_length, _minor_limit; 431 int _curr_length, _minor_count; 432 int _next_arc; 433 434 public: 435 436 /// Constructor CandidateListPivotRule(NetworkSimplex & ns)437 CandidateListPivotRule(NetworkSimplex &ns) : 438 _source(ns._source), _target(ns._target), 439 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 440 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 441 _next_arc(0) 442 { 443 // The main parameters of the pivot rule 444 const double LIST_LENGTH_FACTOR = 0.25; 445 const int MIN_LIST_LENGTH = 10; 446 const double MINOR_LIMIT_FACTOR = 0.1; 447 const int MIN_MINOR_LIMIT = 3; 448 449 _list_length = std::max( int(LIST_LENGTH_FACTOR * 450 std::sqrt(double(_search_arc_num))), 451 MIN_LIST_LENGTH ); 452 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), 453 MIN_MINOR_LIMIT ); 454 _curr_length = _minor_count = 0; 455 _candidates.resize(_list_length); 456 } 457 458 /// Find next entering arc findEnteringArc()459 bool findEnteringArc() { 460 Cost min, c; 461 int e; 462 if (_curr_length > 0 && _minor_count < _minor_limit) { 463 // Minor iteration: select the best eligible arc from the 464 // current candidate list 465 ++_minor_count; 466 min = 0; 467 for (int i = 0; i < _curr_length; ++i) { 468 e = _candidates[i]; 469 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 470 if (c < min) { 471 min = c; 472 _in_arc = e; 473 } 474 else if (c >= 0) { 475 _candidates[i--] = _candidates[--_curr_length]; 476 } 477 } 478 if (min < 0) return true; 479 } 480 481 // Major iteration: build a new candidate list 482 min = 0; 483 _curr_length = 0; 484 for (e = _next_arc; e != _search_arc_num; ++e) { 485 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 486 if (c < 0) { 487 _candidates[_curr_length++] = e; 488 if (c < min) { 489 min = c; 490 _in_arc = e; 491 } 492 if (_curr_length == _list_length) goto search_end; 493 } 494 } 495 for (e = 0; e != _next_arc; ++e) { 496 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 497 if (c < 0) { 498 _candidates[_curr_length++] = e; 499 if (c < min) { 500 min = c; 501 _in_arc = e; 502 } 503 if (_curr_length == _list_length) goto search_end; 504 } 505 } 506 if (_curr_length == 0) return false; 507 508 search_end: 509 _minor_count = 1; 510 _next_arc = e; 511 return true; 512 } 513 514 }; //class CandidateListPivotRule 515 516 517 // Implementation of the Altering Candidate List pivot rule 518 class AlteringListPivotRule 519 { 520 private: 521 522 // References to the NetworkSimplex class 523 const IntVector &_source; 524 const IntVector &_target; 525 const CostVector &_cost; 526 const CharVector &_state; 527 const CostVector &_pi; 528 int &_in_arc; 529 int _search_arc_num; 530 531 // Pivot rule data 532 int _block_size, _head_length, _curr_length; 533 int _next_arc; 534 IntVector _candidates; 535 CostVector _cand_cost; 536 537 // Functor class to compare arcs during sort of the candidate list 538 class SortFunc 539 { 540 private: 541 const CostVector &_map; 542 public: SortFunc(const CostVector & map)543 SortFunc(const CostVector &map) : _map(map) {} operator()544 bool operator()(int left, int right) { 545 return _map[left] < _map[right]; 546 } 547 }; 548 549 SortFunc _sort_func; 550 551 public: 552 553 // Constructor AlteringListPivotRule(NetworkSimplex & ns)554 AlteringListPivotRule(NetworkSimplex &ns) : 555 _source(ns._source), _target(ns._target), 556 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 557 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 558 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) 559 { 560 // The main parameters of the pivot rule 561 const double BLOCK_SIZE_FACTOR = 1.0; 562 const int MIN_BLOCK_SIZE = 10; 563 const double HEAD_LENGTH_FACTOR = 0.01; 564 const int MIN_HEAD_LENGTH = 3; 565 566 _block_size = std::max( int(BLOCK_SIZE_FACTOR * 567 std::sqrt(double(_search_arc_num))), 568 MIN_BLOCK_SIZE ); 569 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), 570 MIN_HEAD_LENGTH ); 571 _candidates.resize(_head_length + _block_size); 572 _curr_length = 0; 573 } 574 575 // Find next entering arc findEnteringArc()576 bool findEnteringArc() { 577 // Check the current candidate list 578 int e; 579 Cost c; 580 for (int i = 0; i != _curr_length; ++i) { 581 e = _candidates[i]; 582 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 583 if (c < 0) { 584 _cand_cost[e] = c; 585 } else { 586 _candidates[i--] = _candidates[--_curr_length]; 587 } 588 } 589 590 // Extend the list 591 int cnt = _block_size; 592 int limit = _head_length; 593 594 for (e = _next_arc; e != _search_arc_num; ++e) { 595 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 596 if (c < 0) { 597 _cand_cost[e] = c; 598 _candidates[_curr_length++] = e; 599 } 600 if (--cnt == 0) { 601 if (_curr_length > limit) goto search_end; 602 limit = 0; 603 cnt = _block_size; 604 } 605 } 606 for (e = 0; e != _next_arc; ++e) { 607 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 608 if (c < 0) { 609 _cand_cost[e] = c; 610 _candidates[_curr_length++] = e; 611 } 612 if (--cnt == 0) { 613 if (_curr_length > limit) goto search_end; 614 limit = 0; 615 cnt = _block_size; 616 } 617 } 618 if (_curr_length == 0) return false; 619 620 search_end: 621 622 // Perform partial sort operation on the candidate list 623 int new_length = std::min(_head_length + 1, _curr_length); 624 std::partial_sort(_candidates.begin(), _candidates.begin() + new_length, 625 _candidates.begin() + _curr_length, _sort_func); 626 627 // Select the entering arc and remove it from the list 628 _in_arc = _candidates[0]; 629 _next_arc = e; 630 _candidates[0] = _candidates[new_length - 1]; 631 _curr_length = new_length - 1; 632 return true; 633 } 634 635 }; //class AlteringListPivotRule 636 637 public: 638 639 /// \brief Constructor. 640 /// 641 /// The constructor of the class. 642 /// 643 /// \param graph The digraph the algorithm runs on. 644 /// \param arc_mixing Indicate if the arcs will be stored in a 645 /// mixed order in the internal data structure. 646 /// In general, it leads to similar performance as using the original 647 /// arc order, but it makes the algorithm more robust and in special 648 /// cases, even significantly faster. Therefore, it is enabled by default. 649 NetworkSimplex(const GR& graph, bool arc_mixing = true) : _graph(graph)650 _graph(graph), _node_id(graph), _arc_id(graph), 651 _arc_mixing(arc_mixing), 652 MAX_VALUE(std::numeric_limits<Value>::max()), 653 INF(std::numeric_limits<Value>::has_infinity ? 654 std::numeric_limits<Value>::infinity() : MAX_VALUE) 655 { 656 // Check the number types 657 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, 658 "The flow type of NetworkSimplex must be signed"); 659 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 660 "The cost type of NetworkSimplex must be signed"); 661 662 // Reset data structures 663 reset(); 664 } 665 666 /// \name Parameters 667 /// The parameters of the algorithm can be specified using these 668 /// functions. 669 670 /// @{ 671 672 /// \brief Set the lower bounds on the arcs. 673 /// 674 /// This function sets the lower bounds on the arcs. 675 /// If it is not used before calling \ref run(), the lower bounds 676 /// will be set to zero on all arcs. 677 /// 678 /// \param map An arc map storing the lower bounds. 679 /// Its \c Value type must be convertible to the \c Value type 680 /// of the algorithm. 681 /// 682 /// \return <tt>(*this)</tt> 683 template <typename LowerMap> lowerMap(const LowerMap & map)684 NetworkSimplex& lowerMap(const LowerMap& map) { 685 _has_lower = true; 686 for (ArcIt a(_graph); a != INVALID; ++a) { 687 _lower[_arc_id[a]] = map[a]; 688 } 689 return *this; 690 } 691 692 /// \brief Set the upper bounds (capacities) on the arcs. 693 /// 694 /// This function sets the upper bounds (capacities) on the arcs. 695 /// If it is not used before calling \ref run(), the upper bounds 696 /// will be set to \ref INF on all arcs (i.e. the flow value will be 697 /// unbounded from above). 698 /// 699 /// \param map An arc map storing the upper bounds. 700 /// Its \c Value type must be convertible to the \c Value type 701 /// of the algorithm. 702 /// 703 /// \return <tt>(*this)</tt> 704 template<typename UpperMap> upperMap(const UpperMap & map)705 NetworkSimplex& upperMap(const UpperMap& map) { 706 for (ArcIt a(_graph); a != INVALID; ++a) { 707 _upper[_arc_id[a]] = map[a]; 708 } 709 return *this; 710 } 711 712 /// \brief Set the costs of the arcs. 713 /// 714 /// This function sets the costs of the arcs. 715 /// If it is not used before calling \ref run(), the costs 716 /// will be set to \c 1 on all arcs. 717 /// 718 /// \param map An arc map storing the costs. 719 /// Its \c Value type must be convertible to the \c Cost type 720 /// of the algorithm. 721 /// 722 /// \return <tt>(*this)</tt> 723 template<typename CostMap> costMap(const CostMap & map)724 NetworkSimplex& costMap(const CostMap& map) { 725 for (ArcIt a(_graph); a != INVALID; ++a) { 726 _cost[_arc_id[a]] = map[a]; 727 } 728 return *this; 729 } 730 731 /// \brief Set the supply values of the nodes. 732 /// 733 /// This function sets the supply values of the nodes. 734 /// If neither this function nor \ref stSupply() is used before 735 /// calling \ref run(), the supply of each node will be set to zero. 736 /// 737 /// \param map A node map storing the supply values. 738 /// Its \c Value type must be convertible to the \c Value type 739 /// of the algorithm. 740 /// 741 /// \return <tt>(*this)</tt> 742 /// 743 /// \sa supplyType() 744 template<typename SupplyMap> supplyMap(const SupplyMap & map)745 NetworkSimplex& supplyMap(const SupplyMap& map) { 746 for (NodeIt n(_graph); n != INVALID; ++n) { 747 _supply[_node_id[n]] = map[n]; 748 } 749 return *this; 750 } 751 752 /// \brief Set single source and target nodes and a supply value. 753 /// 754 /// This function sets a single source node and a single target node 755 /// and the required flow value. 756 /// If neither this function nor \ref supplyMap() is used before 757 /// calling \ref run(), the supply of each node will be set to zero. 758 /// 759 /// Using this function has the same effect as using \ref supplyMap() 760 /// with a map in which \c k is assigned to \c s, \c -k is 761 /// assigned to \c t and all other nodes have zero supply value. 762 /// 763 /// \param s The source node. 764 /// \param t The target node. 765 /// \param k The required amount of flow from node \c s to node \c t 766 /// (i.e. the supply of \c s and the demand of \c t). 767 /// 768 /// \return <tt>(*this)</tt> stSupply(const Node & s,const Node & t,Value k)769 NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { 770 for (int i = 0; i != _node_num; ++i) { 771 _supply[i] = 0; 772 } 773 _supply[_node_id[s]] = k; 774 _supply[_node_id[t]] = -k; 775 return *this; 776 } 777 778 /// \brief Set the type of the supply constraints. 779 /// 780 /// This function sets the type of the supply/demand constraints. 781 /// If it is not used before calling \ref run(), the \ref GEQ supply 782 /// type will be used. 783 /// 784 /// For more information, see \ref SupplyType. 785 /// 786 /// \return <tt>(*this)</tt> supplyType(SupplyType supply_type)787 NetworkSimplex& supplyType(SupplyType supply_type) { 788 _stype = supply_type; 789 return *this; 790 } 791 792 /// @} 793 794 /// \name Execution Control 795 /// The algorithm can be executed using \ref run(). 796 797 /// @{ 798 799 /// \brief Run the algorithm. 800 /// 801 /// This function runs the algorithm. 802 /// The paramters can be specified using functions \ref lowerMap(), 803 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 804 /// \ref supplyType(). 805 /// For example, 806 /// \code 807 /// NetworkSimplex<ListDigraph> ns(graph); 808 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 809 /// .supplyMap(sup).run(); 810 /// \endcode 811 /// 812 /// This function can be called more than once. All the given parameters 813 /// are kept for the next call, unless \ref resetParams() or \ref reset() 814 /// is used, thus only the modified parameters have to be set again. 815 /// If the underlying digraph was also modified after the construction 816 /// of the class (or the last \ref reset() call), then the \ref reset() 817 /// function must be called. 818 /// 819 /// \param pivot_rule The pivot rule that will be used during the 820 /// algorithm. For more information, see \ref PivotRule. 821 /// 822 /// \return \c INFEASIBLE if no feasible flow exists, 823 /// \n \c OPTIMAL if the problem has optimal solution 824 /// (i.e. it is feasible and bounded), and the algorithm has found 825 /// optimal flow and node potentials (primal and dual solutions), 826 /// \n \c UNBOUNDED if the objective function of the problem is 827 /// unbounded, i.e. there is a directed cycle having negative total 828 /// cost and infinite upper bound. 829 /// 830 /// \see ProblemType, PivotRule 831 /// \see resetParams(), reset() 832 ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { 833 if (!init()) return INFEASIBLE; 834 return start(pivot_rule); 835 } 836 837 /// \brief Reset all the parameters that have been given before. 838 /// 839 /// This function resets all the paramaters that have been given 840 /// before using functions \ref lowerMap(), \ref upperMap(), 841 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). 842 /// 843 /// It is useful for multiple \ref run() calls. Basically, all the given 844 /// parameters are kept for the next \ref run() call, unless 845 /// \ref resetParams() or \ref reset() is used. 846 /// If the underlying digraph was also modified after the construction 847 /// of the class or the last \ref reset() call, then the \ref reset() 848 /// function must be used, otherwise \ref resetParams() is sufficient. 849 /// 850 /// For example, 851 /// \code 852 /// NetworkSimplex<ListDigraph> ns(graph); 853 /// 854 /// // First run 855 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 856 /// .supplyMap(sup).run(); 857 /// 858 /// // Run again with modified cost map (resetParams() is not called, 859 /// // so only the cost map have to be set again) 860 /// cost[e] += 100; 861 /// ns.costMap(cost).run(); 862 /// 863 /// // Run again from scratch using resetParams() 864 /// // (the lower bounds will be set to zero on all arcs) 865 /// ns.resetParams(); 866 /// ns.upperMap(capacity).costMap(cost) 867 /// .supplyMap(sup).run(); 868 /// \endcode 869 /// 870 /// \return <tt>(*this)</tt> 871 /// 872 /// \see reset(), run() resetParams()873 NetworkSimplex& resetParams() { 874 for (int i = 0; i != _node_num; ++i) { 875 _supply[i] = 0; 876 } 877 for (int i = 0; i != _arc_num; ++i) { 878 _lower[i] = 0; 879 _upper[i] = INF; 880 _cost[i] = 1; 881 } 882 _has_lower = false; 883 _stype = GEQ; 884 return *this; 885 } 886 887 /// \brief Reset the internal data structures and all the parameters 888 /// that have been given before. 889 /// 890 /// This function resets the internal data structures and all the 891 /// paramaters that have been given before using functions \ref lowerMap(), 892 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 893 /// \ref supplyType(). 894 /// 895 /// It is useful for multiple \ref run() calls. Basically, all the given 896 /// parameters are kept for the next \ref run() call, unless 897 /// \ref resetParams() or \ref reset() is used. 898 /// If the underlying digraph was also modified after the construction 899 /// of the class or the last \ref reset() call, then the \ref reset() 900 /// function must be used, otherwise \ref resetParams() is sufficient. 901 /// 902 /// See \ref resetParams() for examples. 903 /// 904 /// \return <tt>(*this)</tt> 905 /// 906 /// \see resetParams(), run() reset()907 NetworkSimplex& reset() { 908 // Resize vectors 909 _node_num = countNodes(_graph); 910 _arc_num = countArcs(_graph); 911 int all_node_num = _node_num + 1; 912 int max_arc_num = _arc_num + 2 * _node_num; 913 914 _source.resize(max_arc_num); 915 _target.resize(max_arc_num); 916 917 _lower.resize(_arc_num); 918 _upper.resize(_arc_num); 919 _cap.resize(max_arc_num); 920 _cost.resize(max_arc_num); 921 _supply.resize(all_node_num); 922 _flow.resize(max_arc_num); 923 _pi.resize(all_node_num); 924 925 _parent.resize(all_node_num); 926 _pred.resize(all_node_num); 927 _pred_dir.resize(all_node_num); 928 _thread.resize(all_node_num); 929 _rev_thread.resize(all_node_num); 930 _succ_num.resize(all_node_num); 931 _last_succ.resize(all_node_num); 932 _state.resize(max_arc_num); 933 934 // Copy the graph 935 int i = 0; 936 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 937 _node_id[n] = i; 938 } 939 if (_arc_mixing && _node_num > 1) { 940 // Store the arcs in a mixed order 941 const int skip = std::max(_arc_num / _node_num, 3); 942 int i = 0, j = 0; 943 for (ArcIt a(_graph); a != INVALID; ++a) { 944 _arc_id[a] = i; 945 _source[i] = _node_id[_graph.source(a)]; 946 _target[i] = _node_id[_graph.target(a)]; 947 if ((i += skip) >= _arc_num) i = ++j; 948 } 949 } else { 950 // Store the arcs in the original order 951 int i = 0; 952 for (ArcIt a(_graph); a != INVALID; ++a, ++i) { 953 _arc_id[a] = i; 954 _source[i] = _node_id[_graph.source(a)]; 955 _target[i] = _node_id[_graph.target(a)]; 956 } 957 } 958 959 // Reset parameters 960 resetParams(); 961 return *this; 962 } 963 964 /// @} 965 966 /// \name Query Functions 967 /// The results of the algorithm can be obtained using these 968 /// functions.\n 969 /// The \ref run() function must be called before using them. 970 971 /// @{ 972 973 /// \brief Return the total cost of the found flow. 974 /// 975 /// This function returns the total cost of the found flow. 976 /// Its complexity is O(m). 977 /// 978 /// \note The return type of the function can be specified as a 979 /// template parameter. For example, 980 /// \code 981 /// ns.totalCost<double>(); 982 /// \endcode 983 /// It is useful if the total cost cannot be stored in the \c Cost 984 /// type of the algorithm, which is the default return type of the 985 /// function. 986 /// 987 /// \pre \ref run() must be called before using this function. 988 template <typename Number> totalCost()989 Number totalCost() const { 990 Number c = 0; 991 for (ArcIt a(_graph); a != INVALID; ++a) { 992 int i = _arc_id[a]; 993 c += Number(_flow[i]) * Number(_cost[i]); 994 } 995 return c; 996 } 997 998 #ifndef DOXYGEN totalCost()999 Cost totalCost() const { 1000 return totalCost<Cost>(); 1001 } 1002 #endif 1003 1004 /// \brief Return the flow on the given arc. 1005 /// 1006 /// This function returns the flow on the given arc. 1007 /// 1008 /// \pre \ref run() must be called before using this function. flow(const Arc & a)1009 Value flow(const Arc& a) const { 1010 return _flow[_arc_id[a]]; 1011 } 1012 1013 /// \brief Copy the flow values (the primal solution) into the 1014 /// given map. 1015 /// 1016 /// This function copies the flow value on each arc into the given 1017 /// map. The \c Value type of the algorithm must be convertible to 1018 /// the \c Value type of the map. 1019 /// 1020 /// \pre \ref run() must be called before using this function. 1021 template <typename FlowMap> flowMap(FlowMap & map)1022 void flowMap(FlowMap &map) const { 1023 for (ArcIt a(_graph); a != INVALID; ++a) { 1024 map.set(a, _flow[_arc_id[a]]); 1025 } 1026 } 1027 1028 /// \brief Return the potential (dual value) of the given node. 1029 /// 1030 /// This function returns the potential (dual value) of the 1031 /// given node. 1032 /// 1033 /// \pre \ref run() must be called before using this function. potential(const Node & n)1034 Cost potential(const Node& n) const { 1035 return _pi[_node_id[n]]; 1036 } 1037 1038 /// \brief Copy the potential values (the dual solution) into the 1039 /// given map. 1040 /// 1041 /// This function copies the potential (dual value) of each node 1042 /// into the given map. 1043 /// The \c Cost type of the algorithm must be convertible to the 1044 /// \c Value type of the map. 1045 /// 1046 /// \pre \ref run() must be called before using this function. 1047 template <typename PotentialMap> potentialMap(PotentialMap & map)1048 void potentialMap(PotentialMap &map) const { 1049 for (NodeIt n(_graph); n != INVALID; ++n) { 1050 map.set(n, _pi[_node_id[n]]); 1051 } 1052 } 1053 1054 /// @} 1055 1056 private: 1057 1058 // Initialize internal data structures init()1059 bool init() { 1060 if (_node_num == 0) return false; 1061 1062 // Check the sum of supply values 1063 _sum_supply = 0; 1064 for (int i = 0; i != _node_num; ++i) { 1065 _sum_supply += _supply[i]; 1066 } 1067 if ( !((_stype == GEQ && _sum_supply <= 0) || 1068 (_stype == LEQ && _sum_supply >= 0)) ) return false; 1069 1070 // Check lower and upper bounds 1071 LEMON_DEBUG(checkBoundMaps(), 1072 "Upper bounds must be greater or equal to the lower bounds"); 1073 1074 // Remove non-zero lower bounds 1075 if (_has_lower) { 1076 for (int i = 0; i != _arc_num; ++i) { 1077 Value c = _lower[i]; 1078 if (c >= 0) { 1079 _cap[i] = _upper[i] < MAX_VALUE ? _upper[i] - c : INF; 1080 } else { 1081 _cap[i] = _upper[i] < MAX_VALUE + c ? _upper[i] - c : INF; 1082 } 1083 _supply[_source[i]] -= c; 1084 _supply[_target[i]] += c; 1085 } 1086 } else { 1087 for (int i = 0; i != _arc_num; ++i) { 1088 _cap[i] = _upper[i]; 1089 } 1090 } 1091 1092 // Initialize artifical cost 1093 Cost ART_COST; 1094 if (std::numeric_limits<Cost>::is_exact) { 1095 ART_COST = std::numeric_limits<Cost>::max() / 2 + 1; 1096 } else { 1097 ART_COST = 0; 1098 for (int i = 0; i != _arc_num; ++i) { 1099 if (_cost[i] > ART_COST) ART_COST = _cost[i]; 1100 } 1101 ART_COST = (ART_COST + 1) * _node_num; 1102 } 1103 1104 // Initialize arc maps 1105 for (int i = 0; i != _arc_num; ++i) { 1106 _flow[i] = 0; 1107 _state[i] = STATE_LOWER; 1108 } 1109 1110 // Set data for the artificial root node 1111 _root = _node_num; 1112 _parent[_root] = -1; 1113 _pred[_root] = -1; 1114 _thread[_root] = 0; 1115 _rev_thread[0] = _root; 1116 _succ_num[_root] = _node_num + 1; 1117 _last_succ[_root] = _root - 1; 1118 _supply[_root] = -_sum_supply; 1119 _pi[_root] = 0; 1120 1121 // Add artificial arcs and initialize the spanning tree data structure 1122 if (_sum_supply == 0) { 1123 // EQ supply constraints 1124 _search_arc_num = _arc_num; 1125 _all_arc_num = _arc_num + _node_num; 1126 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1127 _parent[u] = _root; 1128 _pred[u] = e; 1129 _thread[u] = u + 1; 1130 _rev_thread[u + 1] = u; 1131 _succ_num[u] = 1; 1132 _last_succ[u] = u; 1133 _cap[e] = INF; 1134 _state[e] = STATE_TREE; 1135 if (_supply[u] >= 0) { 1136 _pred_dir[u] = DIR_UP; 1137 _pi[u] = 0; 1138 _source[e] = u; 1139 _target[e] = _root; 1140 _flow[e] = _supply[u]; 1141 _cost[e] = 0; 1142 } else { 1143 _pred_dir[u] = DIR_DOWN; 1144 _pi[u] = ART_COST; 1145 _source[e] = _root; 1146 _target[e] = u; 1147 _flow[e] = -_supply[u]; 1148 _cost[e] = ART_COST; 1149 } 1150 } 1151 } 1152 else if (_sum_supply > 0) { 1153 // LEQ supply constraints 1154 _search_arc_num = _arc_num + _node_num; 1155 int f = _arc_num + _node_num; 1156 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1157 _parent[u] = _root; 1158 _thread[u] = u + 1; 1159 _rev_thread[u + 1] = u; 1160 _succ_num[u] = 1; 1161 _last_succ[u] = u; 1162 if (_supply[u] >= 0) { 1163 _pred_dir[u] = DIR_UP; 1164 _pi[u] = 0; 1165 _pred[u] = e; 1166 _source[e] = u; 1167 _target[e] = _root; 1168 _cap[e] = INF; 1169 _flow[e] = _supply[u]; 1170 _cost[e] = 0; 1171 _state[e] = STATE_TREE; 1172 } else { 1173 _pred_dir[u] = DIR_DOWN; 1174 _pi[u] = ART_COST; 1175 _pred[u] = f; 1176 _source[f] = _root; 1177 _target[f] = u; 1178 _cap[f] = INF; 1179 _flow[f] = -_supply[u]; 1180 _cost[f] = ART_COST; 1181 _state[f] = STATE_TREE; 1182 _source[e] = u; 1183 _target[e] = _root; 1184 _cap[e] = INF; 1185 _flow[e] = 0; 1186 _cost[e] = 0; 1187 _state[e] = STATE_LOWER; 1188 ++f; 1189 } 1190 } 1191 _all_arc_num = f; 1192 } 1193 else { 1194 // GEQ supply constraints 1195 _search_arc_num = _arc_num + _node_num; 1196 int f = _arc_num + _node_num; 1197 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1198 _parent[u] = _root; 1199 _thread[u] = u + 1; 1200 _rev_thread[u + 1] = u; 1201 _succ_num[u] = 1; 1202 _last_succ[u] = u; 1203 if (_supply[u] <= 0) { 1204 _pred_dir[u] = DIR_DOWN; 1205 _pi[u] = 0; 1206 _pred[u] = e; 1207 _source[e] = _root; 1208 _target[e] = u; 1209 _cap[e] = INF; 1210 _flow[e] = -_supply[u]; 1211 _cost[e] = 0; 1212 _state[e] = STATE_TREE; 1213 } else { 1214 _pred_dir[u] = DIR_UP; 1215 _pi[u] = -ART_COST; 1216 _pred[u] = f; 1217 _source[f] = u; 1218 _target[f] = _root; 1219 _cap[f] = INF; 1220 _flow[f] = _supply[u]; 1221 _state[f] = STATE_TREE; 1222 _cost[f] = ART_COST; 1223 _source[e] = _root; 1224 _target[e] = u; 1225 _cap[e] = INF; 1226 _flow[e] = 0; 1227 _cost[e] = 0; 1228 _state[e] = STATE_LOWER; 1229 ++f; 1230 } 1231 } 1232 _all_arc_num = f; 1233 } 1234 1235 return true; 1236 } 1237 1238 // Check if the upper bound is greater than or equal to the lower bound 1239 // on each arc. checkBoundMaps()1240 bool checkBoundMaps() { 1241 for (int j = 0; j != _arc_num; ++j) { 1242 if (_upper[j] < _lower[j]) return false; 1243 } 1244 return true; 1245 } 1246 1247 // Find the join node findJoinNode()1248 void findJoinNode() { 1249 int u = _source[in_arc]; 1250 int v = _target[in_arc]; 1251 while (u != v) { 1252 if (_succ_num[u] < _succ_num[v]) { 1253 u = _parent[u]; 1254 } else { 1255 v = _parent[v]; 1256 } 1257 } 1258 join = u; 1259 } 1260 1261 // Find the leaving arc of the cycle and returns true if the 1262 // leaving arc is not the same as the entering arc findLeavingArc()1263 bool findLeavingArc() { 1264 // Initialize first and second nodes according to the direction 1265 // of the cycle 1266 int first, second; 1267 if (_state[in_arc] == STATE_LOWER) { 1268 first = _source[in_arc]; 1269 second = _target[in_arc]; 1270 } else { 1271 first = _target[in_arc]; 1272 second = _source[in_arc]; 1273 } 1274 delta = _cap[in_arc]; 1275 int result = 0; 1276 Value c, d; 1277 int e; 1278 1279 // Search the cycle form the first node to the join node 1280 for (int u = first; u != join; u = _parent[u]) { 1281 e = _pred[u]; 1282 d = _flow[e]; 1283 if (_pred_dir[u] == DIR_DOWN) { 1284 c = _cap[e]; 1285 d = c >= MAX_VALUE ? INF : c - d; 1286 } 1287 if (d < delta) { 1288 delta = d; 1289 u_out = u; 1290 result = 1; 1291 } 1292 } 1293 1294 // Search the cycle form the second node to the join node 1295 for (int u = second; u != join; u = _parent[u]) { 1296 e = _pred[u]; 1297 d = _flow[e]; 1298 if (_pred_dir[u] == DIR_UP) { 1299 c = _cap[e]; 1300 d = c >= MAX_VALUE ? INF : c - d; 1301 } 1302 if (d <= delta) { 1303 delta = d; 1304 u_out = u; 1305 result = 2; 1306 } 1307 } 1308 1309 if (result == 1) { 1310 u_in = first; 1311 v_in = second; 1312 } else { 1313 u_in = second; 1314 v_in = first; 1315 } 1316 return result != 0; 1317 } 1318 1319 // Change _flow and _state vectors changeFlow(bool change)1320 void changeFlow(bool change) { 1321 // Augment along the cycle 1322 if (delta > 0) { 1323 Value val = _state[in_arc] * delta; 1324 _flow[in_arc] += val; 1325 for (int u = _source[in_arc]; u != join; u = _parent[u]) { 1326 _flow[_pred[u]] -= _pred_dir[u] * val; 1327 } 1328 for (int u = _target[in_arc]; u != join; u = _parent[u]) { 1329 _flow[_pred[u]] += _pred_dir[u] * val; 1330 } 1331 } 1332 // Update the state of the entering and leaving arcs 1333 if (change) { 1334 _state[in_arc] = STATE_TREE; 1335 _state[_pred[u_out]] = 1336 (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; 1337 } else { 1338 _state[in_arc] = -_state[in_arc]; 1339 } 1340 } 1341 1342 // Update the tree structure updateTreeStructure()1343 void updateTreeStructure() { 1344 int old_rev_thread = _rev_thread[u_out]; 1345 int old_succ_num = _succ_num[u_out]; 1346 int old_last_succ = _last_succ[u_out]; 1347 v_out = _parent[u_out]; 1348 1349 // Check if u_in and u_out coincide 1350 if (u_in == u_out) { 1351 // Update _parent, _pred, _pred_dir 1352 _parent[u_in] = v_in; 1353 _pred[u_in] = in_arc; 1354 _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN; 1355 1356 // Update _thread and _rev_thread 1357 if (_thread[v_in] != u_out) { 1358 int after = _thread[old_last_succ]; 1359 _thread[old_rev_thread] = after; 1360 _rev_thread[after] = old_rev_thread; 1361 after = _thread[v_in]; 1362 _thread[v_in] = u_out; 1363 _rev_thread[u_out] = v_in; 1364 _thread[old_last_succ] = after; 1365 _rev_thread[after] = old_last_succ; 1366 } 1367 } else { 1368 // Handle the case when old_rev_thread equals to v_in 1369 // (it also means that join and v_out coincide) 1370 int thread_continue = old_rev_thread == v_in ? 1371 _thread[old_last_succ] : _thread[v_in]; 1372 1373 // Update _thread and _parent along the stem nodes (i.e. the nodes 1374 // between u_in and u_out, whose parent have to be changed) 1375 int stem = u_in; // the current stem node 1376 int par_stem = v_in; // the new parent of stem 1377 int next_stem; // the next stem node 1378 int last = _last_succ[u_in]; // the last successor of stem 1379 int before, after = _thread[last]; 1380 _thread[v_in] = u_in; 1381 _dirty_revs.clear(); 1382 _dirty_revs.push_back(v_in); 1383 while (stem != u_out) { 1384 // Insert the next stem node into the thread list 1385 next_stem = _parent[stem]; 1386 _thread[last] = next_stem; 1387 _dirty_revs.push_back(last); 1388 1389 // Remove the subtree of stem from the thread list 1390 before = _rev_thread[stem]; 1391 _thread[before] = after; 1392 _rev_thread[after] = before; 1393 1394 // Change the parent node and shift stem nodes 1395 _parent[stem] = par_stem; 1396 par_stem = stem; 1397 stem = next_stem; 1398 1399 // Update last and after 1400 last = _last_succ[stem] == _last_succ[par_stem] ? 1401 _rev_thread[par_stem] : _last_succ[stem]; 1402 after = _thread[last]; 1403 } 1404 _parent[u_out] = par_stem; 1405 _thread[last] = thread_continue; 1406 _rev_thread[thread_continue] = last; 1407 _last_succ[u_out] = last; 1408 1409 // Remove the subtree of u_out from the thread list except for 1410 // the case when old_rev_thread equals to v_in 1411 if (old_rev_thread != v_in) { 1412 _thread[old_rev_thread] = after; 1413 _rev_thread[after] = old_rev_thread; 1414 } 1415 1416 // Update _rev_thread using the new _thread values 1417 for (int i = 0; i != int(_dirty_revs.size()); ++i) { 1418 int u = _dirty_revs[i]; 1419 _rev_thread[_thread[u]] = u; 1420 } 1421 1422 // Update _pred, _pred_dir, _last_succ and _succ_num for the 1423 // stem nodes from u_out to u_in 1424 int tmp_sc = 0, tmp_ls = _last_succ[u_out]; 1425 for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) { 1426 _pred[u] = _pred[p]; 1427 _pred_dir[u] = -_pred_dir[p]; 1428 tmp_sc += _succ_num[u] - _succ_num[p]; 1429 _succ_num[u] = tmp_sc; 1430 _last_succ[p] = tmp_ls; 1431 } 1432 _pred[u_in] = in_arc; 1433 _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN; 1434 _succ_num[u_in] = old_succ_num; 1435 } 1436 1437 // Update _last_succ from v_in towards the root 1438 int up_limit_out = _last_succ[join] == v_in ? join : -1; 1439 int last_succ_out = _last_succ[u_out]; 1440 for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) { 1441 _last_succ[u] = last_succ_out; 1442 } 1443 1444 // Update _last_succ from v_out towards the root 1445 if (join != old_rev_thread && v_in != old_rev_thread) { 1446 for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; 1447 u = _parent[u]) { 1448 _last_succ[u] = old_rev_thread; 1449 } 1450 } 1451 else if (last_succ_out != old_last_succ) { 1452 for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; 1453 u = _parent[u]) { 1454 _last_succ[u] = last_succ_out; 1455 } 1456 } 1457 1458 // Update _succ_num from v_in to join 1459 for (int u = v_in; u != join; u = _parent[u]) { 1460 _succ_num[u] += old_succ_num; 1461 } 1462 // Update _succ_num from v_out to join 1463 for (int u = v_out; u != join; u = _parent[u]) { 1464 _succ_num[u] -= old_succ_num; 1465 } 1466 } 1467 1468 // Update potentials in the subtree that has been moved updatePotential()1469 void updatePotential() { 1470 Cost sigma = _pi[v_in] - _pi[u_in] - 1471 _pred_dir[u_in] * _cost[in_arc]; 1472 int end = _thread[_last_succ[u_in]]; 1473 for (int u = u_in; u != end; u = _thread[u]) { 1474 _pi[u] += sigma; 1475 } 1476 } 1477 1478 // Heuristic initial pivots initialPivots()1479 bool initialPivots() { 1480 Value curr, total = 0; 1481 std::vector<Node> supply_nodes, demand_nodes; 1482 for (NodeIt u(_graph); u != INVALID; ++u) { 1483 curr = _supply[_node_id[u]]; 1484 if (curr > 0) { 1485 total += curr; 1486 supply_nodes.push_back(u); 1487 } 1488 else if (curr < 0) { 1489 demand_nodes.push_back(u); 1490 } 1491 } 1492 if (_sum_supply > 0) total -= _sum_supply; 1493 if (total <= 0) return true; 1494 1495 IntVector arc_vector; 1496 if (_sum_supply >= 0) { 1497 if (supply_nodes.size() == 1 && demand_nodes.size() == 1) { 1498 // Perform a reverse graph search from the sink to the source 1499 typename GR::template NodeMap<bool> reached(_graph, false); 1500 Node s = supply_nodes[0], t = demand_nodes[0]; 1501 std::vector<Node> stack; 1502 reached[t] = true; 1503 stack.push_back(t); 1504 while (!stack.empty()) { 1505 Node u, v = stack.back(); 1506 stack.pop_back(); 1507 if (v == s) break; 1508 for (InArcIt a(_graph, v); a != INVALID; ++a) { 1509 if (reached[u = _graph.source(a)]) continue; 1510 int j = _arc_id[a]; 1511 if (_cap[j] >= total) { 1512 arc_vector.push_back(j); 1513 reached[u] = true; 1514 stack.push_back(u); 1515 } 1516 } 1517 } 1518 } else { 1519 // Find the min. cost incoming arc for each demand node 1520 for (int i = 0; i != int(demand_nodes.size()); ++i) { 1521 Node v = demand_nodes[i]; 1522 Cost c, min_cost = std::numeric_limits<Cost>::max(); 1523 Arc min_arc = INVALID; 1524 for (InArcIt a(_graph, v); a != INVALID; ++a) { 1525 c = _cost[_arc_id[a]]; 1526 if (c < min_cost) { 1527 min_cost = c; 1528 min_arc = a; 1529 } 1530 } 1531 if (min_arc != INVALID) { 1532 arc_vector.push_back(_arc_id[min_arc]); 1533 } 1534 } 1535 } 1536 } else { 1537 // Find the min. cost outgoing arc for each supply node 1538 for (int i = 0; i != int(supply_nodes.size()); ++i) { 1539 Node u = supply_nodes[i]; 1540 Cost c, min_cost = std::numeric_limits<Cost>::max(); 1541 Arc min_arc = INVALID; 1542 for (OutArcIt a(_graph, u); a != INVALID; ++a) { 1543 c = _cost[_arc_id[a]]; 1544 if (c < min_cost) { 1545 min_cost = c; 1546 min_arc = a; 1547 } 1548 } 1549 if (min_arc != INVALID) { 1550 arc_vector.push_back(_arc_id[min_arc]); 1551 } 1552 } 1553 } 1554 1555 // Perform heuristic initial pivots 1556 for (int i = 0; i != int(arc_vector.size()); ++i) { 1557 in_arc = arc_vector[i]; 1558 if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] - 1559 _pi[_target[in_arc]]) >= 0) continue; 1560 findJoinNode(); 1561 bool change = findLeavingArc(); 1562 if (delta >= MAX_VALUE) return false; 1563 changeFlow(change); 1564 if (change) { 1565 updateTreeStructure(); 1566 updatePotential(); 1567 } 1568 } 1569 return true; 1570 } 1571 1572 // Execute the algorithm start(PivotRule pivot_rule)1573 ProblemType start(PivotRule pivot_rule) { 1574 // Select the pivot rule implementation 1575 switch (pivot_rule) { 1576 case FIRST_ELIGIBLE: 1577 return start<FirstEligiblePivotRule>(); 1578 case BEST_ELIGIBLE: 1579 return start<BestEligiblePivotRule>(); 1580 case BLOCK_SEARCH: 1581 return start<BlockSearchPivotRule>(); 1582 case CANDIDATE_LIST: 1583 return start<CandidateListPivotRule>(); 1584 case ALTERING_LIST: 1585 return start<AlteringListPivotRule>(); 1586 } 1587 return INFEASIBLE; // avoid warning 1588 } 1589 1590 template <typename PivotRuleImpl> start()1591 ProblemType start() { 1592 PivotRuleImpl pivot(*this); 1593 1594 // Perform heuristic initial pivots 1595 if (!initialPivots()) return UNBOUNDED; 1596 1597 // Execute the Network Simplex algorithm 1598 while (pivot.findEnteringArc()) { 1599 findJoinNode(); 1600 bool change = findLeavingArc(); 1601 if (delta >= MAX_VALUE) return UNBOUNDED; 1602 changeFlow(change); 1603 if (change) { 1604 updateTreeStructure(); 1605 updatePotential(); 1606 } 1607 } 1608 1609 // Check feasibility 1610 for (int e = _search_arc_num; e != _all_arc_num; ++e) { 1611 if (_flow[e] != 0) return INFEASIBLE; 1612 } 1613 1614 // Transform the solution and the supply map to the original form 1615 if (_has_lower) { 1616 for (int i = 0; i != _arc_num; ++i) { 1617 Value c = _lower[i]; 1618 if (c != 0) { 1619 _flow[i] += c; 1620 _supply[_source[i]] += c; 1621 _supply[_target[i]] -= c; 1622 } 1623 } 1624 } 1625 1626 // Shift potentials to meet the requirements of the GEQ/LEQ type 1627 // optimality conditions 1628 if (_sum_supply == 0) { 1629 if (_stype == GEQ) { 1630 Cost max_pot = -std::numeric_limits<Cost>::max(); 1631 for (int i = 0; i != _node_num; ++i) { 1632 if (_pi[i] > max_pot) max_pot = _pi[i]; 1633 } 1634 if (max_pot > 0) { 1635 for (int i = 0; i != _node_num; ++i) 1636 _pi[i] -= max_pot; 1637 } 1638 } else { 1639 Cost min_pot = std::numeric_limits<Cost>::max(); 1640 for (int i = 0; i != _node_num; ++i) { 1641 if (_pi[i] < min_pot) min_pot = _pi[i]; 1642 } 1643 if (min_pot < 0) { 1644 for (int i = 0; i != _node_num; ++i) 1645 _pi[i] -= min_pot; 1646 } 1647 } 1648 } 1649 1650 return OPTIMAL; 1651 } 1652 1653 }; //class NetworkSimplex 1654 1655 ///@} 1656 1657 } //namespace lemon 1658 1659 #endif //LEMON_NETWORK_SIMPLEX_H 1660