# HG changeset patch
# User kpeter
# Date 1223213803 0
# Node ID 30fb4d68b0e8d4c6b552e91ee18b2160a77aabf3
# Parent  6aa6fcaeaea51036b337104e3c3aa85b2c7b4573
Improve network simplex algorithm

 - Remove "Limited Search" and "Combined" pivot rules.
 - Add a new pivot rule "Altering Candidate List".
 - Make the edge selection faster in every pivot rule.
 - Set the default rule to "Block Search".
 - Doc improvements.

The algorithm became about 15-35 percent faster on various input files.
"Block Search" pivot rule proved to be by far the fastest on all inputs.

diff -r 6aa6fcaeaea5 -r 30fb4d68b0e8 lemon/network_simplex.h
--- a/lemon/network_simplex.h	Fri Sep 19 15:14:41 2008 +0000
+++ b/lemon/network_simplex.h	Sun Oct 05 13:36:43 2008 +0000
@@ -32,16 +32,18 @@
 #include <lemon/smart_graph.h>
 #include <lemon/math.h>
 
+#define _DEBUG_
+
 namespace lemon {
 
   /// \addtogroup min_cost_flow
   /// @{
 
-  /// \brief Implementation of the network simplex algorithm for
-  /// finding a minimum cost flow.
+  /// \brief Implementation of the primal network simplex algorithm
+  /// for finding a minimum cost flow.
   ///
-  /// \ref NetworkSimplex implements the network simplex algorithm for
-  /// finding a minimum cost flow.
+  /// \ref NetworkSimplex implements the primal network simplex algorithm
+  /// for finding a minimum cost flow.
   ///
   /// \tparam Graph The directed graph type the algorithm runs on.
   /// \tparam LowerMap The type of the lower bound map.
@@ -55,16 +57,15 @@
   /// - The value types of the maps should be convertible to each other.
   /// - \c CostMap::Value must be signed type.
   ///
-  /// \note \ref NetworkSimplex provides six different pivot rule
+  /// \note \ref NetworkSimplex provides five different pivot rule
   /// implementations that significantly affect the efficiency of the
   /// algorithm.
-  /// By default a combined pivot rule is used, which is the fastest
-  /// implementation according to our benchmark tests.
-  /// Another pivot rule can be selected using \ref run() function
-  /// with the proper parameter.
+  /// By default "Block Search" pivot rule is used, which proved to be
+  /// by far the most efficient according to our benchmark tests.
+  /// However another pivot rule can be selected using \ref run()
+  /// function with the proper parameter.
   ///
   /// \author Peter Kovacs
-
   template < typename Graph,
              typename LowerMap = typename Graph::template EdgeMap<int>,
              typename CapacityMap = typename Graph::template EdgeMap<int>,
@@ -89,10 +90,13 @@
     typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
     typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
     typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
+    typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;
 
     typedef typename Graph::template NodeMap<Node> NodeRefMap;
     typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
 
+    typedef std::vector<Edge> EdgeVector;
+
   public:
 
     /// The type of the flow map.
@@ -107,9 +111,8 @@
       FIRST_ELIGIBLE_PIVOT,
       BEST_ELIGIBLE_PIVOT,
       BLOCK_SEARCH_PIVOT,
-      LIMITED_SEARCH_PIVOT,
       CANDIDATE_LIST_PIVOT,
-      COMBINED_PIVOT
+      ALTERING_LIST_PIVOT
     };
 
   private:
@@ -148,32 +151,40 @@
     ///
     /// This class implements the "First Eligible" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
     class FirstEligiblePivotRule
     {
     private:
 
+      // References to the NetworkSimplex class
       NetworkSimplex &_ns;
-      EdgeIt _next_edge;
+      EdgeVector &_edges;
+
+      int _next_edge;
 
     public:
 
-      /// Constructor.
-      FirstEligiblePivotRule(NetworkSimplex &ns) :
-        _ns(ns), _next_edge(ns._graph) {}
+      /// Constructor
+      FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges), _next_edge(0) {}
 
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
-        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
+        Edge e;
+        for (int i = _next_edge; i < int(_edges.size()); ++i) {
+          e = _edges[i];
           if (_ns._state[e] * _ns._red_cost[e] < 0) {
             _ns._in_edge = e;
-            _next_edge = ++e;
+            _next_edge = i + 1;
             return true;
           }
         }
-        for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
+        for (int i = 0; i < _next_edge; ++i) {
+          e = _edges[i];
           if (_ns._state[e] * _ns._red_cost[e] < 0) {
             _ns._in_edge = e;
-            _next_edge = ++e;
+            _next_edge = i + 1;
             return true;
           }
         }
@@ -186,21 +197,28 @@
     ///
     /// This class implements the "Best Eligible" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
     class BestEligiblePivotRule
     {
     private:
 
+      // References to the NetworkSimplex class
       NetworkSimplex &_ns;
+      EdgeVector &_edges;
 
     public:
 
-      /// Constructor.
-      BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
+      /// Constructor
+      BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges) {}
 
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
         Cost min = 0;
-        for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
+        Edge e;
+        for (int i = 0; i < int(_edges.size()); ++i) {
+          e = _edges[i];
           if (_ns._state[e] * _ns._red_cost[e] < min) {
             min = _ns._state[e] * _ns._red_cost[e];
             _ns._in_edge = e;
@@ -215,206 +233,291 @@
     ///
     /// This class implements the "Block Search" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
     class BlockSearchPivotRule
     {
     private:
 
+      // References to the NetworkSimplex class
       NetworkSimplex &_ns;
-      EdgeIt _next_edge, _min_edge;
+      EdgeVector &_edges;
+
       int _block_size;
-
-      static const int MIN_BLOCK_SIZE = 10;
+      int _next_edge, _min_edge;
 
     public:
 
-      /// Constructor.
-      BlockSearchPivotRule(NetworkSimplex &ns) :
-        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
+      /// Constructor
+      BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
       {
-        _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
-        if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_SIZE;
+        // The main parameters of the pivot rule
+        const double BLOCK_SIZE_FACTOR = 2.0;
+        const int MIN_BLOCK_SIZE = 10;
+
+        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
+                                MIN_BLOCK_SIZE );
       }
 
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
         Cost curr, min = 0;
-        int cnt = 0;
-        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
+        Edge e;
+        int cnt = _block_size;
+        int i;
+        for (i = _next_edge; i < int(_edges.size()); ++i) {
+          e = _edges[i];
           if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
             min = curr;
-            _min_edge = e;
+            _min_edge = i;
           }
-          if (++cnt == _block_size) {
+          if (--cnt == 0) {
             if (min < 0) break;
-            cnt = 0;
+            cnt = _block_size;
           }
         }
-        if (min == 0) {
-          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
+        if (min == 0 || cnt > 0) {
+          for (i = 0; i < _next_edge; ++i) {
+            e = _edges[i];
             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
               min = curr;
-              _min_edge = e;
+              _min_edge = i;
             }
-            if (++cnt == _block_size) {
+            if (--cnt == 0) {
               if (min < 0) break;
-              cnt = 0;
+              cnt = _block_size;
             }
           }
         }
-        _ns._in_edge = _min_edge;
-        _next_edge = ++_min_edge;
-        return min < 0;
+        if (min >= 0) return false;
+        _ns._in_edge = _edges[_min_edge];
+        _next_edge = i;
+        return true;
       }
     }; //class BlockSearchPivotRule
 
-    /// \brief Implementation of the "Limited Search" pivot rule for the
-    /// \ref NetworkSimplex "network simplex" algorithm.
-    ///
-    /// This class implements the "Limited Search" pivot rule
-    /// for the \ref NetworkSimplex "network simplex" algorithm.
-    class LimitedSearchPivotRule
-    {
-    private:
-
-      NetworkSimplex &_ns;
-      EdgeIt _next_edge, _min_edge;
-      int _sample_size;
-
-      static const int SAMPLE_SIZE_FACTOR = 15;
-      static const int MIN_SAMPLE_SIZE = 10;
-
-    public:
-
-      /// Constructor.
-      LimitedSearchPivotRule(NetworkSimplex &ns) :
-        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
-      {
-        _sample_size = countEdges(_ns._graph) *
-                       SAMPLE_SIZE_FACTOR / 10000;
-        if (_sample_size < MIN_SAMPLE_SIZE)
-          _sample_size = MIN_SAMPLE_SIZE;
-      }
-
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
-        Cost curr, min = 0;
-        int cnt = 0;
-        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
-          if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
-            min = curr;
-            _min_edge = e;
-          }
-          if (curr < 0 && ++cnt == _sample_size) break;
-        }
-        if (min == 0) {
-          for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
-            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
-              min = curr;
-              _min_edge = e;
-            }
-            if (curr < 0 && ++cnt == _sample_size) break;
-          }
-        }
-        _ns._in_edge = _min_edge;
-        _next_edge = ++_min_edge;
-        return min < 0;
-      }
-    }; //class LimitedSearchPivotRule
-
     /// \brief Implementation of the "Candidate List" pivot rule for the
     /// \ref NetworkSimplex "network simplex" algorithm.
     ///
     /// This class implements the "Candidate List" pivot rule
     /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
     class CandidateListPivotRule
     {
     private:
 
+      // References to the NetworkSimplex class
       NetworkSimplex &_ns;
+      EdgeVector &_edges;
 
-      // The list of candidate edges.
-      std::vector<Edge> _candidates;
-      // The maximum length of the edge list.
-      int _list_length;
-      // The maximum number of minor iterations between two major
-      // itarations.
-      int _minor_limit;
-
-      int _minor_count;
-      EdgeIt _next_edge;
-
-      static const int LIST_LENGTH_FACTOR = 20;
-      static const int MINOR_LIMIT_FACTOR = 10;
-      static const int MIN_LIST_LENGTH = 10;
-      static const int MIN_MINOR_LIMIT = 2;
+      EdgeVector _candidates;
+      int _list_length, _minor_limit;
+      int _curr_length, _minor_count;
+      int _next_edge, _min_edge;
 
     public:
 
-      /// Constructor.
-      CandidateListPivotRule(NetworkSimplex &ns) :
-        _ns(ns), _next_edge(ns._graph)
+      /// Constructor
+      CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
       {
-        int edge_num = countEdges(_ns._graph);
-        _minor_count = 0;
-        _list_length = edge_num * LIST_LENGTH_FACTOR / 10000;
-        if (_list_length < MIN_LIST_LENGTH)
-          _list_length = MIN_LIST_LENGTH;
-        _minor_limit = _list_length * MINOR_LIMIT_FACTOR / 100;
-        if (_minor_limit < MIN_MINOR_LIMIT)
-          _minor_limit = MIN_MINOR_LIMIT;
+        // The main parameters of the pivot rule
+        const double LIST_LENGTH_FACTOR = 1.0;
+        const int MIN_LIST_LENGTH = 10;
+        const double MINOR_LIMIT_FACTOR = 0.1;
+        const int MIN_MINOR_LIMIT = 3;
+
+        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())),
+                                 MIN_LIST_LENGTH );
+        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
+                                 MIN_MINOR_LIMIT );
+        _curr_length = _minor_count = 0;
+        _candidates.resize(_list_length);
       }
 
-      /// Finds the next entering edge.
-      bool findEnteringEdge() {
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
         Cost min, curr;
-        if (_minor_count < _minor_limit && _candidates.size() > 0) {
-          // Minor iteration
+        if (_curr_length > 0 && _minor_count < _minor_limit) {
+          // Minor iteration: selecting the best eligible edge from
+          // the current candidate list
           ++_minor_count;
           Edge e;
           min = 0;
-          for (int i = 0; i < int(_candidates.size()); ++i) {
+          for (int i = 0; i < _curr_length; ++i) {
             e = _candidates[i];
-            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+            curr = _ns._state[e] * _ns._red_cost[e];
+            if (curr < min) {
               min = curr;
               _ns._in_edge = e;
             }
+            if (curr >= 0) {
+              _candidates[i--] = _candidates[--_curr_length];
+            }
           }
           if (min < 0) return true;
         }
 
-        // Major iteration
-        _candidates.clear();
-        EdgeIt e = _next_edge;
+        // Major iteration: building a new candidate list
+        Edge e;
         min = 0;
-        for ( ; e != INVALID; ++e) {
+        _curr_length = 0;
+        int i;
+        for (i = _next_edge; i < int(_edges.size()); ++i) {
+          e = _edges[i];
           if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
-            _candidates.push_back(e);
+            _candidates[_curr_length++] = e;
             if (curr < min) {
               min = curr;
-              _ns._in_edge = e;
+              _min_edge = i;
             }
-            if (int(_candidates.size()) == _list_length) break;
+            if (_curr_length == _list_length) break;
           }
         }
-        if (int(_candidates.size()) < _list_length) {
-          for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
+        if (_curr_length < _list_length) {
+          for (i = 0; i < _next_edge; ++i) {
+            e = _edges[i];
             if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
-              _candidates.push_back(e);
+              _candidates[_curr_length++] = e;
               if (curr < min) {
                 min = curr;
-                _ns._in_edge = e;
+                _min_edge = i;
               }
-              if (int(_candidates.size()) == _list_length) break;
+              if (_curr_length == _list_length) break;
             }
           }
         }
-        if (_candidates.size() == 0) return false;
+        if (_curr_length == 0) return false;
         _minor_count = 1;
-        _next_edge = ++e;
+        _ns._in_edge = _edges[_min_edge];
+        _next_edge = i;
         return true;
       }
     }; //class CandidateListPivotRule
 
+    /// \brief Implementation of the "Altering Candidate List" pivot rule
+    /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// This class implements the "Altering Candidate List" pivot rule
+    /// for the \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// For more information see \ref NetworkSimplex::run().
+    class AlteringListPivotRule
+    {
+    private:
+
+      // References to the NetworkSimplex class
+      NetworkSimplex &_ns;
+      EdgeVector &_edges;
+
+      EdgeVector _candidates;
+      SCostMap _cand_cost;
+      int _block_size, _head_length, _curr_length;
+      int _next_edge;
+
+      // Functor class to compare edges during sort of the candidate list
+      class SortFunc
+      {
+      private:
+        const SCostMap &_map;
+      public:
+        SortFunc(const SCostMap &map) : _map(map) {}
+        bool operator()(const Edge &e1, const Edge &e2) {
+	  return _map[e1] < _map[e2];
+        }
+      };
+
+      SortFunc _sort_func;
+
+    public:
+
+      /// Constructor
+      AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
+        _ns(ns), _edges(edges), _cand_cost(_ns._graph),
+        _next_edge(0), _sort_func(_cand_cost)
+      {
+        // The main parameters of the pivot rule
+        const double BLOCK_SIZE_FACTOR = 1.0;
+        const int MIN_BLOCK_SIZE = 10;
+        const double HEAD_LENGTH_FACTOR = 0.1;
+        const int MIN_HEAD_LENGTH = 5;
+
+        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
+                                MIN_BLOCK_SIZE );
+        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
+                                 MIN_HEAD_LENGTH );
+        _candidates.resize(_head_length + _block_size);
+        _curr_length = 0;
+      }
+
+      /// Find next entering edge
+      inline bool findEnteringEdge() {
+        // Checking the current candidate list
+        Edge e;
+        for (int idx = 0; idx < _curr_length; ++idx) {
+          e = _candidates[idx];
+          if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
+            _candidates[idx--] = _candidates[--_curr_length];
+          }
+        }
+
+        // Extending the list
+        int cnt = _block_size;
+        int last_edge = 0;
+        int limit = _head_length;
+        for (int i = _next_edge; i < int(_edges.size()); ++i) {
+          e = _edges[i];
+          if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+            _candidates[_curr_length++] = e;
+            last_edge = i;
+          }
+          if (--cnt == 0) {
+            if (_curr_length > limit) break;
+            limit = 0;
+            cnt = _block_size;
+          }
+        }
+        if (_curr_length <= limit) {
+          for (int i = 0; i < _next_edge; ++i) {
+            e = _edges[i];
+            if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
+              _candidates[_curr_length++] = e;
+              last_edge = i;
+            }
+            if (--cnt == 0) {
+              if (_curr_length > limit) break;
+              limit = 0;
+              cnt = _block_size;
+            }
+          }
+        }
+        if (_curr_length == 0) return false;
+        _next_edge = last_edge + 1;
+
+        // Sorting the list partially
+        EdgeVector::iterator sort_end = _candidates.begin();
+        EdgeVector::iterator vector_end = _candidates.begin();
+        for (int idx = 0; idx < _curr_length; ++idx) {
+          ++vector_end;
+          if (idx <= _head_length) ++sort_end;
+        }
+        partial_sort(_candidates.begin(), sort_end, vector_end, _sort_func);
+
+        _ns._in_edge = _candidates[0];
+        if (_curr_length > _head_length) {
+          _candidates[0] = _candidates[_head_length - 1];
+          _curr_length = _head_length - 1;
+        } else {
+          _candidates[0] = _candidates[_curr_length - 1];
+          --_curr_length;
+        }
+
+        return true;
+      }
+    }; //class AlteringListPivotRule
+
   private:
 
     // State constants for edges
@@ -424,9 +527,6 @@
       STATE_LOWER =  1
     };
 
-    // Constant for the combined pivot rule.
-    static const int COMBINED_PIVOT_MAX_DEG = 5;
-
   private:
 
     // The directed graph the algorithm runs on
@@ -466,6 +566,9 @@
     // The reduced cost map
     ReducedCostMap _red_cost;
 
+    // The non-artifical edges
+    EdgeVector _edges;
+
     // Members for handling the original graph
     FlowMap *_flow_result;
     PotentialMap *_potential_result;
@@ -672,9 +775,9 @@
       if (_local_potential) delete _potential_result;
     }
 
-    /// \brief Sets the flow map.
+    /// \brief Set the flow map.
     ///
-    /// Sets the flow map.
+    /// Set the flow map.
     ///
     /// \return \c (*this)
     NetworkSimplex& flowMap(FlowMap &map) {
@@ -686,9 +789,9 @@
       return *this;
     }
 
-    /// \brief Sets the potential map.
+    /// \brief Set the potential map.
     ///
-    /// Sets the potential map.
+    /// Set the potential map.
     ///
     /// \return \c (*this)
     NetworkSimplex& potentialMap(PotentialMap &map) {
@@ -701,8 +804,6 @@
     }
 
     /// \name Execution control
-    /// The only way to execute the algorithm is to call the run()
-    /// function.
 
     /// @{
 
@@ -726,50 +827,49 @@
     /// every iteration in a wraparound fashion and the best eligible
     /// edge is selected from this block (\ref BlockSearchPivotRule).
     ///
-    /// - LIMITED_SEARCH_PIVOT A specified number of eligible edges are
-    /// examined in every iteration in a wraparound fashion and the best
-    /// one is selected from them (\ref LimitedSearchPivotRule).
+    /// - CANDIDATE_LIST_PIVOT In a major iteration a candidate list is
+    /// built from eligible edges in a wraparound fashion and in the
+    /// following minor iterations the best eligible edge is selected
+    /// from this list (\ref CandidateListPivotRule).
     ///
-    /// - CANDIDATE_LIST_PIVOT In major iterations a candidate list is
-    /// built from eligible edges and it is used for edge selection in
-    /// the following minor iterations (\ref CandidateListPivotRule).
+    /// - ALTERING_LIST_PIVOT It is a modified version of the
+    /// "Candidate List" pivot rule. It keeps only the several best
+    /// eligible edges from the former candidate list and extends this
+    /// list in every iteration (\ref AlteringListPivotRule).
     ///
-    /// - COMBINED_PIVOT This is a combined version of the two fastest
-    /// pivot rules.
-    /// For rather sparse graphs \ref LimitedSearchPivotRule
-    /// "Limited Search" implementation is used, otherwise
-    /// \ref BlockSearchPivotRule "Block Search" pivot rule is used.
-    /// According to our benchmark tests this combined method is the
-    /// most efficient.
+    /// According to our comprehensive benchmark tests the "Block Search"
+    /// pivot rule proved to be by far the fastest and the most robust
+    /// on various test inputs. Thus it is the default option.
     ///
     /// \return \c true if a feasible flow can be found.
-    bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
+    bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
       return init() && start(pivot_rule);
     }
 
     /// @}
 
     /// \name Query Functions
-    /// The result of the algorithm can be obtained using these
-    /// functions.
-    /// \n run() must be called before using them.
+    /// The results of the algorithm can be obtained using these
+    /// functions.\n
+    /// \ref lemon::NetworkSimplex::run() "run()" must be called before
+    /// using them.
 
     /// @{
 
-    /// \brief Returns a const reference to the edge map storing the
+    /// \brief Return a const reference to the edge map storing the
     /// found flow.
     ///
-    /// Returns a const reference to the edge map storing the found flow.
+    /// Return a const reference to the edge map storing the found flow.
     ///
     /// \pre \ref run() must be called before using this function.
     const FlowMap& flowMap() const {
       return *_flow_result;
     }
 
-    /// \brief Returns a const reference to the node map storing the
+    /// \brief Return a const reference to the node map storing the
     /// found potentials (the dual solution).
     ///
-    /// Returns a const reference to the node map storing the found
+    /// Return a const reference to the node map storing the found
     /// potentials (the dual solution).
     ///
     /// \pre \ref run() must be called before using this function.
@@ -777,27 +877,27 @@
       return *_potential_result;
     }
 
-    /// \brief Returns the flow on the given edge.
+    /// \brief Return the flow on the given edge.
     ///
-    /// Returns the flow on the given edge.
+    /// Return the flow on the given edge.
     ///
     /// \pre \ref run() must be called before using this function.
     Capacity flow(const typename Graph::Edge& edge) const {
       return (*_flow_result)[edge];
     }
 
-    /// \brief Returns the potential of the given node.
+    /// \brief Return the potential of the given node.
     ///
-    /// Returns the potential of the given node.
+    /// Return the potential of the given node.
     ///
     /// \pre \ref run() must be called before using this function.
     Cost potential(const typename Graph::Node& node) const {
       return (*_potential_result)[node];
     }
 
-    /// \brief Returns the total cost of the found flow.
+    /// \brief Return the total cost of the found flow.
     ///
-    /// Returns the total cost of the found flow. The complexity of the
+    /// Return the total cost of the found flow. The complexity of the
     /// function is \f$ O(e) \f$.
     ///
     /// \pre \ref run() must be called before using this function.
@@ -812,7 +912,7 @@
 
   private:
 
-    /// \brief Extends the underlaying graph and initializes all the
+    /// \brief Extend the underlying graph and initialize all the
     /// node and edge maps.
     bool init() {
       if (!_valid_supply) return false;
@@ -827,8 +927,10 @@
         _local_potential = true;
       }
 
-      // Initializing state and flow maps
+      // Initializing the edge vector and the edge maps
+      _edges.reserve(countEdges(_graph));
       for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _edges.push_back(e);
         _flow[e] = 0;
         _state[e] = STATE_LOWER;
       }
@@ -873,8 +975,8 @@
       return true;
     }
 
-    /// Finds the join node.
-    Node findJoinNode() {
+    /// Find the join node.
+    inline Node findJoinNode() {
       Node u = _graph.source(_in_edge);
       Node v = _graph.target(_in_edge);
       while (u != v) {
@@ -888,9 +990,10 @@
       return u;
     }
 
-    /// \brief Finds the leaving edge of the cycle. Returns \c true if
-    /// the leaving edge is not the same as the entering edge.
-    bool findLeavingEdge() {
+    /// \brief Find the leaving edge of the cycle.
+    /// \return \c true if the leaving edge is not the same as the
+    /// entering edge.
+    inline bool findLeavingEdge() {
       // Initializing first and second nodes according to the direction
       // of the cycle
       if (_state[_in_edge] == STATE_LOWER) {
@@ -934,8 +1037,8 @@
       return result;
     }
 
-    /// Changes \c flow and \c state edge maps.
-    void changeFlows(bool change) {
+    /// Change \c flow and \c state edge maps.
+    inline void changeFlows(bool change) {
       // Augmenting along the cycle
       if (delta > 0) {
         Capacity val = _state[_in_edge] * delta;
@@ -957,8 +1060,8 @@
       }
     }
 
-    /// Updates \c thread and \c parent node maps.
-    void updateThreadParent() {
+    /// Update \c thread and \c parent node maps.
+    inline void updateThreadParent() {
       Node u;
       v_out = _parent[u_out];
 
@@ -1016,8 +1119,8 @@
       }
     }
 
-    /// Updates \c pred_edge and \c forward node maps.
-    void updatePredEdge() {
+    /// Update \c pred_edge and \c forward node maps.
+    inline void updatePredEdge() {
       Node u = u_out, v;
       while (u != u_in) {
         v = _parent[u];
@@ -1029,8 +1132,8 @@
       _forward[u_in] = (u_in == _graph.source(_in_edge));
     }
 
-    /// Updates \c depth and \c potential node maps.
-    void updateDepthPotential() {
+    /// Update \c depth and \c potential node maps.
+    inline void updateDepthPotential() {
       _depth[u_in] = _depth[v_in] + 1;
       _potential[u_in] = _forward[u_in] ?
         _potential[v_in] - _cost[_pred_edge[u_in]] :
@@ -1049,8 +1152,9 @@
       }
     }
 
-    /// Executes the algorithm.
+    /// Execute the algorithm.
     bool start(PivotRuleEnum pivot_rule) {
+      // Selecting the pivot rule implementation
       switch (pivot_rule) {
         case FIRST_ELIGIBLE_PIVOT:
           return start<FirstEligiblePivotRule>();
@@ -1058,23 +1162,17 @@
           return start<BestEligiblePivotRule>();
         case BLOCK_SEARCH_PIVOT:
           return start<BlockSearchPivotRule>();
-        case LIMITED_SEARCH_PIVOT:
-          return start<LimitedSearchPivotRule>();
         case CANDIDATE_LIST_PIVOT:
           return start<CandidateListPivotRule>();
-        case COMBINED_PIVOT:
-          if ( countEdges(_graph) / countNodes(_graph) <=
-               COMBINED_PIVOT_MAX_DEG )
-            return start<LimitedSearchPivotRule>();
-          else
-            return start<BlockSearchPivotRule>();
+        case ALTERING_LIST_PIVOT:
+          return start<AlteringListPivotRule>();
       }
       return false;
     }
 
     template<class PivotRuleImplementation>
     bool start() {
-      PivotRuleImplementation pivot(*this);
+      PivotRuleImplementation pivot(*this, _edges);
 
       // Executing the network simplex algorithm
       while (pivot.findEnteringEdge()) {