# HG changeset patch
# User kpeter
# Date 1203305526 0
# Node ID e866e288cba60565198685aed117ae069f86fcb3
# Parent  7058c9690e7d95ecf157d2fd577ca449ff0bd899
Major improvements in NetworkSimplex.

Main changes:
- Use -potenital[] instead of potential[] to conform to the usual
terminology.
- Use function parameter instead of #define commands to select pivot rule.
- Use much faster implementation for the candidate list pivot rule.
It is about 5-20 times faster now.
- Add a new pivot rule called "Limited Search" that is a modified
version of "Block Search". It is about 25 percent faster on rather
sparse graphs.
- By default "Limited Search" is used for sparse graphs and
"Block Search" is used otherwise. This combined method is the most
efficient on every input class.
- Change the name of private members to start with "_".
- Change the name of function parameters not to start with "_".
- Remove unnecessary documentation for private members.
- Many doc improvements.

diff -r 7058c9690e7d -r e866e288cba6 lemon/network_simplex.h
--- a/lemon/network_simplex.h	Mon Feb 18 03:30:53 2008 +0000
+++ b/lemon/network_simplex.h	Mon Feb 18 03:32:06 2008 +0000
@@ -22,47 +22,15 @@
 /// \ingroup min_cost_flow
 ///
 /// \file
-/// \brief The network simplex algorithm for finding a minimum cost flow.
+/// \brief Network simplex algorithm for finding a minimum cost flow.
 
+#include <vector>
 #include <limits>
+
 #include <lemon/graph_adaptor.h>
 #include <lemon/graph_utils.h>
 #include <lemon/smart_graph.h>
-
-/// \brief The pivot rule used in the algorithm.
-//#define FIRST_ELIGIBLE_PIVOT
-//#define BEST_ELIGIBLE_PIVOT
-#define EDGE_BLOCK_PIVOT
-//#define CANDIDATE_LIST_PIVOT
-//#define SORTED_LIST_PIVOT
-
-//#define _DEBUG_ITER_
-
-
-// State constant for edges at their lower bounds.
-#define LOWER   1
-// State constant for edges in the spanning tree.
-#define TREE    0
-// State constant for edges at their upper bounds.
-#define UPPER   -1
-
-#ifdef EDGE_BLOCK_PIVOT
-  #include <lemon/math.h>
-  #define MIN_BLOCK_SIZE        10
-#endif
-
-#ifdef CANDIDATE_LIST_PIVOT
-  #include <vector>
-  #define LIST_LENGTH_DIV       50
-  #define MINOR_LIMIT_DIV       200
-#endif
-
-#ifdef SORTED_LIST_PIVOT
-  #include <vector>
-  #include <algorithm>
-  #define LIST_LENGTH_DIV       100
-  #define LOWER_DIV             4
-#endif
+#include <lemon/math.h>
 
 namespace lemon {
 
@@ -75,31 +43,38 @@
   /// \ref NetworkSimplex implements the network simplex algorithm for
   /// finding a minimum cost flow.
   ///
-  /// \param Graph The directed graph type the algorithm runs on.
-  /// \param LowerMap The type of the lower bound map.
-  /// \param CapacityMap The type of the capacity (upper bound) map.
-  /// \param CostMap The type of the cost (length) map.
-  /// \param SupplyMap The type of the supply map.
+  /// \tparam Graph The directed graph type the algorithm runs on.
+  /// \tparam LowerMap The type of the lower bound map.
+  /// \tparam CapacityMap The type of the capacity (upper bound) map.
+  /// \tparam CostMap The type of the cost (length) map.
+  /// \tparam SupplyMap The type of the supply map.
   ///
   /// \warning
-  /// - Edge capacities and costs should be non-negative integers.
-  ///   However \c CostMap::Value should be signed type.
-  /// - Supply values should be signed integers.
-  /// - \c LowerMap::Value must be convertible to
-  ///   \c CapacityMap::Value and \c CapacityMap::Value must be
-  ///   convertible to \c SupplyMap::Value.
+  /// - Edge capacities and costs should be \e non-negative \e integers.
+  /// - Supply values should be \e signed \e integers.
+  /// - \c LowerMap::Value must be convertible to \c CapacityMap::Value.
+  /// - \c CapacityMap::Value and \c SupplyMap::Value must be
+  ///   convertible to each other.
+  /// - All value types must be convertible to \c CostMap::Value, which
+  ///   must be signed type.
+  ///
+  /// \note \ref NetworkSimplex provides six 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.
   ///
   /// \author Peter Kovacs
 
   template < typename Graph,
              typename LowerMap = typename Graph::template EdgeMap<int>,
-             typename CapacityMap = LowerMap,
+             typename CapacityMap = typename Graph::template EdgeMap<int>,
              typename CostMap = typename Graph::template EdgeMap<int>,
-             typename SupplyMap = typename Graph::template NodeMap
-                                  <typename CapacityMap::Value> >
+             typename SupplyMap = typename Graph::template NodeMap<int> >
   class NetworkSimplex
   {
-    typedef typename LowerMap::Value Lower;
     typedef typename CapacityMap::Value Capacity;
     typedef typename CostMap::Value Cost;
     typedef typename SupplyMap::Value Supply;
@@ -107,7 +82,6 @@
     typedef SmartGraph SGraph;
     GRAPH_TYPEDEFS(typename SGraph);
 
-    typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
     typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
     typedef typename SGraph::template EdgeMap<Cost> SCostMap;
     typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
@@ -129,77 +103,380 @@
     /// The type of the potential map.
     typedef typename Graph::template NodeMap<Cost> PotentialMap;
 
-  protected:
+  public:
 
+    /// Enum type to select the pivot rule used by \ref run().
+    enum PivotRuleEnum {
+      FIRST_ELIGIBLE_PIVOT,
+      BEST_ELIGIBLE_PIVOT,
+      BLOCK_SEARCH_PIVOT,
+      LIMITED_SEARCH_PIVOT,
+      CANDIDATE_LIST_PIVOT,
+      COMBINED_PIVOT
+    };
+
+  private:
+
+    /// \brief Map adaptor class for handling reduced edge costs.
+    ///
     /// Map adaptor class for handling reduced edge costs.
     class ReducedCostMap : public MapBase<Edge, Cost>
     {
     private:
 
-      const SGraph &gr;
-      const SCostMap &cost_map;
-      const SPotentialMap &pot_map;
+      const SGraph &_gr;
+      const SCostMap &_cost_map;
+      const SPotentialMap &_pot_map;
 
     public:
 
-      ReducedCostMap( const SGraph &_gr,
-                      const SCostMap &_cm,
-                      const SPotentialMap &_pm ) :
-        gr(_gr), cost_map(_cm), pot_map(_pm) {}
+      ///\e
+      ReducedCostMap( const SGraph &gr,
+                      const SCostMap &cost_map,
+                      const SPotentialMap &pot_map ) :
+        _gr(gr), _cost_map(cost_map), _pot_map(pm) {}
 
+      ///\e
       Cost operator[](const Edge &e) const {
-        return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
+        return _cost_map[e] + _pot_map[_gr.source(e)]
+                           - _pot_map[_gr.target(e)];
       }
 
     }; //class ReducedCostMap
 
-  protected:
+  private:
 
-    /// The directed graph the algorithm runs on.
-    SGraph graph;
-    /// The original graph.
-    const Graph &graph_ref;
-    /// The original lower bound map.
-    const LowerMap *lower;
-    /// The capacity map.
-    SCapacityMap capacity;
-    /// The cost map.
-    SCostMap cost;
-    /// The supply map.
-    SSupplyMap supply;
-    /// The reduced cost map.
-    ReducedCostMap red_cost;
-    bool valid_supply;
+    /// \brief Implementation of the "First Eligible" pivot rule for the
+    /// \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// This class implements the "First Eligible" pivot rule
+    /// for the \ref NetworkSimplex "network simplex" algorithm.
+    class FirstEligiblePivotRule
+    {
+    private:
 
-    /// The edge map of the current flow.
-    SCapacityMap flow;
-    /// The potential node map.
-    SPotentialMap potential;
-    FlowMap flow_result;
-    PotentialMap potential_result;
+      NetworkSimplex &_ns;
+      EdgeIt _next_edge;
 
-    /// Node reference for the original graph.
-    NodeRefMap node_ref;
-    /// Edge reference for the original graph.
-    EdgeRefMap edge_ref;
+    public:
 
-    /// The \c depth node map of the spanning tree structure.
-    IntNodeMap depth;
-    /// The \c parent node map of the spanning tree structure.
-    NodeNodeMap parent;
-    /// The \c pred_edge node map of the spanning tree structure.
-    EdgeNodeMap pred_edge;
-    /// The \c thread node map of the spanning tree structure.
-    NodeNodeMap thread;
-    /// The \c forward node map of the spanning tree structure.
-    BoolNodeMap forward;
-    /// The state edge map.
-    IntEdgeMap state;
-    /// The root node of the starting spanning tree.
-    Node root;
+      /// Constructor.
+      FirstEligiblePivotRule(NetworkSimplex &ns) :
+        _ns(ns), _next_edge(ns._graph) {}
+
+      /// Finds the next entering edge.
+      bool findEnteringEdge() {
+        for (EdgeIt e = _next_edge; e != INVALID; ++e) {
+          if (_ns._state[e] * _ns._red_cost[e] < 0) {
+            _ns._in_edge = e;
+            _next_edge = ++e;
+            return true;
+          }
+        }
+        for (EdgeIt e(_ns._graph); e != _next_edge; ++e) {
+          if (_ns._state[e] * _ns._red_cost[e] < 0) {
+            _ns._in_edge = e;
+            _next_edge = ++e;
+            return true;
+          }
+        }
+        return false;
+      }
+    }; //class FirstEligiblePivotRule
+
+    /// \brief Implementation of the "Best Eligible" pivot rule for the
+    /// \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// This class implements the "Best Eligible" pivot rule
+    /// for the \ref NetworkSimplex "network simplex" algorithm.
+    class BestEligiblePivotRule
+    {
+    private:
+
+      NetworkSimplex &_ns;
+
+    public:
+
+      /// Constructor.
+      BestEligiblePivotRule(NetworkSimplex &ns) : _ns(ns) {}
+
+      /// Finds the next entering edge.
+      bool findEnteringEdge() {
+        Cost min = 0;
+        for (EdgeIt e(_ns._graph); e != INVALID; ++e) {
+          if (_ns._state[e] * _ns._red_cost[e] < min) {
+            min = _ns._state[e] * _ns._red_cost[e];
+            _ns._in_edge = e;
+          }
+        }
+        return min < 0;
+      }
+    }; //class BestEligiblePivotRule
+
+    /// \brief Implementation of the "Block Search" pivot rule for the
+    /// \ref NetworkSimplex "network simplex" algorithm.
+    ///
+    /// This class implements the "Block Search" pivot rule
+    /// for the \ref NetworkSimplex "network simplex" algorithm.
+    class BlockSearchPivotRule
+    {
+    private:
+
+      NetworkSimplex &_ns;
+      EdgeIt _next_edge, _min_edge;
+      int _block_size;
+
+      static const int MIN_BLOCK_SIZE = 10;
+
+    public:
+
+      /// Constructor.
+      BlockSearchPivotRule(NetworkSimplex &ns) :
+        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
+      {
+        _block_size = 2 * int(sqrt(countEdges(_ns._graph)));
+        if (_block_size < MIN_BLOCK_SIZE) _block_size = MIN_BLOCK_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 (++cnt == _block_size) {
+            if (min < 0) break;
+            cnt = 0;
+          }
+        }
+        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 (++cnt == _block_size) {
+              if (min < 0) break;
+              cnt = 0;
+            }
+          }
+        }
+        _ns._in_edge = _min_edge;
+        _next_edge = ++_min_edge;
+        return min < 0;
+      }
+    }; //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 MIN_SAMPLE_SIZE = 10;
+      static const double SAMPLE_SIZE_FACTOR = 0.0015;
+
+    public:
+
+      /// Constructor.
+      LimitedSearchPivotRule(NetworkSimplex &ns) :
+        _ns(ns), _next_edge(ns._graph), _min_edge(ns._graph)
+      {
+        _sample_size = int(SAMPLE_SIZE_FACTOR * countEdges(_ns._graph));
+        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.
+    class CandidateListPivotRule
+    {
+    private:
+
+      NetworkSimplex &_ns;
+
+      // 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 double LIST_LENGTH_FACTOR = 0.002;
+      static const double MINOR_LIMIT_FACTOR = 0.1;
+      static const int MIN_LIST_LENGTH = 10;
+      static const int MIN_MINOR_LIMIT = 2;
+
+    public:
+
+      /// Constructor.
+      CandidateListPivotRule(NetworkSimplex &ns) :
+        _ns(ns), _next_edge(ns._graph)
+      {
+        int edge_num = countEdges(_ns._graph);
+        _minor_count = 0;
+        _list_length = int(edge_num * LIST_LENGTH_FACTOR);
+        if (_list_length < MIN_LIST_LENGTH)
+          _list_length = MIN_LIST_LENGTH;
+        _minor_limit = int(_list_length * MINOR_LIMIT_FACTOR);
+        if (_minor_limit < MIN_MINOR_LIMIT)
+          _minor_limit = MIN_MINOR_LIMIT;
+      }
+
+      /// Finds the next entering edge.
+      bool findEnteringEdge() {
+        Cost min, curr;
+        if (_minor_count < _minor_limit && _candidates.size() > 0) {
+          // Minor iteration
+          ++_minor_count;
+          Edge e;
+          min = 0;
+          for (int i = 0; i < int(_candidates.size()); ++i) {
+            e = _candidates[i];
+            if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
+              min = curr;
+              _ns._in_edge = e;
+            }
+          }
+          if (min < 0) return true;
+        }
+
+        // Major iteration
+        _candidates.clear();
+        EdgeIt e = _next_edge;
+        min = 0;
+        for ( ; e != INVALID; ++e) {
+          if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+            _candidates.push_back(e);
+            if (curr < min) {
+              min = curr;
+              _ns._in_edge = e;
+            }
+            if (int(_candidates.size()) == _list_length) break;
+          }
+        }
+        if (int(_candidates.size()) < _list_length) {
+          for (e = EdgeIt(_ns._graph); e != _next_edge; ++e) {
+            if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
+              _candidates.push_back(e);
+              if (curr < min) {
+                min = curr;
+                _ns._in_edge = e;
+              }
+              if (int(_candidates.size()) == _list_length) break;
+            }
+          }
+        }
+        if (_candidates.size() == 0) return false;
+        _minor_count = 1;
+        _next_edge = ++e;
+        return true;
+      }
+    }; //class CandidateListPivotRule
+
+  private:
+
+    // State constant for edges at their lower bounds.
+    static const int STATE_LOWER =  1;
+    // State constant for edges in the spanning tree.
+    static const int STATE_TREE  =  0;
+    // State constant for edges at their upper bounds.
+    static const int STATE_UPPER = -1;
+
+    // Constant for the combined pivot rule.
+    static const int COMBINED_PIVOT_MAX_DEG = 5;
+
+  private:
+
+    // The directed graph the algorithm runs on
+    SGraph _graph;
+    // The original graph
+    const Graph &_graph_ref;
+    // The original lower bound map
+    const LowerMap *_lower;
+    // The capacity map
+    SCapacityMap _capacity;
+    // The cost map
+    SCostMap _cost;
+    // The supply map
+    SSupplyMap _supply;
+    bool _valid_supply;
+
+    // Edge map of the current flow
+    SCapacityMap _flow;
+    // Node map of the current potentials
+    SPotentialMap _potential;
+
+    // The depth node map of the spanning tree structure
+    IntNodeMap _depth;
+    // The parent node map of the spanning tree structure
+    NodeNodeMap _parent;
+    // The pred_edge node map of the spanning tree structure
+    EdgeNodeMap _pred_edge;
+    // The thread node map of the spanning tree structure
+    NodeNodeMap _thread;
+    // The forward node map of the spanning tree structure
+    BoolNodeMap _forward;
+    // The state edge map
+    IntEdgeMap _state;
+    // The root node of the starting spanning tree
+    Node _root;
+
+    // The reduced cost map
+    ReducedCostMap _red_cost;
+
+    // Members for handling the original graph
+    FlowMap _flow_result;
+    PotentialMap _potential_result;
+    NodeRefMap _node_ref;
+    EdgeRefMap _edge_ref;
 
     // The entering edge of the current pivot iteration.
-    Edge in_edge;
+    Edge _in_edge;
+
     // Temporary nodes used in the current pivot iteration.
     Node join, u_in, v_in, u_out, v_out;
     Node right, first, second, last;
@@ -208,82 +485,56 @@
     // pivot iteration.
     Capacity delta;
 
-#ifdef EDGE_BLOCK_PIVOT
-    /// The size of blocks for the "Edge Block" pivot rule.
-    int block_size;
-#endif
-#ifdef CANDIDATE_LIST_PIVOT
-    /// \brief The list of candidate edges for the "Candidate List"
-    /// pivot rule.
-    std::vector<Edge> candidates;
-    /// \brief The maximum length of the edge list for the
-    /// "Candidate List" pivot rule.
-    int list_length;
-    /// \brief The maximum number of minor iterations between two major
-    /// itarations.
-    int minor_limit;
-    /// \brief The number of minor iterations.
-    int minor_count;
-#endif
-#ifdef SORTED_LIST_PIVOT
-    /// \brief The list of candidate edges for the "Sorted List"
-    /// pivot rule.
-    std::vector<Edge> candidates;
-    /// \brief The maximum length of the edge list for the
-    /// "Sorted List" pivot rule.
-    int list_length;
-    int list_index;
-#endif
-
   public :
 
     /// \brief General constructor of the class (with lower bounds).
     ///
     /// General constructor of the class (with lower bounds).
     ///
-    /// \param _graph The directed graph the algorithm runs on.
-    /// \param _lower The lower bounds of the edges.
-    /// \param _capacity The capacities (upper bounds) of the edges.
-    /// \param _cost The cost (length) values of the edges.
-    /// \param _supply The supply values of the nodes (signed).
-    NetworkSimplex( const Graph &_graph,
-                    const LowerMap &_lower,
-                    const CapacityMap &_capacity,
-                    const CostMap &_cost,
-                    const SupplyMap &_supply ) :
-      graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
-      supply(graph), flow(graph), flow_result(_graph), potential(graph),
-      potential_result(_graph), depth(graph), parent(graph),
-      pred_edge(graph), thread(graph), forward(graph), state(graph),
-      node_ref(graph_ref), edge_ref(graph_ref),
-      red_cost(graph, cost, potential)
+    /// \param graph The directed graph the algorithm runs on.
+    /// \param lower The lower bounds of the edges.
+    /// \param capacity The capacities (upper bounds) of the edges.
+    /// \param cost The cost (length) values of the edges.
+    /// \param supply The supply values of the nodes (signed).
+    NetworkSimplex( const Graph &graph,
+                    const LowerMap &lower,
+                    const CapacityMap &capacity,
+                    const CostMap &cost,
+                    const SupplyMap &supply ) :
+      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
+      _cost(_graph), _supply(_graph), _flow(_graph),
+      _potential(_graph), _depth(_graph), _parent(_graph),
+      _pred_edge(_graph), _thread(_graph), _forward(_graph),
+      _state(_graph), _red_cost(_graph, _cost, _potential),
+      _flow_result(graph), _potential_result(graph),
+      _node_ref(graph), _edge_ref(graph)
     {
       // Checking the sum of supply values
       Supply sum = 0;
-      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
-        sum += _supply[n];
-      if (!(valid_supply = sum == 0)) return;
+      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
+        sum += supply[n];
+      if (!(_valid_supply = sum == 0)) return;
 
-      // Copying graph_ref to graph
-      graph.reserveNode(countNodes(graph_ref) + 1);
-      graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));
-      copyGraph(graph, graph_ref)
-        .edgeMap(cost, _cost)
-        .nodeRef(node_ref)
-        .edgeRef(edge_ref)
+      // Copying _graph_ref to _graph
+      _graph.reserveNode(countNodes(_graph_ref) + 1);
+      _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
+      copyGraph(_graph, _graph_ref)
+        .edgeMap(_cost, cost)
+        .nodeRef(_node_ref)
+        .edgeRef(_edge_ref)
         .run();
 
       // Removing non-zero lower bounds
-      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
-        capacity[edge_ref[e]] = _capacity[e] - _lower[e];
+      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
+        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
       }
-      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
-        Supply s = _supply[n];
-        for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
-          s += _lower[e];
-        for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
-          s -= _lower[e];
-        supply[node_ref[n]] = s;
+      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
+        Supply s = supply[n];
+        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
+          s += lower[e];
+        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
+          s -= lower[e];
+        _supply[_node_ref[n]] = s;
       }
     }
 
@@ -291,34 +542,35 @@
     ///
     /// General constructor of the class (without lower bounds).
     ///
-    /// \param _graph The directed graph the algorithm runs on.
-    /// \param _capacity The capacities (upper bounds) of the edges.
-    /// \param _cost The cost (length) values of the edges.
-    /// \param _supply The supply values of the nodes (signed).
-    NetworkSimplex( const Graph &_graph,
-                    const CapacityMap &_capacity,
-                    const CostMap &_cost,
-                    const SupplyMap &_supply ) :
-      graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
-      supply(graph), flow(graph), flow_result(_graph), potential(graph),
-      potential_result(_graph), depth(graph), parent(graph),
-      pred_edge(graph), thread(graph), forward(graph), state(graph),
-      node_ref(graph_ref), edge_ref(graph_ref),
-      red_cost(graph, cost, potential)
+    /// \param graph The directed graph the algorithm runs on.
+    /// \param capacity The capacities (upper bounds) of the edges.
+    /// \param cost The cost (length) values of the edges.
+    /// \param supply The supply values of the nodes (signed).
+    NetworkSimplex( const Graph &graph,
+                    const CapacityMap &capacity,
+                    const CostMap &cost,
+                    const SupplyMap &supply ) :
+      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
+      _cost(_graph), _supply(_graph), _flow(_graph),
+      _potential(_graph), _depth(_graph), _parent(_graph),
+      _pred_edge(_graph), _thread(_graph), _forward(_graph),
+      _state(_graph), _red_cost(_graph, _cost, _potential),
+      _flow_result(graph), _potential_result(graph),
+      _node_ref(graph), _edge_ref(graph)
     {
       // Checking the sum of supply values
       Supply sum = 0;
-      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
-        sum += _supply[n];
-      if (!(valid_supply = sum == 0)) return;
+      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
+        sum += supply[n];
+      if (!(_valid_supply = sum == 0)) return;
 
-      // Copying graph_ref to graph
-      copyGraph(graph, graph_ref)
-        .edgeMap(capacity, _capacity)
-        .edgeMap(cost, _cost)
-        .nodeMap(supply, _supply)
-        .nodeRef(node_ref)
-        .edgeRef(edge_ref)
+      // Copying _graph_ref to graph
+      copyGraph(_graph, _graph_ref)
+        .edgeMap(_capacity, capacity)
+        .edgeMap(_cost, cost)
+        .nodeMap(_supply, supply)
+        .nodeRef(_node_ref)
+        .edgeRef(_edge_ref)
         .run();
     }
 
@@ -326,115 +578,150 @@
     ///
     /// Simple constructor of the class (with lower bounds).
     ///
-    /// \param _graph The directed graph the algorithm runs on.
-    /// \param _lower The lower bounds of the edges.
-    /// \param _capacity The capacities (upper bounds) of the edges.
-    /// \param _cost The cost (length) values of the edges.
-    /// \param _s The source node.
-    /// \param _t The target node.
-    /// \param _flow_value The required amount of flow from node \c _s
-    /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
-    NetworkSimplex( const Graph &_graph,
-                    const LowerMap &_lower,
-                    const CapacityMap &_capacity,
-                    const CostMap &_cost,
-                    typename Graph::Node _s,
-                    typename Graph::Node _t,
-                    typename SupplyMap::Value _flow_value ) :
-      graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
-      supply(graph), flow(graph), flow_result(_graph), potential(graph),
-      potential_result(_graph), depth(graph), parent(graph),
-      pred_edge(graph), thread(graph), forward(graph), state(graph),
-      node_ref(graph_ref), edge_ref(graph_ref),
-      red_cost(graph, cost, potential)
+    /// \param graph The directed graph the algorithm runs on.
+    /// \param lower The lower bounds of the edges.
+    /// \param capacity The capacities (upper bounds) of the edges.
+    /// \param cost The cost (length) values of the edges.
+    /// \param s The source node.
+    /// \param t The target node.
+    /// \param flow_value The required amount of flow from node \c s
+    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
+    NetworkSimplex( const Graph &graph,
+                    const LowerMap &lower,
+                    const CapacityMap &capacity,
+                    const CostMap &cost,
+                    typename Graph::Node s,
+                    typename Graph::Node t,
+                    typename SupplyMap::Value flow_value ) :
+      _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
+      _cost(_graph), _supply(_graph), _flow(_graph),
+      _potential(_graph), _depth(_graph), _parent(_graph),
+      _pred_edge(_graph), _thread(_graph), _forward(_graph),
+      _state(_graph), _red_cost(_graph, _cost, _potential),
+      _flow_result(graph), _potential_result(graph),
+      _node_ref(graph), _edge_ref(graph)
     {
-      // Copying graph_ref to graph
-      copyGraph(graph, graph_ref)
-        .edgeMap(cost, _cost)
-        .nodeRef(node_ref)
-        .edgeRef(edge_ref)
+      // Copying _graph_ref to graph
+      copyGraph(_graph, _graph_ref)
+        .edgeMap(_cost, cost)
+        .nodeRef(_node_ref)
+        .edgeRef(_edge_ref)
         .run();
 
       // Removing non-zero lower bounds
-      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
-        capacity[edge_ref[e]] = _capacity[e] - _lower[e];
+      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
+        _capacity[_edge_ref[e]] = capacity[e] - lower[e];
       }
-      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
-        Supply s = 0;
-        if (n == _s) s =  _flow_value;
-        if (n == _t) s = -_flow_value;
-        for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
-          s += _lower[e];
-        for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
-          s -= _lower[e];
-        supply[node_ref[n]] = s;
+      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n) {
+        Supply sum = 0;
+        if (n == s) sum =  flow_value;
+        if (n == t) sum = -flow_value;
+        for (typename Graph::InEdgeIt e(_graph_ref, n); e != INVALID; ++e)
+          sum += lower[e];
+        for (typename Graph::OutEdgeIt e(_graph_ref, n); e != INVALID; ++e)
+          sum -= lower[e];
+        _supply[_node_ref[n]] = sum;
       }
-      valid_supply = true;
+      _valid_supply = true;
     }
 
     /// \brief Simple constructor of the class (without lower bounds).
     ///
     /// Simple constructor of the class (without lower bounds).
     ///
-    /// \param _graph The directed graph the algorithm runs on.
-    /// \param _capacity The capacities (upper bounds) of the edges.
-    /// \param _cost The cost (length) values of the edges.
-    /// \param _s The source node.
-    /// \param _t The target node.
-    /// \param _flow_value The required amount of flow from node \c _s
-    /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
-    NetworkSimplex( const Graph &_graph,
-                    const CapacityMap &_capacity,
-                    const CostMap &_cost,
-                    typename Graph::Node _s,
-                    typename Graph::Node _t,
-                    typename SupplyMap::Value _flow_value ) :
-      graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
-      supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),
-      potential_result(_graph), depth(graph), parent(graph),
-      pred_edge(graph), thread(graph), forward(graph), state(graph),
-      node_ref(graph_ref), edge_ref(graph_ref),
-      red_cost(graph, cost, potential)
+    /// \param graph The directed graph the algorithm runs on.
+    /// \param capacity The capacities (upper bounds) of the edges.
+    /// \param cost The cost (length) values of the edges.
+    /// \param s The source node.
+    /// \param t The target node.
+    /// \param flow_value The required amount of flow from node \c s
+    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
+    NetworkSimplex( const Graph &graph,
+                    const CapacityMap &capacity,
+                    const CostMap &cost,
+                    typename Graph::Node s,
+                    typename Graph::Node t,
+                    typename SupplyMap::Value flow_value ) :
+      _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
+      _cost(_graph), _supply(_graph, 0), _flow(_graph),
+      _potential(_graph), _depth(_graph), _parent(_graph),
+      _pred_edge(_graph), _thread(_graph), _forward(_graph),
+      _state(_graph), _red_cost(_graph, _cost, _potential),
+      _flow_result(graph), _potential_result(graph),
+      _node_ref(graph), _edge_ref(graph)
     {
-      // Copying graph_ref to graph
-      copyGraph(graph, graph_ref)
-        .edgeMap(capacity, _capacity)
-        .edgeMap(cost, _cost)
-        .nodeRef(node_ref)
-        .edgeRef(edge_ref)
+      // Copying _graph_ref to graph
+      copyGraph(_graph, _graph_ref)
+        .edgeMap(_capacity, capacity)
+        .edgeMap(_cost, cost)
+        .nodeRef(_node_ref)
+        .edgeRef(_edge_ref)
         .run();
-      supply[node_ref[_s]] =  _flow_value;
-      supply[node_ref[_t]] = -_flow_value;
-      valid_supply = true;
+      _supply[_node_ref[s]] =  flow_value;
+      _supply[_node_ref[t]] = -flow_value;
+      _valid_supply = true;
     }
 
     /// \brief Runs the algorithm.
     ///
     /// Runs the algorithm.
     ///
+    /// \param pivot_rule The pivot rule that is used during the
+    /// algorithm.
+    ///
+    /// The available pivot rules:
+    ///
+    /// - FIRST_ELIGIBLE_PIVOT The next eligible edge is selected in
+    /// a wraparound fashion in every iteration
+    /// (\ref FirstEligiblePivotRule).
+    ///
+    /// - BEST_ELIGIBLE_PIVOT The best eligible edge is selected in
+    /// every iteration (\ref BestEligiblePivotRule).
+    ///
+    /// - BLOCK_SEARCH_PIVOT A specified number of edges are examined in
+    /// 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 major iterations a candidate list is
+    /// built from eligible edges and it is used for edge selection in
+    /// the following minor iterations (\ref CandidateListPivotRule).
+    ///
+    /// - 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.
+    ///
     /// \return \c true if a feasible flow can be found.
-    bool run() {
-      return init() && start();
+    bool run(PivotRuleEnum pivot_rule = COMBINED_PIVOT) {
+      return init() && start(pivot_rule);
     }
 
-    /// \brief Returns a const reference to the flow map.
+    /// \brief Returns a const reference to the edge map storing the
+    /// found flow.
     ///
-    /// Returns a const reference to the flow map.
+    /// Returns 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;
+      return _flow_result;
     }
 
-    /// \brief Returns a const reference to the potential map (the dual
-    /// solution).
+    /// \brief Returns a const reference to the node map storing the
+    /// found potentials (the dual solution).
     ///
-    /// Returns a const reference to the potential map (the dual
-    /// solution).
+    /// Returns a const reference to the node map storing the found
+    /// potentials (the dual solution).
     ///
     /// \pre \ref run() must be called before using this function.
     const PotentialMap& potentialMap() const {
-      return potential_result;
+      return _potential_result;
     }
 
     /// \brief Returns the total cost of the found flow.
@@ -445,248 +732,75 @@
     /// \pre \ref run() must be called before using this function.
     Cost totalCost() const {
       Cost c = 0;
-      for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
-        c += flow_result[e] * cost[edge_ref[e]];
+      for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
+        c += _flow_result[e] * _cost[_edge_ref[e]];
       return c;
     }
 
-  protected:
+  private:
 
     /// \brief Extends the underlaying graph and initializes all the
     /// node and edge maps.
     bool init() {
-      if (!valid_supply) return false;
+      if (!_valid_supply) return false;
 
       // Initializing state and flow maps
-      for (EdgeIt e(graph); e != INVALID; ++e) {
-        flow[e] = 0;
-        state[e] = LOWER;
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        _flow[e] = 0;
+        _state[e] = STATE_LOWER;
       }
 
       // Adding an artificial root node to the graph
-      root = graph.addNode();
-      parent[root] = INVALID;
-      pred_edge[root] = INVALID;
-      depth[root] = 0;
-      supply[root] = 0;
-      potential[root] = 0;
+      _root = _graph.addNode();
+      _parent[_root] = INVALID;
+      _pred_edge[_root] = INVALID;
+      _depth[_root] = 0;
+      _supply[_root] = 0;
+      _potential[_root] = 0;
 
       // Adding artificial edges to the graph and initializing the node
       // maps of the spanning tree data structure
-      Supply sum = 0;
-      Node last = root;
+      Node last = _root;
       Edge e;
       Cost max_cost = std::numeric_limits<Cost>::max() / 4;
-      for (NodeIt u(graph); u != INVALID; ++u) {
-        if (u == root) continue;
-        thread[last] = u;
+      for (NodeIt u(_graph); u != INVALID; ++u) {
+        if (u == _root) continue;
+        _thread[last] = u;
         last = u;
-        parent[u] = root;
-        depth[u] = 1;
-        sum += supply[u];
-        if (supply[u] >= 0) {
-          e = graph.addEdge(u, root);
-          flow[e] = supply[u];
-          forward[u] = true;
-          potential[u] = max_cost;
+        _parent[u] = _root;
+        _depth[u] = 1;
+        if (_supply[u] >= 0) {
+          e = _graph.addEdge(u, _root);
+          _flow[e] = _supply[u];
+          _forward[u] = true;
+          _potential[u] = -max_cost;
         } else {
-          e = graph.addEdge(root, u);
-          flow[e] = -supply[u];
-          forward[u] = false;
-          potential[u] = -max_cost;
+          e = _graph.addEdge(_root, u);
+          _flow[e] = -_supply[u];
+          _forward[u] = false;
+          _potential[u] = max_cost;
         }
-        cost[e] = max_cost;
-        capacity[e] = std::numeric_limits<Capacity>::max();
-        state[e] = TREE;
-        pred_edge[u] = e;
+        _cost[e] = max_cost;
+        _capacity[e] = std::numeric_limits<Capacity>::max();
+        _state[e] = STATE_TREE;
+        _pred_edge[u] = e;
       }
-      thread[last] = root;
+      _thread[last] = _root;
 
-#ifdef EDGE_BLOCK_PIVOT
-      // Initializing block_size for the edge block pivot rule
-      int edge_num = countEdges(graph);
-      block_size = 2 * int(sqrt(countEdges(graph)));
-      if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE;
-#endif
-#ifdef CANDIDATE_LIST_PIVOT
-      int edge_num = countEdges(graph);
-      minor_count = 0;
-      list_length = edge_num / LIST_LENGTH_DIV;
-      minor_limit = edge_num / MINOR_LIMIT_DIV;
-#endif
-#ifdef SORTED_LIST_PIVOT
-      int edge_num = countEdges(graph);
-      list_index = 0;
-      list_length = edge_num / LIST_LENGTH_DIV;
-#endif
-
-      return sum == 0;
+      return true;
     }
 
-#ifdef FIRST_ELIGIBLE_PIVOT
-    /// \brief Finds entering edge according to the "First Eligible"
-    /// pivot rule.
-    bool findEnteringEdge(EdgeIt &next_edge) {
-      for (EdgeIt e = next_edge; e != INVALID; ++e) {
-        if (state[e] * red_cost[e] < 0) {
-          in_edge = e;
-          next_edge = ++e;
-          return true;
+    /// Finds the join node.
+    Node findJoinNode() {
+      Node u = _graph.source(_in_edge);
+      Node v = _graph.target(_in_edge);
+      while (u != v) {
+        if (_depth[u] == _depth[v]) {
+          u = _parent[u];
+          v = _parent[v];
         }
-      }
-      for (EdgeIt e(graph); e != next_edge; ++e) {
-        if (state[e] * red_cost[e] < 0) {
-          in_edge = e;
-          next_edge = ++e;
-          return true;
-        }
-      }
-      return false;
-    }
-#endif
-
-#ifdef BEST_ELIGIBLE_PIVOT
-    /// \brief Finds entering edge according to the "Best Eligible"
-    /// pivot rule.
-    bool findEnteringEdge() {
-      Cost min = 0;
-      for (EdgeIt e(graph); e != INVALID; ++e) {
-        if (state[e] * red_cost[e] < min) {
-          min = state[e] * red_cost[e];
-          in_edge = e;
-        }
-      }
-      return min < 0;
-    }
-#endif
-
-#ifdef EDGE_BLOCK_PIVOT
-    /// \brief Finds entering edge according to the "Edge Block"
-    /// pivot rule.
-    bool findEnteringEdge(EdgeIt &next_edge) {
-      // Performing edge block selection
-      Cost curr, min = 0;
-      EdgeIt min_edge(graph);
-      int cnt = 0;
-      for (EdgeIt e = next_edge; e != INVALID; ++e) {
-        if ((curr = state[e] * red_cost[e]) < min) {
-          min = curr;
-          min_edge = e;
-        }
-        if (++cnt == block_size) {
-          if (min < 0) break;
-          cnt = 0;
-        }
-      }
-      if (!(min < 0)) {
-        for (EdgeIt e(graph); e != next_edge; ++e) {
-          if ((curr = state[e] * red_cost[e]) < min) {
-            min = curr;
-            min_edge = e;
-          }
-          if (++cnt == block_size) {
-            if (min < 0) break;
-            cnt = 0;
-          }
-        }
-      }
-      in_edge = min_edge;
-      if ((next_edge = ++min_edge) == INVALID)
-        next_edge = EdgeIt(graph);
-      return min < 0;
-    }
-#endif
-
-#ifdef CANDIDATE_LIST_PIVOT
-    /// \brief Finds entering edge according to the "Candidate List"
-    /// pivot rule.
-    bool findEnteringEdge() {
-      typedef typename std::vector<Edge>::iterator ListIt;
-
-      if (minor_count >= minor_limit || candidates.size() == 0) {
-        // Major iteration
-        candidates.clear();
-        for (EdgeIt e(graph); e != INVALID; ++e) {
-          if (state[e] * red_cost[e] < 0) {
-            candidates.push_back(e);
-            if (candidates.size() == list_length) break;
-          }
-        }
-        if (candidates.size() == 0) return false;
-      }
-
-      // Minor iteration
-      ++minor_count;
-      Cost min = 0;
-      Edge e;
-      for (int i = 0; i < candidates.size(); ++i) {
-        e = candidates[i];
-        if (state[e] * red_cost[e] < min) {
-          min = state[e] * red_cost[e];
-          in_edge = e;
-        }
-      }
-      return true;
-    }
-#endif
-
-#ifdef SORTED_LIST_PIVOT
-    /// \brief Functor class to compare edges during sort of the
-    /// candidate list.
-    class SortFunc
-    {
-    private:
-      const IntEdgeMap &st;
-      const ReducedCostMap &rc;
-    public:
-      SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
-        st(_st), rc(_rc) {}
-      bool operator()(const Edge &e1, const Edge &e2) {
-        return st[e1] * rc[e1] < st[e2] * rc[e2];
-      }
-    };
-
-    /// \brief Finds entering edge according to the "Sorted List"
-    /// pivot rule.
-    bool findEnteringEdge() {
-      static SortFunc sort_func(state, red_cost);
-
-      // Minor iteration
-      while (list_index < candidates.size()) {
-        in_edge = candidates[list_index++];
-        if (state[in_edge] * red_cost[in_edge] < 0) return true;
-      }
-
-      // Major iteration
-      candidates.clear();
-      Cost curr, min = 0;
-      for (EdgeIt e(graph); e != INVALID; ++e) {
-        if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
-          candidates.push_back(e);
-          if (curr < min) min = curr;
-          if (candidates.size() == list_length) break;
-        }
-      }
-      if (candidates.size() == 0) return false;
-      sort(candidates.begin(), candidates.end(), sort_func);
-      in_edge = candidates[0];
-      list_index = 1;
-      return true;
-    }
-#endif
-
-    /// \brief Finds the join node.
-    Node findJoinNode() {
-      // Finding the join node
-      Node u = graph.source(in_edge);
-      Node v = graph.target(in_edge);
-      while (u != v) {
-        if (depth[u] == depth[v]) {
-          u = parent[u];
-          v = parent[v];
-        }
-        else if (depth[u] > depth[v]) u = parent[u];
-        else v = parent[v];
+        else if (_depth[u] > _depth[v]) u = _parent[u];
+        else v = _parent[v];
       }
       return u;
     }
@@ -696,23 +810,23 @@
     bool findLeavingEdge() {
       // Initializing first and second nodes according to the direction
       // of the cycle
-      if (state[in_edge] == LOWER) {
-        first = graph.source(in_edge);
-        second  = graph.target(in_edge);
+      if (_state[_in_edge] == STATE_LOWER) {
+        first  = _graph.source(_in_edge);
+        second = _graph.target(_in_edge);
       } else {
-        first = graph.target(in_edge);
-        second  = graph.source(in_edge);
+        first  = _graph.target(_in_edge);
+        second = _graph.source(_in_edge);
       }
-      delta = capacity[in_edge];
+      delta = _capacity[_in_edge];
       bool result = false;
       Capacity d;
       Edge e;
 
       // Searching the cycle along the path form the first node to the
       // root node
-      for (Node u = first; u != join; u = parent[u]) {
-        e = pred_edge[u];
-        d = forward[u] ? flow[e] : capacity[e] - flow[e];
+      for (Node u = first; u != join; u = _parent[u]) {
+        e = _pred_edge[u];
+        d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
         if (d < delta) {
           delta = d;
           u_out = u;
@@ -723,9 +837,9 @@
       }
       // Searching the cycle along the path form the second node to the
       // root node
-      for (Node u = second; u != join; u = parent[u]) {
-        e = pred_edge[u];
-        d = forward[u] ? capacity[e] - flow[e] : flow[e];
+      for (Node u = second; u != join; u = _parent[u]) {
+        e = _pred_edge[u];
+        d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
         if (d <= delta) {
           delta = d;
           u_out = u;
@@ -737,134 +851,150 @@
       return result;
     }
 
-    /// \brief Changes \c flow and \c state edge maps.
+    /// Changes \c flow and \c state edge maps.
     void changeFlows(bool change) {
       // Augmenting along the cycle
       if (delta > 0) {
-        Capacity val = state[in_edge] * delta;
-        flow[in_edge] += val;
-        for (Node u = graph.source(in_edge); u != join; u = parent[u]) {
-          flow[pred_edge[u]] += forward[u] ? -val : val;
+        Capacity val = _state[_in_edge] * delta;
+        _flow[_in_edge] += val;
+        for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
+          _flow[_pred_edge[u]] += _forward[u] ? -val : val;
         }
-        for (Node u = graph.target(in_edge); u != join; u = parent[u]) {
-          flow[pred_edge[u]] += forward[u] ? val : -val;
+        for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
+          _flow[_pred_edge[u]] += _forward[u] ? val : -val;
         }
       }
       // Updating the state of the entering and leaving edges
       if (change) {
-        state[in_edge] = TREE;
-        state[pred_edge[u_out]] =
-          (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
+        _state[_in_edge] = STATE_TREE;
+        _state[_pred_edge[u_out]] =
+          (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
       } else {
-        state[in_edge] = -state[in_edge];
+        _state[_in_edge] = -_state[_in_edge];
       }
     }
 
-    /// \brief Updates \c thread and \c parent node maps.
+    /// Updates \c thread and \c parent node maps.
     void updateThreadParent() {
       Node u;
-      v_out = parent[u_out];
+      v_out = _parent[u_out];
 
       // Handling the case when join and v_out coincide
       bool par_first = false;
       if (join == v_out) {
-        for (u = join; u != u_in && u != v_in; u = thread[u]) ;
+        for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
         if (u == v_in) {
           par_first = true;
-          while (thread[u] != u_out) u = thread[u];
+          while (_thread[u] != u_out) u = _thread[u];
           first = u;
         }
       }
 
       // Finding the last successor of u_in (u) and the node after it
       // (right) according to the thread index
-      for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;
-      right = thread[u];
-      if (thread[v_in] == u_out) {
-        for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
-        if (last == u_out) last = thread[last];
+      for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
+      right = _thread[u];
+      if (_thread[v_in] == u_out) {
+        for (last = u; _depth[last] > _depth[u_out]; last = _thread[last]) ;
+        if (last == u_out) last = _thread[last];
       }
-      else last = thread[v_in];
+      else last = _thread[v_in];
 
       // Updating stem nodes
-      thread[v_in] = stem = u_in;
+      _thread[v_in] = stem = u_in;
       par_stem = v_in;
       while (stem != u_out) {
-        thread[u] = new_stem = parent[stem];
+        _thread[u] = new_stem = _parent[stem];
 
         // Finding the node just before the stem node (u) according to
         // the original thread index
-        for (u = new_stem; thread[u] != stem; u = thread[u]) ;
-        thread[u] = right;
+        for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
+        _thread[u] = right;
 
         // Changing the parent node of stem and shifting stem and
         // par_stem nodes
-        parent[stem] = par_stem;
+        _parent[stem] = par_stem;
         par_stem = stem;
         stem = new_stem;
 
         // Finding the last successor of stem (u) and the node after it
         // (right) according to the thread index
-        for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
-        right = thread[u];
+        for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
+        right = _thread[u];
       }
-      parent[u_out] = par_stem;
-      thread[u] = last;
+      _parent[u_out] = par_stem;
+      _thread[u] = last;
 
       if (join == v_out && par_first) {
-        if (first != v_in) thread[first] = right;
+        if (first != v_in) _thread[first] = right;
       } else {
-        for (u = v_out; thread[u] != u_out; u = thread[u]) ;
-        thread[u] = right;
+        for (u = v_out; _thread[u] != u_out; u = _thread[u]) ;
+        _thread[u] = right;
       }
     }
 
-    /// \brief Updates \c pred_edge and \c forward node maps.
+    /// Updates \c pred_edge and \c forward node maps.
     void updatePredEdge() {
       Node u = u_out, v;
       while (u != u_in) {
-        v = parent[u];
-        pred_edge[u] = pred_edge[v];
-        forward[u] = !forward[v];
+        v = _parent[u];
+        _pred_edge[u] = _pred_edge[v];
+        _forward[u] = !_forward[v];
         u = v;
       }
-      pred_edge[u_in] = in_edge;
-      forward[u_in] = (u_in == graph.source(in_edge));
+      _pred_edge[u_in] = _in_edge;
+      _forward[u_in] = (u_in == _graph.source(_in_edge));
     }
 
-    /// \brief Updates \c depth and \c potential node maps.
+    /// Updates \c depth and \c potential node maps.
     void updateDepthPotential() {
-      depth[u_in] = depth[v_in] + 1;
-      potential[u_in] = forward[u_in] ?
-        potential[v_in] + cost[pred_edge[u_in]] :
-        potential[v_in] - cost[pred_edge[u_in]];
+      _depth[u_in] = _depth[v_in] + 1;
+      _potential[u_in] = _forward[u_in] ?
+        _potential[v_in] - _cost[_pred_edge[u_in]] :
+        _potential[v_in] + _cost[_pred_edge[u_in]];
 
-      Node u = thread[u_in], v;
+      Node u = _thread[u_in], v;
       while (true) {
-        v = parent[u];
+        v = _parent[u];
         if (v == INVALID) break;
-        depth[u] = depth[v] + 1;
-        potential[u] = forward[u] ?
-          potential[v] + cost[pred_edge[u]] :
-          potential[v] - cost[pred_edge[u]];
-        if (depth[u] <= depth[v_in]) break;
-        u = thread[u];
+        _depth[u] = _depth[v] + 1;
+        _potential[u] = _forward[u] ?
+          _potential[v] - _cost[_pred_edge[u]] :
+          _potential[v] + _cost[_pred_edge[u]];
+        if (_depth[u] <= _depth[v_in]) break;
+        u = _thread[u];
       }
     }
 
-    /// \brief Executes the algorithm.
+    /// Executes the algorithm.
+    bool start(PivotRuleEnum pivot_rule) {
+      switch (pivot_rule) {
+        case FIRST_ELIGIBLE_PIVOT:
+          return start<FirstEligiblePivotRule>();
+        case BEST_ELIGIBLE_PIVOT:
+          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>();
+      }
+      return false;
+    }
+
+    template<class PivotRuleImplementation>
     bool start() {
-      // Processing pivots
-#ifdef _DEBUG_ITER_
-      int iter_num = 0;
-#endif
-#if defined(FIRST_ELIGIBLE_PIVOT) || defined(EDGE_BLOCK_PIVOT)
-      EdgeIt next_edge(graph);
-      while (findEnteringEdge(next_edge))
-#else
-      while (findEnteringEdge())
-#endif
-      {
+      PivotRuleImplementation pivot(*this);
+
+      // Executing the network simplex algorithm
+      while (pivot.findEnteringEdge()) {
         join = findJoinNode();
         bool change = findLeavingEdge();
         changeFlows(change);
@@ -873,34 +1003,26 @@
           updatePredEdge();
           updateDepthPotential();
         }
-#ifdef _DEBUG_ITER_
-        ++iter_num;
-#endif
       }
 
-#ifdef _DEBUG_ITER_
-      std::cout << "Network Simplex algorithm finished. " << iter_num
-                << " pivot iterations performed." << std::endl;
-#endif
+      // Checking if the flow amount equals zero on all the artificial
+      // edges
+      for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
+        if (_flow[e] > 0) return false;
+      for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
+        if (_flow[e] > 0) return false;
 
-      // Checking if the flow amount equals zero on all the
-      // artificial edges
-      for (InEdgeIt e(graph, root); e != INVALID; ++e)
-        if (flow[e] > 0) return false;
-      for (OutEdgeIt e(graph, root); e != INVALID; ++e)
-        if (flow[e] > 0) return false;
-
-      // Copying flow values to flow_result
-      if (lower) {
-        for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
-          flow_result[e] = (*lower)[e] + flow[edge_ref[e]];
+      // Copying flow values to _flow_result
+      if (_lower) {
+        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
+          _flow_result[e] = (*_lower)[e] + _flow[_edge_ref[e]];
       } else {
-        for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
-          flow_result[e] = flow[edge_ref[e]];
+        for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
+          _flow_result[e] = _flow[_edge_ref[e]];
       }
-      // Copying potential values to potential_result
-      for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
-        potential_result[n] = potential[node_ref[n]];
+      // Copying potential values to _potential_result
+      for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
+        _potential_result[n] = _potential[_node_ref[n]];
 
       return true;
     }