gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Improve the Altering List pivot rule for NetworkSimplex (#435) Much less candidate arcs are preserved from an iteration to the next one and partial_sort() is used instead of heap operations.
0 1 0
default
1 file changed with 23 insertions and 20 deletions:
↑ Collapse diff ↑
Ignore white space 64 line context
... ...
@@ -93,98 +93,101 @@
93 93
      /// bounded), and the algorithm has found optimal flow and node
94 94
      /// potentials (primal and dual solutions).
95 95
      OPTIMAL,
96 96
      /// The objective function of the problem is unbounded, i.e.
97 97
      /// there is a directed cycle having negative total cost and
98 98
      /// infinite upper bound.
99 99
      UNBOUNDED
100 100
    };
101 101

	
102 102
    /// \brief Constants for selecting the type of the supply constraints.
103 103
    ///
104 104
    /// Enum type containing constants for selecting the supply type,
105 105
    /// i.e. the direction of the inequalities in the supply/demand
106 106
    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
107 107
    ///
108 108
    /// The default supply type is \c GEQ, the \c LEQ type can be
109 109
    /// selected using \ref supplyType().
110 110
    /// The equality form is a special case of both supply types.
111 111
    enum SupplyType {
112 112
      /// This option means that there are <em>"greater or equal"</em>
113 113
      /// supply/demand constraints in the definition of the problem.
114 114
      GEQ,
115 115
      /// This option means that there are <em>"less or equal"</em>
116 116
      /// supply/demand constraints in the definition of the problem.
117 117
      LEQ
118 118
    };
119 119

	
120 120
    /// \brief Constants for selecting the pivot rule.
121 121
    ///
122 122
    /// Enum type containing constants for selecting the pivot rule for
123 123
    /// the \ref run() function.
124 124
    ///
125
    /// \ref NetworkSimplex provides five different pivot rule
126
    /// implementations that significantly affect the running time
125
    /// \ref NetworkSimplex provides five different implementations for
126
    /// the pivot strategy that significantly affects the running time
127 127
    /// of the algorithm.
128
    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
129
    /// turend out to be the most efficient and the most robust on various
130
    /// test inputs.
131
    /// However, another pivot rule can be selected using the \ref run()
132
    /// function with the proper parameter.
128
    /// According to experimental tests conducted on various problem
129
    /// instances, \ref BLOCK_SEARCH "Block Search" and
130
    /// \ref ALTERING_LIST "Altering Candidate List" rules turned out
131
    /// to be the most efficient.
132
    /// Since \ref BLOCK_SEARCH "Block Search" is a simpler strategy that
133
    /// seemed to be slightly more robust, it is used by default.
134
    /// However, another pivot rule can easily be selected using the
135
    /// \ref run() function with the proper parameter.
133 136
    enum PivotRule {
134 137

	
135 138
      /// The \e First \e Eligible pivot rule.
136 139
      /// The next eligible arc is selected in a wraparound fashion
137 140
      /// in every iteration.
138 141
      FIRST_ELIGIBLE,
139 142

	
140 143
      /// The \e Best \e Eligible pivot rule.
141 144
      /// The best eligible arc is selected in every iteration.
142 145
      BEST_ELIGIBLE,
143 146

	
144 147
      /// The \e Block \e Search pivot rule.
145 148
      /// A specified number of arcs are examined in every iteration
146 149
      /// in a wraparound fashion and the best eligible arc is selected
147 150
      /// from this block.
148 151
      BLOCK_SEARCH,
149 152

	
150 153
      /// The \e Candidate \e List pivot rule.
151 154
      /// In a major iteration a candidate list is built from eligible arcs
152 155
      /// in a wraparound fashion and in the following minor iterations
153 156
      /// the best eligible arc is selected from this list.
154 157
      CANDIDATE_LIST,
155 158

	
156 159
      /// The \e Altering \e Candidate \e List pivot rule.
157 160
      /// It is a modified version of the Candidate List method.
158
      /// It keeps only the several best eligible arcs from the former
161
      /// It keeps only a few of the best eligible arcs from the former
159 162
      /// candidate list and extends this list in every iteration.
160 163
      ALTERING_LIST
161 164
    };
162 165

	
163 166
  private:
164 167

	
165 168
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
166 169

	
167 170
    typedef std::vector<int> IntVector;
168 171
    typedef std::vector<Value> ValueVector;
169 172
    typedef std::vector<Cost> CostVector;
170 173
    typedef std::vector<signed char> CharVector;
171 174
    // Note: vector<signed char> is used instead of vector<ArcState> and
172 175
    // vector<ArcDirection> for efficiency reasons
173 176

	
174 177
    // State constants for arcs
175 178
    enum ArcState {
176 179
      STATE_UPPER = -1,
177 180
      STATE_TREE  =  0,
178 181
      STATE_LOWER =  1
179 182
    };
180 183

	
181 184
    // Direction constants for tree arcs
182 185
    enum ArcDirection {
183 186
      DIR_DOWN = -1,
184 187
      DIR_UP   =  1
185 188
    };
186 189

	
187 190
  private:
188 191

	
189 192
    // Data related to the underlying digraph
190 193
    const GR &_graph;
... ...
@@ -509,151 +512,151 @@
509 512

	
510 513
    }; //class CandidateListPivotRule
511 514

	
512 515

	
513 516
    // Implementation of the Altering Candidate List pivot rule
514 517
    class AlteringListPivotRule
515 518
    {
516 519
    private:
517 520

	
518 521
      // References to the NetworkSimplex class
519 522
      const IntVector  &_source;
520 523
      const IntVector  &_target;
521 524
      const CostVector &_cost;
522 525
      const CharVector &_state;
523 526
      const CostVector &_pi;
524 527
      int &_in_arc;
525 528
      int _search_arc_num;
526 529

	
527 530
      // Pivot rule data
528 531
      int _block_size, _head_length, _curr_length;
529 532
      int _next_arc;
530 533
      IntVector _candidates;
531 534
      CostVector _cand_cost;
532 535

	
533 536
      // Functor class to compare arcs during sort of the candidate list
534 537
      class SortFunc
535 538
      {
536 539
      private:
537 540
        const CostVector &_map;
538 541
      public:
539 542
        SortFunc(const CostVector &map) : _map(map) {}
540 543
        bool operator()(int left, int right) {
541
          return _map[left] > _map[right];
544
          return _map[left] < _map[right];
542 545
        }
543 546
      };
544 547

	
545 548
      SortFunc _sort_func;
546 549

	
547 550
    public:
548 551

	
549 552
      // Constructor
550 553
      AlteringListPivotRule(NetworkSimplex &ns) :
551 554
        _source(ns._source), _target(ns._target),
552 555
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
553 556
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
554 557
        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
555 558
      {
556 559
        // The main parameters of the pivot rule
557 560
        const double BLOCK_SIZE_FACTOR = 1.0;
558 561
        const int MIN_BLOCK_SIZE = 10;
559
        const double HEAD_LENGTH_FACTOR = 0.1;
562
        const double HEAD_LENGTH_FACTOR = 0.01;
560 563
        const int MIN_HEAD_LENGTH = 3;
561 564

	
562 565
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
563 566
                                    std::sqrt(double(_search_arc_num))),
564 567
                                MIN_BLOCK_SIZE );
565 568
        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
566 569
                                 MIN_HEAD_LENGTH );
567 570
        _candidates.resize(_head_length + _block_size);
568 571
        _curr_length = 0;
569 572
      }
570 573

	
571 574
      // Find next entering arc
572 575
      bool findEnteringArc() {
573 576
        // Check the current candidate list
574 577
        int e;
575 578
        Cost c;
576 579
        for (int i = 0; i != _curr_length; ++i) {
577 580
          e = _candidates[i];
578 581
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
579 582
          if (c < 0) {
580 583
            _cand_cost[e] = c;
581 584
          } else {
582 585
            _candidates[i--] = _candidates[--_curr_length];
583 586
          }
584 587
        }
585 588

	
586 589
        // Extend the list
587 590
        int cnt = _block_size;
588 591
        int limit = _head_length;
589 592

	
590 593
        for (e = _next_arc; e != _search_arc_num; ++e) {
591 594
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
592 595
          if (c < 0) {
593 596
            _cand_cost[e] = c;
594 597
            _candidates[_curr_length++] = e;
595 598
          }
596 599
          if (--cnt == 0) {
597 600
            if (_curr_length > limit) goto search_end;
598 601
            limit = 0;
599 602
            cnt = _block_size;
600 603
          }
601 604
        }
602 605
        for (e = 0; e != _next_arc; ++e) {
603
          _cand_cost[e] = _state[e] *
604
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
605
          if (_cand_cost[e] < 0) {
606
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
607
          if (c < 0) {
608
            _cand_cost[e] = c;
606 609
            _candidates[_curr_length++] = e;
607 610
          }
608 611
          if (--cnt == 0) {
609 612
            if (_curr_length > limit) goto search_end;
610 613
            limit = 0;
611 614
            cnt = _block_size;
612 615
          }
613 616
        }
614 617
        if (_curr_length == 0) return false;
615 618

	
616 619
      search_end:
617 620

	
618
        // Make heap of the candidate list (approximating a partial sort)
619
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
620
                   _sort_func );
621
        // Perform partial sort operation on the candidate list
622
        int new_length = std::min(_head_length + 1, _curr_length);
623
        std::partial_sort(_candidates.begin(), _candidates.begin() + new_length,
624
                          _candidates.begin() + _curr_length, _sort_func);
621 625

	
622
        // Pop the first element of the heap
626
        // Select the entering arc and remove it from the list
623 627
        _in_arc = _candidates[0];
624 628
        _next_arc = e;
625
        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
626
                  _sort_func );
627
        _curr_length = std::min(_head_length, _curr_length - 1);
629
        _candidates[0] = _candidates[new_length - 1];
630
        _curr_length = new_length - 1;
628 631
        return true;
629 632
      }
630 633

	
631 634
    }; //class AlteringListPivotRule
632 635

	
633 636
  public:
634 637

	
635 638
    /// \brief Constructor.
636 639
    ///
637 640
    /// The constructor of the class.
638 641
    ///
639 642
    /// \param graph The digraph the algorithm runs on.
640 643
    /// \param arc_mixing Indicate if the arcs will be stored in a
641 644
    /// mixed order in the internal data structure.
642 645
    /// In general, it leads to similar performance as using the original
643 646
    /// arc order, but it makes the algorithm more robust and in special
644 647
    /// cases, even significantly faster. Therefore, it is enabled by default.
645 648
    NetworkSimplex(const GR& graph, bool arc_mixing = true) :
646 649
      _graph(graph), _node_id(graph), _arc_id(graph),
647 650
      _arc_mixing(arc_mixing),
648 651
      MAX(std::numeric_limits<Value>::max()),
649 652
      INF(std::numeric_limits<Value>::has_infinity ?
650 653
          std::numeric_limits<Value>::infinity() : MAX)
651 654
    {
652 655
      // Check the number types
653 656
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
654 657
        "The flow type of NetworkSimplex must be signed");
655 658
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
656 659
        "The cost type of NetworkSimplex must be signed");
657 660

	
658 661
      // Reset data structures
659 662
      reset();
0 comments (0 inline)