1.1 --- a/lemon/network_simplex.h	Thu Nov 13 15:29:04 2008 +0000
     1.2 +++ b/lemon/network_simplex.h	Thu Nov 13 16:17:50 2008 +0000
     1.3 @@ -26,6 +26,7 @@
     1.4  
     1.5  #include <vector>
     1.6  #include <limits>
     1.7 +#include <algorithm>
     1.8  
     1.9  #include <lemon/graph_adaptor.h>
    1.10  #include <lemon/graph_utils.h>
    1.11 @@ -137,7 +138,7 @@
    1.12        ///\e
    1.13        Cost operator[](const Edge &e) const {
    1.14          return _cost_map[e] + _pot_map[_gr.source(e)]
    1.15 -                           - _pot_map[_gr.target(e)];
    1.16 +                            - _pot_map[_gr.target(e)];
    1.17        }
    1.18  
    1.19      }; //class ReducedCostMap
    1.20 @@ -168,7 +169,7 @@
    1.21          _ns(ns), _edges(edges), _next_edge(0) {}
    1.22  
    1.23        /// Find next entering edge
    1.24 -      inline bool findEnteringEdge() {
    1.25 +      bool findEnteringEdge() {
    1.26          Edge e;
    1.27          for (int i = _next_edge; i < int(_edges.size()); ++i) {
    1.28            e = _edges[i];
    1.29 @@ -212,7 +213,7 @@
    1.30          _ns(ns), _edges(edges) {}
    1.31  
    1.32        /// Find next entering edge
    1.33 -      inline bool findEnteringEdge() {
    1.34 +      bool findEnteringEdge() {
    1.35          Cost min = 0;
    1.36          Edge e;
    1.37          for (int i = 0; i < int(_edges.size()); ++i) {
    1.38 @@ -259,7 +260,7 @@
    1.39        }
    1.40  
    1.41        /// Find next entering edge
    1.42 -      inline bool findEnteringEdge() {
    1.43 +      bool findEnteringEdge() {
    1.44          Cost curr, min = 0;
    1.45          Edge e;
    1.46          int cnt = _block_size;
    1.47 @@ -336,11 +337,11 @@
    1.48        }
    1.49  
    1.50        /// Find next entering edge
    1.51 -      inline bool findEnteringEdge() {
    1.52 +      bool findEnteringEdge() {
    1.53          Cost min, curr;
    1.54          if (_curr_length > 0 && _minor_count < _minor_limit) {
    1.55 -          // Minor iteration: selecting the best eligible edge from
    1.56 -          // the current candidate list
    1.57 +          // Minor iteration: select the best eligible edge from the
    1.58 +          // current candidate list
    1.59            ++_minor_count;
    1.60            Edge e;
    1.61            min = 0;
    1.62 @@ -358,7 +359,7 @@
    1.63            if (min < 0) return true;
    1.64          }
    1.65  
    1.66 -        // Major iteration: building a new candidate list
    1.67 +        // Major iteration: build a new candidate list
    1.68          Edge e;
    1.69          min = 0;
    1.70          _curr_length = 0;
    1.71 @@ -423,7 +424,7 @@
    1.72        public:
    1.73          SortFunc(const SCostMap &map) : _map(map) {}
    1.74          bool operator()(const Edge &e1, const Edge &e2) {
    1.75 -	  return _map[e1] < _map[e2];
    1.76 +          return _map[e1] > _map[e2];
    1.77          }
    1.78        };
    1.79  
    1.80 @@ -437,10 +438,10 @@
    1.81          _next_edge(0), _sort_func(_cand_cost)
    1.82        {
    1.83          // The main parameters of the pivot rule
    1.84 -        const double BLOCK_SIZE_FACTOR = 1.0;
    1.85 +        const double BLOCK_SIZE_FACTOR = 1.5;
    1.86          const int MIN_BLOCK_SIZE = 10;
    1.87          const double HEAD_LENGTH_FACTOR = 0.1;
    1.88 -        const int MIN_HEAD_LENGTH = 5;
    1.89 +        const int MIN_HEAD_LENGTH = 3;
    1.90  
    1.91          _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
    1.92                                  MIN_BLOCK_SIZE );
    1.93 @@ -451,17 +452,17 @@
    1.94        }
    1.95  
    1.96        /// Find next entering edge
    1.97 -      inline bool findEnteringEdge() {
    1.98 -        // Checking the current candidate list
    1.99 +      bool findEnteringEdge() {
   1.100 +        // Check the current candidate list
   1.101          Edge e;
   1.102 -        for (int idx = 0; idx < _curr_length; ++idx) {
   1.103 -          e = _candidates[idx];
   1.104 +        for (int ix = 0; ix < _curr_length; ++ix) {
   1.105 +          e = _candidates[ix];
   1.106            if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
   1.107 -            _candidates[idx--] = _candidates[--_curr_length];
   1.108 +            _candidates[ix--] = _candidates[--_curr_length];
   1.109            }
   1.110          }
   1.111  
   1.112 -        // Extending the list
   1.113 +        // Extend the list
   1.114          int cnt = _block_size;
   1.115          int last_edge = 0;
   1.116          int limit = _head_length;
   1.117 @@ -494,24 +495,15 @@
   1.118          if (_curr_length == 0) return false;
   1.119          _next_edge = last_edge + 1;
   1.120  
   1.121 -        // Sorting the list partially
   1.122 -        EdgeVector::iterator sort_end = _candidates.begin();
   1.123 -        EdgeVector::iterator vector_end = _candidates.begin();
   1.124 -        for (int idx = 0; idx < _curr_length; ++idx) {
   1.125 -          ++vector_end;
   1.126 -          if (idx <= _head_length) ++sort_end;
   1.127 -        }
   1.128 -        partial_sort(_candidates.begin(), sort_end, vector_end, _sort_func);
   1.129 +        // Make heap of the candidate list (approximating a partial sort)
   1.130 +        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   1.131 +                   _sort_func );
   1.132  
   1.133 +        // Pop the first element of the heap
   1.134          _ns._in_edge = _candidates[0];
   1.135 -        if (_curr_length > _head_length) {
   1.136 -          _candidates[0] = _candidates[_head_length - 1];
   1.137 -          _curr_length = _head_length - 1;
   1.138 -        } else {
   1.139 -          _candidates[0] = _candidates[_curr_length - 1];
   1.140 -          --_curr_length;
   1.141 -        }
   1.142 -
   1.143 +        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   1.144 +                  _sort_func );
   1.145 +        _curr_length = std::min(_head_length, _curr_length - 1);
   1.146          return true;
   1.147        }
   1.148      }; //class AlteringListPivotRule
   1.149 @@ -558,7 +550,7 @@
   1.150      BoolNodeMap _forward;
   1.151      // The state edge map
   1.152      IntEdgeMap _state;
   1.153 -    // The root node of the starting spanning tree
   1.154 +    // The artificial root node of the spanning tree structure
   1.155      Node _root;
   1.156  
   1.157      // The reduced cost map
   1.158 @@ -575,15 +567,15 @@
   1.159      NodeRefMap _node_ref;
   1.160      EdgeRefMap _edge_ref;
   1.161  
   1.162 -    // The entering edge of the current pivot iteration.
   1.163 +    // The entering edge of the current pivot iteration
   1.164      Edge _in_edge;
   1.165  
   1.166 -    // Temporary nodes used in the current pivot iteration.
   1.167 +    // Temporary nodes used in the current pivot iteration
   1.168      Node join, u_in, v_in, u_out, v_out;
   1.169      Node right, first, second, last;
   1.170      Node stem, par_stem, new_stem;
   1.171      // The maximum augment amount along the found cycle in the current
   1.172 -    // pivot iteration.
   1.173 +    // pivot iteration
   1.174      Capacity delta;
   1.175  
   1.176    public :
   1.177 @@ -837,8 +829,8 @@
   1.178      /// list in every iteration (\ref AlteringListPivotRule).
   1.179      ///
   1.180      /// According to our comprehensive benchmark tests the "Block Search"
   1.181 -    /// pivot rule proved to be by far the fastest and the most robust
   1.182 -    /// on various test inputs. Thus it is the default option.
   1.183 +    /// pivot rule proved to be the fastest and the most robust on
   1.184 +    /// various test inputs. Thus it is the default option.
   1.185      ///
   1.186      /// \return \c true if a feasible flow can be found.
   1.187      bool run(PivotRuleEnum pivot_rule = BLOCK_SEARCH_PIVOT) {
   1.188 @@ -911,12 +903,11 @@
   1.189  
   1.190    private:
   1.191  
   1.192 -    /// \brief Extend the underlying graph and initialize all the
   1.193 -    /// node and edge maps.
   1.194 +    // Extend the underlying graph and initialize all the node and edge maps
   1.195      bool init() {
   1.196        if (!_valid_supply) return false;
   1.197  
   1.198 -      // Initializing result maps
   1.199 +      // Initialize result maps
   1.200        if (!_flow_result) {
   1.201          _flow_result = new FlowMap(_graph_ref);
   1.202          _local_flow = true;
   1.203 @@ -926,7 +917,7 @@
   1.204          _local_potential = true;
   1.205        }
   1.206  
   1.207 -      // Initializing the edge vector and the edge maps
   1.208 +      // Initialize the edge vector and the edge maps
   1.209        _edges.reserve(countEdges(_graph));
   1.210        for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.211          _edges.push_back(e);
   1.212 @@ -934,7 +925,7 @@
   1.213          _state[e] = STATE_LOWER;
   1.214        }
   1.215  
   1.216 -      // Adding an artificial root node to the graph
   1.217 +      // Add an artificial root node to the graph
   1.218        _root = _graph.addNode();
   1.219        _parent[_root] = INVALID;
   1.220        _pred_edge[_root] = INVALID;
   1.221 @@ -942,8 +933,8 @@
   1.222        _supply[_root] = 0;
   1.223        _potential[_root] = 0;
   1.224  
   1.225 -      // Adding artificial edges to the graph and initializing the node
   1.226 -      // maps of the spanning tree data structure
   1.227 +      // Add artificial edges to the graph and initialize the node maps
   1.228 +      // of the spanning tree data structure
   1.229        Node last = _root;
   1.230        Edge e;
   1.231        Cost max_cost = std::numeric_limits<Cost>::max() / 4;
   1.232 @@ -974,8 +965,8 @@
   1.233        return true;
   1.234      }
   1.235  
   1.236 -    /// Find the join node.
   1.237 -    inline Node findJoinNode() {
   1.238 +    // Find the join node
   1.239 +    void findJoinNode() {
   1.240        Node u = _graph.source(_in_edge);
   1.241        Node v = _graph.target(_in_edge);
   1.242        while (u != v) {
   1.243 @@ -986,14 +977,13 @@
   1.244          else if (_depth[u] > _depth[v]) u = _parent[u];
   1.245          else v = _parent[v];
   1.246        }
   1.247 -      return u;
   1.248 +      join = u;
   1.249      }
   1.250  
   1.251 -    /// \brief Find the leaving edge of the cycle.
   1.252 -    /// \return \c true if the leaving edge is not the same as the
   1.253 -    /// entering edge.
   1.254 -    inline bool findLeavingEdge() {
   1.255 -      // Initializing first and second nodes according to the direction
   1.256 +    // Find the leaving edge of the cycle and returns true if the 
   1.257 +    // leaving edge is not the same as the entering edge
   1.258 +    bool findLeavingEdge() {
   1.259 +      // Initialize first and second nodes according to the direction
   1.260        // of the cycle
   1.261        if (_state[_in_edge] == STATE_LOWER) {
   1.262          first  = _graph.source(_in_edge);
   1.263 @@ -1007,8 +997,7 @@
   1.264        Capacity d;
   1.265        Edge e;
   1.266  
   1.267 -      // Searching the cycle along the path form the first node to the
   1.268 -      // root node
   1.269 +      // Search the cycle along the path form the first node to the root
   1.270        for (Node u = first; u != join; u = _parent[u]) {
   1.271          e = _pred_edge[u];
   1.272          d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
   1.273 @@ -1020,8 +1009,7 @@
   1.274            result = true;
   1.275          }
   1.276        }
   1.277 -      // Searching the cycle along the path form the second node to the
   1.278 -      // root node
   1.279 +      // Search the cycle along the path form the second node to the root
   1.280        for (Node u = second; u != join; u = _parent[u]) {
   1.281          e = _pred_edge[u];
   1.282          d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
   1.283 @@ -1036,9 +1024,9 @@
   1.284        return result;
   1.285      }
   1.286  
   1.287 -    /// Change \c flow and \c state edge maps.
   1.288 -    inline void changeFlows(bool change) {
   1.289 -      // Augmenting along the cycle
   1.290 +    // Change _flow and _state edge maps
   1.291 +    void changeFlows(bool change) {
   1.292 +      // Augment along the cycle
   1.293        if (delta > 0) {
   1.294          Capacity val = _state[_in_edge] * delta;
   1.295          _flow[_in_edge] += val;
   1.296 @@ -1049,7 +1037,7 @@
   1.297            _flow[_pred_edge[u]] += _forward[u] ? val : -val;
   1.298          }
   1.299        }
   1.300 -      // Updating the state of the entering and leaving edges
   1.301 +      // Update the state of the entering and leaving edges
   1.302        if (change) {
   1.303          _state[_in_edge] = STATE_TREE;
   1.304          _state[_pred_edge[u_out]] =
   1.305 @@ -1059,12 +1047,12 @@
   1.306        }
   1.307      }
   1.308  
   1.309 -    /// Update \c thread and \c parent node maps.
   1.310 -    inline void updateThreadParent() {
   1.311 +    // Update _thread and _parent node maps
   1.312 +    void updateThreadParent() {
   1.313        Node u;
   1.314        v_out = _parent[u_out];
   1.315  
   1.316 -      // Handling the case when join and v_out coincide
   1.317 +      // Handle the case when join and v_out coincide
   1.318        bool par_first = false;
   1.319        if (join == v_out) {
   1.320          for (u = join; u != u_in && u != v_in; u = _thread[u]) ;
   1.321 @@ -1075,8 +1063,8 @@
   1.322          }
   1.323        }
   1.324  
   1.325 -      // Finding the last successor of u_in (u) and the node after it
   1.326 -      // (right) according to the thread index
   1.327 +      // Find the last successor of u_in (u) and the node after it (right)
   1.328 +      // according to the thread index
   1.329        for (u = u_in; _depth[_thread[u]] > _depth[u_in]; u = _thread[u]) ;
   1.330        right = _thread[u];
   1.331        if (_thread[v_in] == u_out) {
   1.332 @@ -1085,25 +1073,24 @@
   1.333        }
   1.334        else last = _thread[v_in];
   1.335  
   1.336 -      // Updating stem nodes
   1.337 +      // Update stem nodes
   1.338        _thread[v_in] = stem = u_in;
   1.339        par_stem = v_in;
   1.340        while (stem != u_out) {
   1.341          _thread[u] = new_stem = _parent[stem];
   1.342  
   1.343 -        // Finding the node just before the stem node (u) according to
   1.344 +        // Find the node just before the stem node (u) according to
   1.345          // the original thread index
   1.346          for (u = new_stem; _thread[u] != stem; u = _thread[u]) ;
   1.347          _thread[u] = right;
   1.348  
   1.349 -        // Changing the parent node of stem and shifting stem and
   1.350 -        // par_stem nodes
   1.351 +        // Change the parent node of stem and shift stem and par_stem nodes
   1.352          _parent[stem] = par_stem;
   1.353          par_stem = stem;
   1.354          stem = new_stem;
   1.355  
   1.356 -        // Finding the last successor of stem (u) and the node after it
   1.357 -        // (right) according to the thread index
   1.358 +        // Find the last successor of stem (u) and the node after it (right)
   1.359 +        // according to the thread index
   1.360          for (u = stem; _depth[_thread[u]] > _depth[stem]; u = _thread[u]) ;
   1.361          right = _thread[u];
   1.362        }
   1.363 @@ -1118,8 +1105,8 @@
   1.364        }
   1.365      }
   1.366  
   1.367 -    /// Update \c pred_edge and \c forward node maps.
   1.368 -    inline void updatePredEdge() {
   1.369 +    // Update _pred_edge and _forward node maps.
   1.370 +    void updatePredEdge() {
   1.371        Node u = u_out, v;
   1.372        while (u != u_in) {
   1.373          v = _parent[u];
   1.374 @@ -1131,8 +1118,8 @@
   1.375        _forward[u_in] = (u_in == _graph.source(_in_edge));
   1.376      }
   1.377  
   1.378 -    /// Update \c depth and \c potential node maps.
   1.379 -    inline void updateDepthPotential() {
   1.380 +    // Update _depth and _potential node maps
   1.381 +    void updateDepthPotential() {
   1.382        _depth[u_in] = _depth[v_in] + 1;
   1.383        Cost sigma = _forward[u_in] ?
   1.384          _potential[v_in] - _potential[u_in] - _cost[_pred_edge[u_in]] :
   1.385 @@ -1145,9 +1132,9 @@
   1.386        }
   1.387      }
   1.388  
   1.389 -    /// Execute the algorithm.
   1.390 +    // Execute the algorithm
   1.391      bool start(PivotRuleEnum pivot_rule) {
   1.392 -      // Selecting the pivot rule implementation
   1.393 +      // Select the pivot rule implementation
   1.394        switch (pivot_rule) {
   1.395          case FIRST_ELIGIBLE_PIVOT:
   1.396            return start<FirstEligiblePivotRule>();
   1.397 @@ -1167,9 +1154,9 @@
   1.398      bool start() {
   1.399        PivotRuleImplementation pivot(*this, _edges);
   1.400  
   1.401 -      // Executing the network simplex algorithm
   1.402 +      // Execute the network simplex algorithm
   1.403        while (pivot.findEnteringEdge()) {
   1.404 -        join = findJoinNode();
   1.405 +        findJoinNode();
   1.406          bool change = findLeavingEdge();
   1.407          changeFlows(change);
   1.408          if (change) {
   1.409 @@ -1179,14 +1166,13 @@
   1.410          }
   1.411        }
   1.412  
   1.413 -      // Checking if the flow amount equals zero on all the artificial
   1.414 -      // edges
   1.415 +      // Check if the flow amount equals zero on all the artificial edges
   1.416        for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
   1.417          if (_flow[e] > 0) return false;
   1.418        for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
   1.419          if (_flow[e] > 0) return false;
   1.420  
   1.421 -      // Copying flow values to _flow_result
   1.422 +      // Copy flow values to _flow_result
   1.423        if (_lower) {
   1.424          for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
   1.425            (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
   1.426 @@ -1194,7 +1180,7 @@
   1.427          for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
   1.428            (*_flow_result)[e] = _flow[_edge_ref[e]];
   1.429        }
   1.430 -      // Copying potential values to _potential_result
   1.431 +      // Copy potential values to _potential_result
   1.432        for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
   1.433          (*_potential_result)[n] = _potential[_node_ref[n]];
   1.434