# HG changeset patch
# User kpeter
# Date 1226593070 0
# Node ID d239741cfd44c0410b4560fc400294531d17010f
# Parent  84354c78b06851ec727f4d0d1fbabd828b48e5b1
Various improvements in NetworkSimplex.

- Faster variant of "Altering Candidate List" pivot rule using make_heap
  instead of partial_sort.
- Doc improvements.
- Removing unecessary inline keywords.

diff -r 84354c78b068 -r d239741cfd44 lemon/network_simplex.h
--- a/lemon/network_simplex.h	Thu Nov 13 15:29:04 2008 +0000
+++ b/lemon/network_simplex.h	Thu Nov 13 16:17:50 2008 +0000
@@ -26,6 +26,7 @@
 
 #include <vector>
 #include <limits>
+#include <algorithm>
 
 #include <lemon/graph_adaptor.h>
 #include <lemon/graph_utils.h>
@@ -137,7 +138,7 @@
       ///\e
       Cost operator[](const Edge &e) const {
         return _cost_map[e] + _pot_map[_gr.source(e)]
-                           - _pot_map[_gr.target(e)];
+                            - _pot_map[_gr.target(e)];
       }
 
     }; //class ReducedCostMap
@@ -168,7 +169,7 @@
         _ns(ns), _edges(edges), _next_edge(0) {}
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
+      bool findEnteringEdge() {
         Edge e;
         for (int i = _next_edge; i < int(_edges.size()); ++i) {
           e = _edges[i];
@@ -212,7 +213,7 @@
         _ns(ns), _edges(edges) {}
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
+      bool findEnteringEdge() {
         Cost min = 0;
         Edge e;
         for (int i = 0; i < int(_edges.size()); ++i) {
@@ -259,7 +260,7 @@
       }
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
+      bool findEnteringEdge() {
         Cost curr, min = 0;
         Edge e;
         int cnt = _block_size;
@@ -336,11 +337,11 @@
       }
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
+      bool findEnteringEdge() {
         Cost min, curr;
         if (_curr_length > 0 && _minor_count < _minor_limit) {
-          // Minor iteration: selecting the best eligible edge from
-          // the current candidate list
+          // Minor iteration: select the best eligible edge from the
+          // current candidate list
           ++_minor_count;
           Edge e;
           min = 0;
@@ -358,7 +359,7 @@
           if (min < 0) return true;
         }
 
-        // Major iteration: building a new candidate list
+        // Major iteration: build a new candidate list
         Edge e;
         min = 0;
         _curr_length = 0;
@@ -423,7 +424,7 @@
       public:
         SortFunc(const SCostMap &map) : _map(map) {}
         bool operator()(const Edge &e1, const Edge &e2) {
-	  return _map[e1] < _map[e2];
+          return _map[e1] > _map[e2];
         }
       };
 
@@ -437,10 +438,10 @@
         _next_edge(0), _sort_func(_cand_cost)
       {
         // The main parameters of the pivot rule
-        const double BLOCK_SIZE_FACTOR = 1.0;
+        const double BLOCK_SIZE_FACTOR = 1.5;
         const int MIN_BLOCK_SIZE = 10;
         const double HEAD_LENGTH_FACTOR = 0.1;
-        const int MIN_HEAD_LENGTH = 5;
+        const int MIN_HEAD_LENGTH = 3;
 
         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
                                 MIN_BLOCK_SIZE );
@@ -451,17 +452,17 @@
       }
 
       /// Find next entering edge
-      inline bool findEnteringEdge() {
-        // Checking the current candidate list
+      bool findEnteringEdge() {
+        // Check the current candidate list
         Edge e;
-        for (int idx = 0; idx < _curr_length; ++idx) {
-          e = _candidates[idx];
+        for (int ix = 0; ix < _curr_length; ++ix) {
+          e = _candidates[ix];
           if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
-            _candidates[idx--] = _candidates[--_curr_length];
+            _candidates[ix--] = _candidates[--_curr_length];
           }
         }
 
-        // Extending the list
+        // Extend the list
         int cnt = _block_size;
         int last_edge = 0;
         int limit = _head_length;
@@ -494,24 +495,15 @@
         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);
+        // Make heap of the candidate list (approximating a partial sort)
+        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
+                   _sort_func );
 
+        // Pop the first element of the heap
         _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;
-        }
-
+        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
+                  _sort_func );
+        _curr_length = std::min(_head_length, _curr_length - 1);
         return true;
       }
     }; //class AlteringListPivotRule
@@ -558,7 +550,7 @@
     BoolNodeMap _forward;
     // The state edge map
     IntEdgeMap _state;
-    // The root node of the starting spanning tree
+    // The artificial root node of the spanning tree structure
     Node _root;
 
     // The reduced cost map
@@ -575,15 +567,15 @@
     NodeRefMap _node_ref;
     EdgeRefMap _edge_ref;
 
-    // The entering edge of the current pivot iteration.
+    // The entering edge of the current pivot iteration
     Edge _in_edge;
 
-    // Temporary nodes used in the current pivot iteration.
+    // Temporary nodes used in the current pivot iteration
     Node join, u_in, v_in, u_out, v_out;
     Node right, first, second, last;
     Node stem, par_stem, new_stem;
     // The maximum augment amount along the found cycle in the current
-    // pivot iteration.
+    // pivot iteration
     Capacity delta;
 
   public :
@@ -837,8 +829,8 @@
     /// list in every iteration (\ref AlteringListPivotRule).
     ///
     /// 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.
+    /// pivot rule proved to be 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 = BLOCK_SEARCH_PIVOT) {
@@ -911,12 +903,11 @@
 
   private:
 
-    /// \brief Extend the underlying graph and initialize all the
-    /// node and edge maps.
+    // Extend the underlying graph and initialize all the node and edge maps
     bool init() {
       if (!_valid_supply) return false;
 
-      // Initializing result maps
+      // Initialize result maps
       if (!_flow_result) {
         _flow_result = new FlowMap(_graph_ref);
         _local_flow = true;
@@ -926,7 +917,7 @@
         _local_potential = true;
       }
 
-      // Initializing the edge vector and the edge maps
+      // Initialize the edge vector and the edge maps
       _edges.reserve(countEdges(_graph));
       for (EdgeIt e(_graph); e != INVALID; ++e) {
         _edges.push_back(e);
@@ -934,7 +925,7 @@
         _state[e] = STATE_LOWER;
       }
 
-      // Adding an artificial root node to the graph
+      // Add an artificial root node to the graph
       _root = _graph.addNode();
       _parent[_root] = INVALID;
       _pred_edge[_root] = INVALID;
@@ -942,8 +933,8 @@
       _supply[_root] = 0;
       _potential[_root] = 0;
 
-      // Adding artificial edges to the graph and initializing the node
-      // maps of the spanning tree data structure
+      // Add artificial edges to the graph and initialize the node maps
+      // of the spanning tree data structure
       Node last = _root;
       Edge e;
       Cost max_cost = std::numeric_limits<Cost>::max() / 4;
@@ -974,8 +965,8 @@
       return true;
     }
 
-    /// Find the join node.
-    inline Node findJoinNode() {
+    // Find the join node
+    void findJoinNode() {
       Node u = _graph.source(_in_edge);
       Node v = _graph.target(_in_edge);
       while (u != v) {
@@ -986,14 +977,13 @@
         else if (_depth[u] > _depth[v]) u = _parent[u];
         else v = _parent[v];
       }
-      return u;
+      join = u;
     }
 
-    /// \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
+    // Find the leaving edge of the cycle and returns true if the 
+    // leaving edge is not the same as the entering edge
+    bool findLeavingEdge() {
+      // Initialize first and second nodes according to the direction
       // of the cycle
       if (_state[_in_edge] == STATE_LOWER) {
         first  = _graph.source(_in_edge);
@@ -1007,8 +997,7 @@
       Capacity d;
       Edge e;
 
-      // Searching the cycle along the path form the first node to the
-      // root node
+      // Search the cycle along the path form the first node to the root
       for (Node u = first; u != join; u = _parent[u]) {
         e = _pred_edge[u];
         d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
@@ -1020,8 +1009,7 @@
           result = true;
         }
       }
-      // Searching the cycle along the path form the second node to the
-      // root node
+      // Search the cycle along the path form the second node to the root
       for (Node u = second; u != join; u = _parent[u]) {
         e = _pred_edge[u];
         d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
@@ -1036,9 +1024,9 @@
       return result;
     }
 
-    /// Change \c flow and \c state edge maps.
-    inline void changeFlows(bool change) {
-      // Augmenting along the cycle
+    // Change _flow and _state edge maps
+    void changeFlows(bool change) {
+      // Augment along the cycle
       if (delta > 0) {
         Capacity val = _state[_in_edge] * delta;
         _flow[_in_edge] += val;
@@ -1049,7 +1037,7 @@
           _flow[_pred_edge[u]] += _forward[u] ? val : -val;
         }
       }
-      // Updating the state of the entering and leaving edges
+      // Update the state of the entering and leaving edges
       if (change) {
         _state[_in_edge] = STATE_TREE;
         _state[_pred_edge[u_out]] =
@@ -1059,12 +1047,12 @@
       }
     }
 
-    /// Update \c thread and \c parent node maps.
-    inline void updateThreadParent() {
+    // Update _thread and _parent node maps
+    void updateThreadParent() {
       Node u;
       v_out = _parent[u_out];
 
-      // Handling the case when join and v_out coincide
+      // Handle 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]) ;
@@ -1075,8 +1063,8 @@
         }
       }
 
-      // Finding the last successor of u_in (u) and the node after it
-      // (right) according to the thread index
+      // Find 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) {
@@ -1085,25 +1073,24 @@
       }
       else last = _thread[v_in];
 
-      // Updating stem nodes
+      // Update stem nodes
       _thread[v_in] = stem = u_in;
       par_stem = v_in;
       while (stem != u_out) {
         _thread[u] = new_stem = _parent[stem];
 
-        // Finding the node just before the stem node (u) according to
+        // Find 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;
 
-        // Changing the parent node of stem and shifting stem and
-        // par_stem nodes
+        // Change the parent node of stem and shift stem and par_stem nodes
         _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
+        // Find 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];
       }
@@ -1118,8 +1105,8 @@
       }
     }
 
-    /// Update \c pred_edge and \c forward node maps.
-    inline void updatePredEdge() {
+    // Update _pred_edge and _forward node maps.
+    void updatePredEdge() {
       Node u = u_out, v;
       while (u != u_in) {
         v = _parent[u];
@@ -1131,8 +1118,8 @@
       _forward[u_in] = (u_in == _graph.source(_in_edge));
     }
 
-    /// Update \c depth and \c potential node maps.
-    inline void updateDepthPotential() {
+    // Update _depth and _potential node maps
+    void updateDepthPotential() {
       _depth[u_in] = _depth[v_in] + 1;
       Cost sigma = _forward[u_in] ?
         _potential[v_in] - _potential[u_in] - _cost[_pred_edge[u_in]] :
@@ -1145,9 +1132,9 @@
       }
     }
 
-    /// Execute the algorithm.
+    // Execute the algorithm
     bool start(PivotRuleEnum pivot_rule) {
-      // Selecting the pivot rule implementation
+      // Select the pivot rule implementation
       switch (pivot_rule) {
         case FIRST_ELIGIBLE_PIVOT:
           return start<FirstEligiblePivotRule>();
@@ -1167,9 +1154,9 @@
     bool start() {
       PivotRuleImplementation pivot(*this, _edges);
 
-      // Executing the network simplex algorithm
+      // Execute the network simplex algorithm
       while (pivot.findEnteringEdge()) {
-        join = findJoinNode();
+        findJoinNode();
         bool change = findLeavingEdge();
         changeFlows(change);
         if (change) {
@@ -1179,14 +1166,13 @@
         }
       }
 
-      // Checking if the flow amount equals zero on all the artificial
-      // edges
+      // Check 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
+      // Copy 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]];
@@ -1194,7 +1180,7 @@
         for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
           (*_flow_result)[e] = _flow[_edge_ref[e]];
       }
-      // Copying potential values to _potential_result
+      // Copy potential values to _potential_result
       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
         (*_potential_result)[n] = _potential[_node_ref[n]];