Changeset 979:0bca98cbebbb in lemon for lemon
 Timestamp:
 05/03/10 10:56:24 (10 years ago)
 Branch:
 default
 Parents:
 975:ca4059d63236 (diff), 976:5205145fabf6 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.  Phase:
 public
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

lemon/network_simplex.h
r956 r978 1078 1078 ART_COST = std::numeric_limits<Cost>::max() / 2 + 1; 1079 1079 } else { 1080 ART_COST = std::numeric_limits<Cost>::min();1080 ART_COST = 0; 1081 1081 for (int i = 0; i != _arc_num; ++i) { 1082 1082 if (_cost[i] > ART_COST) ART_COST = _cost[i]; … … 1590 1590 if (_sum_supply == 0) { 1591 1591 if (_stype == GEQ) { 1592 Cost max_pot = std::numeric_limits<Cost>::min();1592 Cost max_pot = std::numeric_limits<Cost>::max(); 1593 1593 for (int i = 0; i != _node_num; ++i) { 1594 1594 if (_pi[i] > max_pot) max_pot = _pi[i]; 
lemon/network_simplex.h
r976 r978 3 3 * This file is a part of LEMON, a generic C++ optimization library. 4 4 * 5 * Copyright (C) 200320 095 * Copyright (C) 20032010 6 6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport 7 7 * (Egervary Research Group on Combinatorial Optimization, EGRES). … … 41 41 /// 42 42 /// \ref NetworkSimplex implements the primal Network Simplex algorithm 43 /// for finding a \ref min_cost_flow "minimum cost flow". 44 /// This algorithm is a specialized version of the linear programming 45 /// simplex method directly for the minimum cost flow problem. 46 /// It is one of the most efficient solution methods. 43 /// for finding a \ref min_cost_flow "minimum cost flow" 44 /// \ref amo93networkflows, \ref dantzig63linearprog, 45 /// \ref 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. 47 49 /// 48 /// In general this classis the fastest implementation available49 /// in LEMON for th e minimum cost flowproblem.50 /// Moreover it supports both directions of the supply/demand inequality51 /// constraints. For more information see \ref SupplyType.50 /// In general, %NetworkSimplex is the fastest implementation available 51 /// in LEMON for this problem. 52 /// Moreover, it supports both directions of the supply/demand inequality 53 /// constraints. For more information, see \ref SupplyType. 52 54 /// 53 55 /// Most of the parameters of the problem (except for the digraph) … … 57 59 /// 58 60 /// \tparam GR The digraph type the algorithm runs on. 59 /// \tparam V The valuetype used for flow amounts, capacity bounds60 /// and supply values in the algorithm. By default it is \c int.61 /// \tparam C The valuetype used for costs and potentials in the62 /// algorithm. By default it is the same as \c V.61 /// \tparam V The number type used for flow amounts, capacity bounds 62 /// and supply values in the algorithm. By default, it is \c int. 63 /// \tparam C The number type used for costs and potentials in the 64 /// algorithm. By default, it is the same as \c V. 63 65 /// 64 /// \warning Both valuetypes must be signed and all input data must66 /// \warning Both number types must be signed and all input data must 65 67 /// be integer. 66 68 /// 67 69 /// \note %NetworkSimplex provides five different pivot rule 68 70 /// implementations, from which the most efficient one is used 69 /// by default. For more information see \ref PivotRule.71 /// by default. For more information, see \ref PivotRule. 70 72 template <typename GR, typename V = int, typename C = V> 71 73 class NetworkSimplex … … 96 98 UNBOUNDED 97 99 }; 98 100 99 101 /// \brief Constants for selecting the type of the supply constraints. 100 102 /// … … 114 116 LEQ 115 117 }; 116 118 117 119 /// \brief Constants for selecting the pivot rule. 118 120 /// … … 123 125 /// implementations that significantly affect the running time 124 126 /// of the algorithm. 125 /// By default \ref BLOCK_SEARCH "Block Search" is used, which127 /// By default, \ref BLOCK_SEARCH "Block Search" is used, which 126 128 /// proved to be the most efficient and the most robust on various 127 /// test inputs according to our benchmark tests.128 /// However another pivot rule can be selected using the \ref run()129 /// test inputs. 130 /// However, another pivot rule can be selected using the \ref run() 129 131 /// function with the proper parameter. 130 132 enum PivotRule { 131 133 132 /// The FirstEligible pivot rule.134 /// The \e First \e Eligible pivot rule. 133 135 /// The next eligible arc is selected in a wraparound fashion 134 136 /// in every iteration. 135 137 FIRST_ELIGIBLE, 136 138 137 /// The BestEligible pivot rule.139 /// The \e Best \e Eligible pivot rule. 138 140 /// The best eligible arc is selected in every iteration. 139 141 BEST_ELIGIBLE, 140 142 141 /// The BlockSearch pivot rule.143 /// The \e Block \e Search pivot rule. 142 144 /// A specified number of arcs are examined in every iteration 143 145 /// in a wraparound fashion and the best eligible arc is selected … … 145 147 BLOCK_SEARCH, 146 148 147 /// The Candidate List pivot rule.149 /// The \e Candidate \e List pivot rule. 148 150 /// In a major iteration a candidate list is built from eligible arcs 149 151 /// in a wraparound fashion and in the following minor iterations … … 151 153 CANDIDATE_LIST, 152 154 153 /// The Altering Candidate List pivot rule.155 /// The \e Altering \e Candidate \e List pivot rule. 154 156 /// It is a modified version of the Candidate List method. 155 157 /// It keeps only the several best eligible arcs from the former … … 157 159 ALTERING_LIST 158 160 }; 159 161 160 162 private: 161 163 162 164 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 163 165 164 typedef std::vector<Arc> ArcVector;165 typedef std::vector<Node> NodeVector;166 166 typedef std::vector<int> IntVector; 167 typedef std::vector<bool> BoolVector;168 167 typedef std::vector<Value> ValueVector; 169 168 typedef std::vector<Cost> CostVector; 169 typedef std::vector<char> BoolVector; 170 // Note: vector<char> is used instead of vector<bool> for efficiency reasons 170 171 171 172 // State constants for arcs 172 enum ArcState Enum{173 enum ArcState { 173 174 STATE_UPPER = 1, 174 175 STATE_TREE = 0, 175 176 STATE_LOWER = 1 176 177 }; 178 179 typedef std::vector<signed char> StateVector; 180 // Note: vector<signed char> is used instead of vector<ArcState> for 181 // efficiency reasons 177 182 178 183 private: … … 195 200 IntVector _source; 196 201 IntVector _target; 202 bool _arc_mixing; 197 203 198 204 // Node and arc data … … 214 220 IntVector _dirty_revs; 215 221 BoolVector _forward; 216 IntVector _state;222 StateVector _state; 217 223 int _root; 218 224 … … 223 229 Value delta; 224 230 231 const Value MAX; 232 225 233 public: 226 234 227 235 /// \brief Constant for infinite upper bounds (capacities). 228 236 /// … … 243 251 const IntVector &_target; 244 252 const CostVector &_cost; 245 const IntVector&_state;253 const StateVector &_state; 246 254 const CostVector &_pi; 247 255 int &_in_arc; … … 264 272 bool findEnteringArc() { 265 273 Cost c; 266 for (int e = _next_arc; e <_search_arc_num; ++e) {274 for (int e = _next_arc; e != _search_arc_num; ++e) { 267 275 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 268 276 if (c < 0) { … … 272 280 } 273 281 } 274 for (int e = 0; e <_next_arc; ++e) {282 for (int e = 0; e != _next_arc; ++e) { 275 283 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 276 284 if (c < 0) { … … 295 303 const IntVector &_target; 296 304 const CostVector &_cost; 297 const IntVector&_state;305 const StateVector &_state; 298 306 const CostVector &_pi; 299 307 int &_in_arc; … … 312 320 bool findEnteringArc() { 313 321 Cost c, min = 0; 314 for (int e = 0; e <_search_arc_num; ++e) {322 for (int e = 0; e != _search_arc_num; ++e) { 315 323 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 316 324 if (c < min) { … … 334 342 const IntVector &_target; 335 343 const CostVector &_cost; 336 const IntVector&_state;344 const StateVector &_state; 337 345 const CostVector &_pi; 338 346 int &_in_arc; … … 353 361 { 354 362 // The main parameters of the pivot rule 355 const double BLOCK_SIZE_FACTOR = 0.5;363 const double BLOCK_SIZE_FACTOR = 1.0; 356 364 const int MIN_BLOCK_SIZE = 10; 357 365 … … 365 373 Cost c, min = 0; 366 374 int cnt = _block_size; 367 int e , min_arc = _next_arc;368 for (e = _next_arc; e <_search_arc_num; ++e) {375 int e; 376 for (e = _next_arc; e != _search_arc_num; ++e) { 369 377 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 370 378 if (c < min) { 371 379 min = c; 372 min_arc = e;380 _in_arc = e; 373 381 } 374 382 if (cnt == 0) { 375 if (min < 0) break;383 if (min < 0) goto search_end; 376 384 cnt = _block_size; 377 385 } 378 386 } 379 if (min == 0  cnt > 0) { 380 for (e = 0; e < _next_arc; ++e) { 381 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 382 if (c < min) { 383 min = c; 384 min_arc = e; 385 } 386 if (cnt == 0) { 387 if (min < 0) break; 388 cnt = _block_size; 389 } 387 for (e = 0; e != _next_arc; ++e) { 388 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 389 if (c < min) { 390 min = c; 391 _in_arc = e; 392 } 393 if (cnt == 0) { 394 if (min < 0) goto search_end; 395 cnt = _block_size; 390 396 } 391 397 } 392 398 if (min >= 0) return false; 393 _in_arc = min_arc; 399 400 search_end: 394 401 _next_arc = e; 395 402 return true; … … 408 415 const IntVector &_target; 409 416 const CostVector &_cost; 410 const IntVector&_state;417 const StateVector &_state; 411 418 const CostVector &_pi; 412 419 int &_in_arc; … … 429 436 { 430 437 // The main parameters of the pivot rule 431 const double LIST_LENGTH_FACTOR = 1.0;438 const double LIST_LENGTH_FACTOR = 0.25; 432 439 const int MIN_LIST_LENGTH = 10; 433 440 const double MINOR_LIMIT_FACTOR = 0.1; … … 446 453 bool findEnteringArc() { 447 454 Cost min, c; 448 int e , min_arc = _next_arc;455 int e; 449 456 if (_curr_length > 0 && _minor_count < _minor_limit) { 450 457 // Minor iteration: select the best eligible arc from the … … 457 464 if (c < min) { 458 465 min = c; 459 min_arc = e;466 _in_arc = e; 460 467 } 461 if (c >= 0) {468 else if (c >= 0) { 462 469 _candidates[i] = _candidates[_curr_length]; 463 470 } 464 471 } 465 if (min < 0) { 466 _in_arc = min_arc; 467 return true; 468 } 472 if (min < 0) return true; 469 473 } 470 474 … … 472 476 min = 0; 473 477 _curr_length = 0; 474 for (e = _next_arc; e <_search_arc_num; ++e) {478 for (e = _next_arc; e != _search_arc_num; ++e) { 475 479 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 476 480 if (c < 0) { … … 478 482 if (c < min) { 479 483 min = c; 480 min_arc = e;484 _in_arc = e; 481 485 } 482 if (_curr_length == _list_length) break; 483 } 484 } 485 if (_curr_length < _list_length) { 486 for (e = 0; e < _next_arc; ++e) { 487 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 488 if (c < 0) { 489 _candidates[_curr_length++] = e; 490 if (c < min) { 491 min = c; 492 min_arc = e; 493 } 494 if (_curr_length == _list_length) break; 486 if (_curr_length == _list_length) goto search_end; 487 } 488 } 489 for (e = 0; e != _next_arc; ++e) { 490 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 491 if (c < 0) { 492 _candidates[_curr_length++] = e; 493 if (c < min) { 494 min = c; 495 _in_arc = e; 495 496 } 497 if (_curr_length == _list_length) goto search_end; 496 498 } 497 499 } 498 500 if (_curr_length == 0) return false; 501 502 search_end: 499 503 _minor_count = 1; 500 _in_arc = min_arc;501 504 _next_arc = e; 502 505 return true; … … 515 518 const IntVector &_target; 516 519 const CostVector &_cost; 517 const IntVector&_state;520 const StateVector &_state; 518 521 const CostVector &_pi; 519 522 int &_in_arc; … … 550 553 { 551 554 // The main parameters of the pivot rule 552 const double BLOCK_SIZE_FACTOR = 1. 5;555 const double BLOCK_SIZE_FACTOR = 1.0; 553 556 const int MIN_BLOCK_SIZE = 10; 554 557 const double HEAD_LENGTH_FACTOR = 0.1; … … 568 571 // Check the current candidate list 569 572 int e; 570 for (int i = 0; i <_curr_length; ++i) {573 for (int i = 0; i != _curr_length; ++i) { 571 574 e = _candidates[i]; 572 575 _cand_cost[e] = _state[e] * … … 579 582 // Extend the list 580 583 int cnt = _block_size; 581 int last_arc = 0;582 584 int limit = _head_length; 583 585 584 for ( int e = _next_arc; e <_search_arc_num; ++e) {586 for (e = _next_arc; e != _search_arc_num; ++e) { 585 587 _cand_cost[e] = _state[e] * 586 588 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 587 589 if (_cand_cost[e] < 0) { 588 590 _candidates[_curr_length++] = e; 589 last_arc = e;590 591 } 591 592 if (cnt == 0) { 592 if (_curr_length > limit) break;593 if (_curr_length > limit) goto search_end; 593 594 limit = 0; 594 595 cnt = _block_size; 595 596 } 596 597 } 597 if (_curr_length <= limit) { 598 for (int e = 0; e < _next_arc; ++e) { 599 _cand_cost[e] = _state[e] * 600 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 601 if (_cand_cost[e] < 0) { 602 _candidates[_curr_length++] = e; 603 last_arc = e; 604 } 605 if (cnt == 0) { 606 if (_curr_length > limit) break; 607 limit = 0; 608 cnt = _block_size; 609 } 598 for (e = 0; e != _next_arc; ++e) { 599 _cand_cost[e] = _state[e] * 600 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 601 if (_cand_cost[e] < 0) { 602 _candidates[_curr_length++] = e; 603 } 604 if (cnt == 0) { 605 if (_curr_length > limit) goto search_end; 606 limit = 0; 607 cnt = _block_size; 610 608 } 611 609 } 612 610 if (_curr_length == 0) return false; 613 _next_arc = last_arc + 1; 611 612 search_end: 614 613 615 614 // Make heap of the candidate list (approximating a partial sort) … … 619 618 // Pop the first element of the heap 620 619 _in_arc = _candidates[0]; 620 _next_arc = e; 621 621 pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, 622 622 _sort_func ); … … 634 634 /// 635 635 /// \param graph The digraph the algorithm runs on. 636 NetworkSimplex(const GR& graph) : 636 /// \param arc_mixing Indicate if the arcs have to be stored in a 637 /// mixed order in the internal data structure. 638 /// In special cases, it could lead to better overall performance, 639 /// but it is usually slower. Therefore it is disabled by default. 640 NetworkSimplex(const GR& graph, bool arc_mixing = false) : 637 641 _graph(graph), _node_id(graph), _arc_id(graph), 642 _arc_mixing(arc_mixing), 643 MAX(std::numeric_limits<Value>::max()), 638 644 INF(std::numeric_limits<Value>::has_infinity ? 639 std::numeric_limits<Value>::infinity() : 640 std::numeric_limits<Value>::max()) 645 std::numeric_limits<Value>::infinity() : MAX) 641 646 { 642 // Check the valuetypes647 // Check the number types 643 648 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, 644 649 "The flow type of NetworkSimplex must be signed"); 645 650 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 646 651 "The cost type of NetworkSimplex must be signed"); 647 652 653 // Reset data structures 654 reset(); 655 } 656 657 /// \name Parameters 658 /// The parameters of the algorithm can be specified using these 659 /// functions. 660 661 /// @{ 662 663 /// \brief Set the lower bounds on the arcs. 664 /// 665 /// This function sets the lower bounds on the arcs. 666 /// If it is not used before calling \ref run(), the lower bounds 667 /// will be set to zero on all arcs. 668 /// 669 /// \param map An arc map storing the lower bounds. 670 /// Its \c Value type must be convertible to the \c Value type 671 /// of the algorithm. 672 /// 673 /// \return <tt>(*this)</tt> 674 template <typename LowerMap> 675 NetworkSimplex& lowerMap(const LowerMap& map) { 676 _have_lower = true; 677 for (ArcIt a(_graph); a != INVALID; ++a) { 678 _lower[_arc_id[a]] = map[a]; 679 } 680 return *this; 681 } 682 683 /// \brief Set the upper bounds (capacities) on the arcs. 684 /// 685 /// This function sets the upper bounds (capacities) on the arcs. 686 /// If it is not used before calling \ref run(), the upper bounds 687 /// will be set to \ref INF on all arcs (i.e. the flow value will be 688 /// unbounded from above). 689 /// 690 /// \param map An arc map storing the upper bounds. 691 /// Its \c Value type must be convertible to the \c Value type 692 /// of the algorithm. 693 /// 694 /// \return <tt>(*this)</tt> 695 template<typename UpperMap> 696 NetworkSimplex& upperMap(const UpperMap& map) { 697 for (ArcIt a(_graph); a != INVALID; ++a) { 698 _upper[_arc_id[a]] = map[a]; 699 } 700 return *this; 701 } 702 703 /// \brief Set the costs of the arcs. 704 /// 705 /// This function sets the costs of the arcs. 706 /// If it is not used before calling \ref run(), the costs 707 /// will be set to \c 1 on all arcs. 708 /// 709 /// \param map An arc map storing the costs. 710 /// Its \c Value type must be convertible to the \c Cost type 711 /// of the algorithm. 712 /// 713 /// \return <tt>(*this)</tt> 714 template<typename CostMap> 715 NetworkSimplex& costMap(const CostMap& map) { 716 for (ArcIt a(_graph); a != INVALID; ++a) { 717 _cost[_arc_id[a]] = map[a]; 718 } 719 return *this; 720 } 721 722 /// \brief Set the supply values of the nodes. 723 /// 724 /// This function sets the supply values of the nodes. 725 /// If neither this function nor \ref stSupply() is used before 726 /// calling \ref run(), the supply of each node will be set to zero. 727 /// 728 /// \param map A node map storing the supply values. 729 /// Its \c Value type must be convertible to the \c Value type 730 /// of the algorithm. 731 /// 732 /// \return <tt>(*this)</tt> 733 template<typename SupplyMap> 734 NetworkSimplex& supplyMap(const SupplyMap& map) { 735 for (NodeIt n(_graph); n != INVALID; ++n) { 736 _supply[_node_id[n]] = map[n]; 737 } 738 return *this; 739 } 740 741 /// \brief Set single source and target nodes and a supply value. 742 /// 743 /// This function sets a single source node and a single target node 744 /// and the required flow value. 745 /// If neither this function nor \ref supplyMap() is used before 746 /// calling \ref run(), the supply of each node will be set to zero. 747 /// 748 /// Using this function has the same effect as using \ref supplyMap() 749 /// with such a map in which \c k is assigned to \c s, \c k is 750 /// assigned to \c t and all other nodes have zero supply value. 751 /// 752 /// \param s The source node. 753 /// \param t The target node. 754 /// \param k The required amount of flow from node \c s to node \c t 755 /// (i.e. the supply of \c s and the demand of \c t). 756 /// 757 /// \return <tt>(*this)</tt> 758 NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { 759 for (int i = 0; i != _node_num; ++i) { 760 _supply[i] = 0; 761 } 762 _supply[_node_id[s]] = k; 763 _supply[_node_id[t]] = k; 764 return *this; 765 } 766 767 /// \brief Set the type of the supply constraints. 768 /// 769 /// This function sets the type of the supply/demand constraints. 770 /// If it is not used before calling \ref run(), the \ref GEQ supply 771 /// type will be used. 772 /// 773 /// For more information, see \ref SupplyType. 774 /// 775 /// \return <tt>(*this)</tt> 776 NetworkSimplex& supplyType(SupplyType supply_type) { 777 _stype = supply_type; 778 return *this; 779 } 780 781 /// @} 782 783 /// \name Execution Control 784 /// The algorithm can be executed using \ref run(). 785 786 /// @{ 787 788 /// \brief Run the algorithm. 789 /// 790 /// This function runs the algorithm. 791 /// The paramters can be specified using functions \ref lowerMap(), 792 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 793 /// \ref supplyType(). 794 /// For example, 795 /// \code 796 /// NetworkSimplex<ListDigraph> ns(graph); 797 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 798 /// .supplyMap(sup).run(); 799 /// \endcode 800 /// 801 /// This function can be called more than once. All the given parameters 802 /// are kept for the next call, unless \ref resetParams() or \ref reset() 803 /// is used, thus only the modified parameters have to be set again. 804 /// If the underlying digraph was also modified after the construction 805 /// of the class (or the last \ref reset() call), then the \ref reset() 806 /// function must be called. 807 /// 808 /// \param pivot_rule The pivot rule that will be used during the 809 /// algorithm. For more information, see \ref PivotRule. 810 /// 811 /// \return \c INFEASIBLE if no feasible flow exists, 812 /// \n \c OPTIMAL if the problem has optimal solution 813 /// (i.e. it is feasible and bounded), and the algorithm has found 814 /// optimal flow and node potentials (primal and dual solutions), 815 /// \n \c UNBOUNDED if the objective function of the problem is 816 /// unbounded, i.e. there is a directed cycle having negative total 817 /// cost and infinite upper bound. 818 /// 819 /// \see ProblemType, PivotRule 820 /// \see resetParams(), reset() 821 ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { 822 if (!init()) return INFEASIBLE; 823 return start(pivot_rule); 824 } 825 826 /// \brief Reset all the parameters that have been given before. 827 /// 828 /// This function resets all the paramaters that have been given 829 /// before using functions \ref lowerMap(), \ref upperMap(), 830 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). 831 /// 832 /// It is useful for multiple \ref run() calls. Basically, all the given 833 /// parameters are kept for the next \ref run() call, unless 834 /// \ref resetParams() or \ref reset() is used. 835 /// If the underlying digraph was also modified after the construction 836 /// of the class or the last \ref reset() call, then the \ref reset() 837 /// function must be used, otherwise \ref resetParams() is sufficient. 838 /// 839 /// For example, 840 /// \code 841 /// NetworkSimplex<ListDigraph> ns(graph); 842 /// 843 /// // First run 844 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 845 /// .supplyMap(sup).run(); 846 /// 847 /// // Run again with modified cost map (resetParams() is not called, 848 /// // so only the cost map have to be set again) 849 /// cost[e] += 100; 850 /// ns.costMap(cost).run(); 851 /// 852 /// // Run again from scratch using resetParams() 853 /// // (the lower bounds will be set to zero on all arcs) 854 /// ns.resetParams(); 855 /// ns.upperMap(capacity).costMap(cost) 856 /// .supplyMap(sup).run(); 857 /// \endcode 858 /// 859 /// \return <tt>(*this)</tt> 860 /// 861 /// \see reset(), run() 862 NetworkSimplex& resetParams() { 863 for (int i = 0; i != _node_num; ++i) { 864 _supply[i] = 0; 865 } 866 for (int i = 0; i != _arc_num; ++i) { 867 _lower[i] = 0; 868 _upper[i] = INF; 869 _cost[i] = 1; 870 } 871 _have_lower = false; 872 _stype = GEQ; 873 return *this; 874 } 875 876 /// \brief Reset the internal data structures and all the parameters 877 /// that have been given before. 878 /// 879 /// This function resets the internal data structures and all the 880 /// paramaters that have been given before using functions \ref lowerMap(), 881 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 882 /// \ref supplyType(). 883 /// 884 /// It is useful for multiple \ref run() calls. Basically, all the given 885 /// parameters are kept for the next \ref run() call, unless 886 /// \ref resetParams() or \ref reset() is used. 887 /// If the underlying digraph was also modified after the construction 888 /// of the class or the last \ref reset() call, then the \ref reset() 889 /// function must be used, otherwise \ref resetParams() is sufficient. 890 /// 891 /// See \ref resetParams() for examples. 892 /// 893 /// \return <tt>(*this)</tt> 894 /// 895 /// \see resetParams(), run() 896 NetworkSimplex& reset() { 648 897 // Resize vectors 649 898 _node_num = countNodes(_graph); … … 672 921 _state.resize(max_arc_num); 673 922 674 // Copy the graph (store the arcs in a mixed order)923 // Copy the graph 675 924 int i = 0; 676 925 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 677 926 _node_id[n] = i; 678 927 } 679 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 680 i = 0; 681 for (ArcIt a(_graph); a != INVALID; ++a) { 682 _arc_id[a] = i; 683 _source[i] = _node_id[_graph.source(a)]; 684 _target[i] = _node_id[_graph.target(a)]; 685 if ((i += k) >= _arc_num) i = (i % k) + 1; 686 } 687 688 // Initialize maps 689 for (int i = 0; i != _node_num; ++i) { 690 _supply[i] = 0; 691 } 692 for (int i = 0; i != _arc_num; ++i) { 693 _lower[i] = 0; 694 _upper[i] = INF; 695 _cost[i] = 1; 696 } 697 _have_lower = false; 698 _stype = GEQ; 699 } 700 701 /// \name Parameters 702 /// The parameters of the algorithm can be specified using these 703 /// functions. 704 705 /// @{ 706 707 /// \brief Set the lower bounds on the arcs. 708 /// 709 /// This function sets the lower bounds on the arcs. 710 /// If it is not used before calling \ref run(), the lower bounds 711 /// will be set to zero on all arcs. 712 /// 713 /// \param map An arc map storing the lower bounds. 714 /// Its \c Value type must be convertible to the \c Value type 715 /// of the algorithm. 716 /// 717 /// \return <tt>(*this)</tt> 718 template <typename LowerMap> 719 NetworkSimplex& lowerMap(const LowerMap& map) { 720 _have_lower = true; 721 for (ArcIt a(_graph); a != INVALID; ++a) { 722 _lower[_arc_id[a]] = map[a]; 723 } 724 return *this; 725 } 726 727 /// \brief Set the upper bounds (capacities) on the arcs. 728 /// 729 /// This function sets the upper bounds (capacities) on the arcs. 730 /// If it is not used before calling \ref run(), the upper bounds 731 /// will be set to \ref INF on all arcs (i.e. the flow value will be 732 /// unbounded from above on each arc). 733 /// 734 /// \param map An arc map storing the upper bounds. 735 /// Its \c Value type must be convertible to the \c Value type 736 /// of the algorithm. 737 /// 738 /// \return <tt>(*this)</tt> 739 template<typename UpperMap> 740 NetworkSimplex& upperMap(const UpperMap& map) { 741 for (ArcIt a(_graph); a != INVALID; ++a) { 742 _upper[_arc_id[a]] = map[a]; 743 } 744 return *this; 745 } 746 747 /// \brief Set the costs of the arcs. 748 /// 749 /// This function sets the costs of the arcs. 750 /// If it is not used before calling \ref run(), the costs 751 /// will be set to \c 1 on all arcs. 752 /// 753 /// \param map An arc map storing the costs. 754 /// Its \c Value type must be convertible to the \c Cost type 755 /// of the algorithm. 756 /// 757 /// \return <tt>(*this)</tt> 758 template<typename CostMap> 759 NetworkSimplex& costMap(const CostMap& map) { 760 for (ArcIt a(_graph); a != INVALID; ++a) { 761 _cost[_arc_id[a]] = map[a]; 762 } 763 return *this; 764 } 765 766 /// \brief Set the supply values of the nodes. 767 /// 768 /// This function sets the supply values of the nodes. 769 /// If neither this function nor \ref stSupply() is used before 770 /// calling \ref run(), the supply of each node will be set to zero. 771 /// (It makes sense only if nonzero lower bounds are given.) 772 /// 773 /// \param map A node map storing the supply values. 774 /// Its \c Value type must be convertible to the \c Value type 775 /// of the algorithm. 776 /// 777 /// \return <tt>(*this)</tt> 778 template<typename SupplyMap> 779 NetworkSimplex& supplyMap(const SupplyMap& map) { 780 for (NodeIt n(_graph); n != INVALID; ++n) { 781 _supply[_node_id[n]] = map[n]; 782 } 783 return *this; 784 } 785 786 /// \brief Set single source and target nodes and a supply value. 787 /// 788 /// This function sets a single source node and a single target node 789 /// and the required flow value. 790 /// If neither this function nor \ref supplyMap() is used before 791 /// calling \ref run(), the supply of each node will be set to zero. 792 /// (It makes sense only if nonzero lower bounds are given.) 793 /// 794 /// Using this function has the same effect as using \ref supplyMap() 795 /// with such a map in which \c k is assigned to \c s, \c k is 796 /// assigned to \c t and all other nodes have zero supply value. 797 /// 798 /// \param s The source node. 799 /// \param t The target node. 800 /// \param k The required amount of flow from node \c s to node \c t 801 /// (i.e. the supply of \c s and the demand of \c t). 802 /// 803 /// \return <tt>(*this)</tt> 804 NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { 805 for (int i = 0; i != _node_num; ++i) { 806 _supply[i] = 0; 807 } 808 _supply[_node_id[s]] = k; 809 _supply[_node_id[t]] = k; 810 return *this; 811 } 812 813 /// \brief Set the type of the supply constraints. 814 /// 815 /// This function sets the type of the supply/demand constraints. 816 /// If it is not used before calling \ref run(), the \ref GEQ supply 817 /// type will be used. 818 /// 819 /// For more information see \ref SupplyType. 820 /// 821 /// \return <tt>(*this)</tt> 822 NetworkSimplex& supplyType(SupplyType supply_type) { 823 _stype = supply_type; 824 return *this; 825 } 826 827 /// @} 828 829 /// \name Execution Control 830 /// The algorithm can be executed using \ref run(). 831 832 /// @{ 833 834 /// \brief Run the algorithm. 835 /// 836 /// This function runs the algorithm. 837 /// The paramters can be specified using functions \ref lowerMap(), 838 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 839 /// \ref supplyType(). 840 /// For example, 841 /// \code 842 /// NetworkSimplex<ListDigraph> ns(graph); 843 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 844 /// .supplyMap(sup).run(); 845 /// \endcode 846 /// 847 /// This function can be called more than once. All the parameters 848 /// that have been given are kept for the next call, unless 849 /// \ref reset() is called, thus only the modified parameters 850 /// have to be set again. See \ref reset() for examples. 851 /// However the underlying digraph must not be modified after this 852 /// class have been constructed, since it copies and extends the graph. 853 /// 854 /// \param pivot_rule The pivot rule that will be used during the 855 /// algorithm. For more information see \ref PivotRule. 856 /// 857 /// \return \c INFEASIBLE if no feasible flow exists, 858 /// \n \c OPTIMAL if the problem has optimal solution 859 /// (i.e. it is feasible and bounded), and the algorithm has found 860 /// optimal flow and node potentials (primal and dual solutions), 861 /// \n \c UNBOUNDED if the objective function of the problem is 862 /// unbounded, i.e. there is a directed cycle having negative total 863 /// cost and infinite upper bound. 864 /// 865 /// \see ProblemType, PivotRule 866 ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { 867 if (!init()) return INFEASIBLE; 868 return start(pivot_rule); 869 } 870 871 /// \brief Reset all the parameters that have been given before. 872 /// 873 /// This function resets all the paramaters that have been given 874 /// before using functions \ref lowerMap(), \ref upperMap(), 875 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). 876 /// 877 /// It is useful for multiple run() calls. If this function is not 878 /// used, all the parameters given before are kept for the next 879 /// \ref run() call. 880 /// However the underlying digraph must not be modified after this 881 /// class have been constructed, since it copies and extends the graph. 882 /// 883 /// For example, 884 /// \code 885 /// NetworkSimplex<ListDigraph> ns(graph); 886 /// 887 /// // First run 888 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 889 /// .supplyMap(sup).run(); 890 /// 891 /// // Run again with modified cost map (reset() is not called, 892 /// // so only the cost map have to be set again) 893 /// cost[e] += 100; 894 /// ns.costMap(cost).run(); 895 /// 896 /// // Run again from scratch using reset() 897 /// // (the lower bounds will be set to zero on all arcs) 898 /// ns.reset(); 899 /// ns.upperMap(capacity).costMap(cost) 900 /// .supplyMap(sup).run(); 901 /// \endcode 902 /// 903 /// \return <tt>(*this)</tt> 904 NetworkSimplex& reset() { 905 for (int i = 0; i != _node_num; ++i) { 906 _supply[i] = 0; 907 } 908 for (int i = 0; i != _arc_num; ++i) { 909 _lower[i] = 0; 910 _upper[i] = INF; 911 _cost[i] = 1; 912 } 913 _have_lower = false; 914 _stype = GEQ; 928 if (_arc_mixing) { 929 // Store the arcs in a mixed order 930 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 931 int i = 0, j = 0; 932 for (ArcIt a(_graph); a != INVALID; ++a) { 933 _arc_id[a] = i; 934 _source[i] = _node_id[_graph.source(a)]; 935 _target[i] = _node_id[_graph.target(a)]; 936 if ((i += k) >= _arc_num) i = ++j; 937 } 938 } else { 939 // Store the arcs in the original order 940 int i = 0; 941 for (ArcIt a(_graph); a != INVALID; ++a, ++i) { 942 _arc_id[a] = i; 943 _source[i] = _node_id[_graph.source(a)]; 944 _target[i] = _node_id[_graph.target(a)]; 945 } 946 } 947 948 // Reset parameters 949 resetParams(); 915 950 return *this; 916 951 } … … 1025 1060 Value c = _lower[i]; 1026 1061 if (c >= 0) { 1027 _cap[i] = _upper[i] < INF? _upper[i]  c : INF;1062 _cap[i] = _upper[i] < MAX ? _upper[i]  c : INF; 1028 1063 } else { 1029 _cap[i] = _upper[i] < INF+ c ? _upper[i]  c : INF;1064 _cap[i] = _upper[i] < MAX + c ? _upper[i]  c : INF; 1030 1065 } 1031 1066 _supply[_source[i]] = c; … … 1055 1090 _state[i] = STATE_LOWER; 1056 1091 } 1057 1092 1058 1093 // Set data for the artificial root node 1059 1094 _root = _node_num; … … 1219 1254 e = _pred[u]; 1220 1255 d = _forward[u] ? 1221 _flow[e] : (_cap[e] == INF? INF : _cap[e]  _flow[e]);1256 _flow[e] : (_cap[e] >= MAX ? INF : _cap[e]  _flow[e]); 1222 1257 if (d < delta) { 1223 1258 delta = d; … … 1229 1264 for (int u = second; u != join; u = _parent[u]) { 1230 1265 e = _pred[u]; 1231 d = _forward[u] ? 1232 (_cap[e] == INF? INF : _cap[e]  _flow[e]) : _flow[e];1266 d = _forward[u] ? 1267 (_cap[e] >= MAX ? INF : _cap[e]  _flow[e]) : _flow[e]; 1233 1268 if (d <= delta) { 1234 1269 delta = d; … … 1331 1366 1332 1367 // Update _rev_thread using the new _thread values 1333 for (int i = 0; i <int(_dirty_revs.size()); ++i) {1368 for (int i = 0; i != int(_dirty_revs.size()); ++i) { 1334 1369 u = _dirty_revs[i]; 1335 1370 _rev_thread[_thread[u]] = u; … … 1401 1436 _pi[u] += sigma; 1402 1437 } 1438 } 1439 1440 // Heuristic initial pivots 1441 bool initialPivots() { 1442 Value curr, total = 0; 1443 std::vector<Node> supply_nodes, demand_nodes; 1444 for (NodeIt u(_graph); u != INVALID; ++u) { 1445 curr = _supply[_node_id[u]]; 1446 if (curr > 0) { 1447 total += curr; 1448 supply_nodes.push_back(u); 1449 } 1450 else if (curr < 0) { 1451 demand_nodes.push_back(u); 1452 } 1453 } 1454 if (_sum_supply > 0) total = _sum_supply; 1455 if (total <= 0) return true; 1456 1457 IntVector arc_vector; 1458 if (_sum_supply >= 0) { 1459 if (supply_nodes.size() == 1 && demand_nodes.size() == 1) { 1460 // Perform a reverse graph search from the sink to the source 1461 typename GR::template NodeMap<bool> reached(_graph, false); 1462 Node s = supply_nodes[0], t = demand_nodes[0]; 1463 std::vector<Node> stack; 1464 reached[t] = true; 1465 stack.push_back(t); 1466 while (!stack.empty()) { 1467 Node u, v = stack.back(); 1468 stack.pop_back(); 1469 if (v == s) break; 1470 for (InArcIt a(_graph, v); a != INVALID; ++a) { 1471 if (reached[u = _graph.source(a)]) continue; 1472 int j = _arc_id[a]; 1473 if (_cap[j] >= total) { 1474 arc_vector.push_back(j); 1475 reached[u] = true; 1476 stack.push_back(u); 1477 } 1478 } 1479 } 1480 } else { 1481 // Find the min. cost incomming arc for each demand node 1482 for (int i = 0; i != int(demand_nodes.size()); ++i) { 1483 Node v = demand_nodes[i]; 1484 Cost c, min_cost = std::numeric_limits<Cost>::max(); 1485 Arc min_arc = INVALID; 1486 for (InArcIt a(_graph, v); a != INVALID; ++a) { 1487 c = _cost[_arc_id[a]]; 1488 if (c < min_cost) { 1489 min_cost = c; 1490 min_arc = a; 1491 } 1492 } 1493 if (min_arc != INVALID) { 1494 arc_vector.push_back(_arc_id[min_arc]); 1495 } 1496 } 1497 } 1498 } else { 1499 // Find the min. cost outgoing arc for each supply node 1500 for (int i = 0; i != int(supply_nodes.size()); ++i) { 1501 Node u = supply_nodes[i]; 1502 Cost c, min_cost = std::numeric_limits<Cost>::max(); 1503 Arc min_arc = INVALID; 1504 for (OutArcIt a(_graph, u); a != INVALID; ++a) { 1505 c = _cost[_arc_id[a]]; 1506 if (c < min_cost) { 1507 min_cost = c; 1508 min_arc = a; 1509 } 1510 } 1511 if (min_arc != INVALID) { 1512 arc_vector.push_back(_arc_id[min_arc]); 1513 } 1514 } 1515 } 1516 1517 // Perform heuristic initial pivots 1518 for (int i = 0; i != int(arc_vector.size()); ++i) { 1519 in_arc = arc_vector[i]; 1520 if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]]  1521 _pi[_target[in_arc]]) >= 0) continue; 1522 findJoinNode(); 1523 bool change = findLeavingArc(); 1524 if (delta >= MAX) return false; 1525 changeFlow(change); 1526 if (change) { 1527 updateTreeStructure(); 1528 updatePotential(); 1529 } 1530 } 1531 return true; 1403 1532 } 1404 1533 … … 1425 1554 PivotRuleImpl pivot(*this); 1426 1555 1556 // Perform heuristic initial pivots 1557 if (!initialPivots()) return UNBOUNDED; 1558 1427 1559 // Execute the Network Simplex algorithm 1428 1560 while (pivot.findEnteringArc()) { 1429 1561 findJoinNode(); 1430 1562 bool change = findLeavingArc(); 1431 if (delta >= INF) return UNBOUNDED;1563 if (delta >= MAX) return UNBOUNDED; 1432 1564 changeFlow(change); 1433 1565 if (change) { … … 1436 1568 } 1437 1569 } 1438 1570 1439 1571 // Check feasibility 1440 1572 for (int e = _search_arc_num; e != _all_arc_num; ++e) { … … 1453 1585 } 1454 1586 } 1455 1587 1456 1588 // Shift potentials to meet the requirements of the GEQ/LEQ type 1457 1589 // optimality conditions
Note: See TracChangeset
for help on using the changeset viewer.