Changes in lemon/network_simplex.h [978:cbf32bf95954:1166:d59484d5fc1f] in lemon

File:

• lemon/network_simplex.h (modified) (38 diffs)

lemon/network_simplex.h

@@ -48 +48 @@
   /// flow problem.
   ///
-  /// In general, %NetworkSimplex is the fastest implementation available
-  /// in LEMON for this problem.
-  /// Moreover, it supports both directions of the supply/demand inequality
-  /// constraints. For more information, see \ref SupplyType.
+  /// In general, \ref NetworkSimplex and \ref CostScaling are the fastest
+  /// implementations available in LEMON for solving this problem.
+  /// (For more information, see \ref min_cost_flow_algs "the module page".)
+  /// Furthermore, this class supports both directions of the supply/demand
+  /// inequality constraints. For more information, see \ref SupplyType.
   ///
   /// Most of the parameters of the problem (except for the digraph)
@@ -64 +65 @@
   /// algorithm. By default, it is the same as \c V.
   ///
-  /// \warning Both number types must be signed and all input data must
+  /// \warning Both \c V and \c C must be signed number types.
+  /// \warning All input data (capacities, supply values, and costs) must
   /// be integer.
   ///
@@ -122 +124 @@
     /// the \ref run() function.
     ///
-    /// \ref NetworkSimplex provides five different pivot rule
-    /// implementations that significantly affect the running time
+    /// \ref NetworkSimplex provides five different implementations for
+    /// the pivot strategy that significantly affects the running time
     /// of the algorithm.
-    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
-    /// proved to be the most efficient and the most robust on various
-    /// test inputs.
-    /// However, another pivot rule can be selected using the \ref run()
-    /// function with the proper parameter.
+    /// According to experimental tests conducted on various problem
+    /// instances, \ref BLOCK_SEARCH "Block Search" and
+    /// \ref ALTERING_LIST "Altering Candidate List" rules turned out
+    /// to be the most efficient.
+    /// Since \ref BLOCK_SEARCH "Block Search" is a simpler strategy that
+    /// seemed to be slightly more robust, it is used by default.
+    /// However, another pivot rule can easily be selected using the
+    /// \ref run() function with the proper parameter.
     enum PivotRule {
 
@@ -155 +160 @@
       /// The \e Altering \e Candidate \e List pivot rule.
       /// It is a modified version of the Candidate List method.
-      /// It keeps only the several best eligible arcs from the former
+      /// It keeps only a few of the best eligible arcs from the former
       /// candidate list and extends this list in every iteration.
       ALTERING_LIST
@@ -167 +172 @@
     typedef std::vector<Value> ValueVector;
     typedef std::vector<Cost> CostVector;
-    typedef std::vector<char> BoolVector;
-    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
+    typedef std::vector<signed char> CharVector;
+    // Note: vector<signed char> is used instead of vector<ArcState> and
+    // vector<ArcDirection> for efficiency reasons
 
     // State constants for arcs
@@ -177 +183 @@
     };
 
-    typedef std::vector<signed char> StateVector;
-    // Note: vector<signed char> is used instead of vector<ArcState> for
-    // efficiency reasons
+    // Direction constants for tree arcs
+    enum ArcDirection {
+      DIR_DOWN = -1,
+      DIR_UP   =  1
+    };
 
   private:
@@ -218 +226 @@
     IntVector _succ_num;
     IntVector _last_succ;
+    CharVector _pred_dir;
+    CharVector _state;
     IntVector _dirty_revs;
-    BoolVector _forward;
-    StateVector _state;
     int _root;
 
     // Temporary data used in the current pivot iteration
     int in_arc, join, u_in, v_in, u_out, v_out;
-    int first, second, right, last;
-    int stem, par_stem, new_stem;
     Value delta;
 
@@ -251 +257 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const StateVector &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -303 +309 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const StateVector &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -342 +348 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const StateVector &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -415 +421 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const StateVector &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -518 +524 @@
       const IntVector  &_target;
       const CostVector &_cost;
-      const StateVector &_state;
+      const CharVector &_state;
       const CostVector &_pi;
       int &_in_arc;
@@ -537 +543 @@
         SortFunc(const CostVector &map) : _map(map) {}
         bool operator()(int left, int right) {
-          return _map[left] > _map[right];
+          return _map[left] < _map[right];
         }
       };
@@ -555 +561 @@
         const double BLOCK_SIZE_FACTOR = 1.0;
         const int MIN_BLOCK_SIZE = 10;
-        const double HEAD_LENGTH_FACTOR = 0.1;
+        const double HEAD_LENGTH_FACTOR = 0.01;
         const int MIN_HEAD_LENGTH = 3;
 
@@ -571 +577 @@
         // Check the current candidate list
         int e;
+        Cost c;
         for (int i = 0; i != _curr_length; ++i) {
           e = _candidates[i];
-          _cand_cost[e] = _state[e] *
-            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-          if (_cand_cost[e] >= 0) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < 0) {
+            _cand_cost[e] = c;
+          } else {
             _candidates[i--] = _candidates[--_curr_length];
           }
@@ -585 +593 @@
 
         for (e = _next_arc; e != _search_arc_num; ++e) {
-          _cand_cost[e] = _state[e] *
-            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-          if (_cand_cost[e] < 0) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < 0) {
+            _cand_cost[e] = c;
             _candidates[_curr_length++] = e;
           }
@@ -597 +605 @@
         }
         for (e = 0; e != _next_arc; ++e) {
-          _cand_cost[e] = _state[e] *
-            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
-          if (_cand_cost[e] < 0) {
+          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
+          if (c < 0) {
+            _cand_cost[e] = c;
             _candidates[_curr_length++] = e;
           }
@@ -612 +620 @@
       search_end:
 
-        // 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
+        // Perform partial sort operation on the candidate list
+        int new_length = std::min(_head_length + 1, _curr_length);
+        std::partial_sort(_candidates.begin(), _candidates.begin() + new_length,
+                          _candidates.begin() + _curr_length, _sort_func);
+
+        // Select the entering arc and remove it from the list
         _in_arc = _candidates[0];
         _next_arc = e;
-        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
-                  _sort_func );
-        _curr_length = std::min(_head_length, _curr_length - 1);
+        _candidates[0] = _candidates[new_length - 1];
+        _curr_length = new_length - 1;
         return true;
       }
@@ -634 +642 @@
     ///
     /// \param graph The digraph the algorithm runs on.
-    /// \param arc_mixing Indicate if the arcs have to be stored in a
+    /// \param arc_mixing Indicate if the arcs will be stored in a
     /// mixed order in the internal data structure.
-    /// In special cases, it could lead to better overall performance,
-    /// but it is usually slower. Therefore it is disabled by default.
-    NetworkSimplex(const GR& graph, bool arc_mixing = false) :
+    /// In general, it leads to similar performance as using the original
+    /// arc order, but it makes the algorithm more robust and in special
+    /// cases, even significantly faster. Therefore, it is enabled by default.
+    NetworkSimplex(const GR& graph, bool arc_mixing = true) :
       _graph(graph), _node_id(graph), _arc_id(graph),
       _arc_mixing(arc_mixing),
@@ -731 +740 @@
     ///
     /// \return <tt>(*this)</tt>
+    ///
+    /// \sa supplyType()
     template<typename SupplyMap>
     NetworkSimplex& supplyMap(const SupplyMap& map) {
@@ -747 +758 @@
     ///
     /// Using this function has the same effect as using \ref supplyMap()
-    /// with such a map in which \c k is assigned to \c s, \c -k is
+    /// with a map in which \c k is assigned to \c s, \c -k is
     /// assigned to \c t and all other nodes have zero supply value.
     ///
@@ -914 +925 @@
       _parent.resize(all_node_num);
       _pred.resize(all_node_num);
-      _forward.resize(all_node_num);
+      _pred_dir.resize(all_node_num);
       _thread.resize(all_node_num);
       _rev_thread.resize(all_node_num);
@@ -928 +939 @@
       if (_arc_mixing) {
         // Store the arcs in a mixed order
-        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
+        const int skip = std::max(_arc_num / _node_num, 3);
         int i = 0, j = 0;
         for (ArcIt a(_graph); a != INVALID; ++a) {
@@ -934 +945 @@
           _source[i] = _node_id[_graph.source(a)];
           _target[i] = _node_id[_graph.target(a)];
-          if ((i += k) >= _arc_num) i = ++j;
+          if ((i += skip) >= _arc_num) i = ++j;
         }
       } else {
@@ -1000 +1011 @@
     }
 
-    /// \brief Return the flow map (the primal solution).
+    /// \brief Copy the flow values (the primal solution) into the
+    /// given map.
     ///
     /// This function copies the flow value on each arc into the given
@@ -1024 +1036 @@
     }
 
-    /// \brief Return the potential map (the dual solution).
+    /// \brief Copy the potential values (the dual solution) into the
+    /// given map.
     ///
     /// This function copies the potential (dual value) of each node
@@ -1117 +1130 @@
           _state[e] = STATE_TREE;
           if (_supply[u] >= 0) {
-            _forward[u] = true;
+            _pred_dir[u] = DIR_UP;
             _pi[u] = 0;
             _source[e] = u;
@@ -1124 +1137 @@
             _cost[e] = 0;
           } else {
-            _forward[u] = false;
+            _pred_dir[u] = DIR_DOWN;
             _pi[u] = ART_COST;
             _source[e] = _root;
@@ -1144 +1157 @@
           _last_succ[u] = u;
           if (_supply[u] >= 0) {
-            _forward[u] = true;
+            _pred_dir[u] = DIR_UP;
             _pi[u] = 0;
             _pred[u] = e;
@@ -1154 +1167 @@
             _state[e] = STATE_TREE;
           } else {
-            _forward[u] = false;
+            _pred_dir[u] = DIR_DOWN;
             _pi[u] = ART_COST;
             _pred[u] = f;
@@ -1185 +1198 @@
           _last_succ[u] = u;
           if (_supply[u] <= 0) {
-            _forward[u] = false;
+            _pred_dir[u] = DIR_DOWN;
             _pi[u] = 0;
             _pred[u] = e;
@@ -1195 +1208 @@
             _state[e] = STATE_TREE;
           } else {
-            _forward[u] = true;
+            _pred_dir[u] = DIR_UP;
             _pi[u] = -ART_COST;
             _pred[u] = f;
@@ -1238 +1251 @@
       // Initialize first and second nodes according to the direction
       // of the cycle
+      int first, second;
       if (_state[in_arc] == STATE_LOWER) {
         first  = _source[in_arc];
@@ -1247 +1261 @@
       delta = _cap[in_arc];
       int result = 0;
-      Value d;
+      Value c, d;
       int e;
 
-      // Search the cycle along the path form the first node to the root
+      // Search the cycle form the first node to the join node
       for (int u = first; u != join; u = _parent[u]) {
         e = _pred[u];
-        d = _forward[u] ?
-          _flow[e] : (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]);
+        d = _flow[e];
+        if (_pred_dir[u] == DIR_DOWN) {
+          c = _cap[e];
+          d = c >= MAX ? INF : c - d;
+        }
         if (d < delta) {
           delta = d;
@@ -1261 +1278 @@
         }
       }
-      // Search the cycle along the path form the second node to the root
+
+      // Search the cycle form the second node to the join node
       for (int u = second; u != join; u = _parent[u]) {
         e = _pred[u];
-        d = _forward[u] ?
-          (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
+        d = _flow[e];
+        if (_pred_dir[u] == DIR_UP) {
+          c = _cap[e];
+          d = c >= MAX ? INF : c - d;
+        }
         if (d <= delta) {
           delta = d;
@@ -1290 +1311 @@
         _flow[in_arc] += val;
         for (int u = _source[in_arc]; u != join; u = _parent[u]) {
-          _flow[_pred[u]] += _forward[u] ? -val : val;
+          _flow[_pred[u]] -= _pred_dir[u] * val;
         }
         for (int u = _target[in_arc]; u != join; u = _parent[u]) {
-          _flow[_pred[u]] += _forward[u] ? val : -val;
+          _flow[_pred[u]] += _pred_dir[u] * val;
         }
       }
@@ -1308 +1329 @@
     // Update the tree structure
     void updateTreeStructure() {
-      int u, w;
       int old_rev_thread = _rev_thread[u_out];
       int old_succ_num = _succ_num[u_out];
@@ -1314 +1334 @@
       v_out = _parent[u_out];
 
-      u = _last_succ[u_in];  // the last successor of u_in
-      right = _thread[u];    // the node after it
-
-      // Handle the case when old_rev_thread equals to v_in
-      // (it also means that join and v_out coincide)
-      if (old_rev_thread == v_in) {
-        last = _thread[_last_succ[u_out]];
+      // Check if u_in and u_out coincide
+      if (u_in == u_out) {
+        // Update _parent, _pred, _pred_dir
+        _parent[u_in] = v_in;
+        _pred[u_in] = in_arc;
+        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
+
+        // Update _thread and _rev_thread
+        if (_thread[v_in] != u_out) {
+          int after = _thread[old_last_succ];
+          _thread[old_rev_thread] = after;
+          _rev_thread[after] = old_rev_thread;
+          after = _thread[v_in];
+          _thread[v_in] = u_out;
+          _rev_thread[u_out] = v_in;
+          _thread[old_last_succ] = after;
+          _rev_thread[after] = old_last_succ;
+        }
       } else {
-        last = _thread[v_in];
-      }
-
-      // Update _thread and _parent along the stem nodes (i.e. the nodes
-      // between u_in and u_out, whose parent have to be changed)
-      _thread[v_in] = stem = u_in;
-      _dirty_revs.clear();
-      _dirty_revs.push_back(v_in);
-      par_stem = v_in;
-      while (stem != u_out) {
-        // Insert the next stem node into the thread list
-        new_stem = _parent[stem];
-        _thread[u] = new_stem;
-        _dirty_revs.push_back(u);
-
-        // Remove the subtree of stem from the thread list
-        w = _rev_thread[stem];
-        _thread[w] = right;
-        _rev_thread[right] = w;
-
-        // Change the parent node and shift stem nodes
-        _parent[stem] = par_stem;
-        par_stem = stem;
-        stem = new_stem;
-
-        // Update u and right
-        u = _last_succ[stem] == _last_succ[par_stem] ?
-          _rev_thread[par_stem] : _last_succ[stem];
-        right = _thread[u];
-      }
-      _parent[u_out] = par_stem;
-      _thread[u] = last;
-      _rev_thread[last] = u;
-      _last_succ[u_out] = u;
-
-      // Remove the subtree of u_out from the thread list except for
-      // the case when old_rev_thread equals to v_in
-      // (it also means that join and v_out coincide)
-      if (old_rev_thread != v_in) {
-        _thread[old_rev_thread] = right;
-        _rev_thread[right] = old_rev_thread;
-      }
-
-      // Update _rev_thread using the new _thread values
-      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
-        u = _dirty_revs[i];
-        _rev_thread[_thread[u]] = u;
-      }
-
-      // Update _pred, _forward, _last_succ and _succ_num for the
-      // stem nodes from u_out to u_in
-      int tmp_sc = 0, tmp_ls = _last_succ[u_out];
-      u = u_out;
-      while (u != u_in) {
-        w = _parent[u];
-        _pred[u] = _pred[w];
-        _forward[u] = !_forward[w];
-        tmp_sc += _succ_num[u] - _succ_num[w];
-        _succ_num[u] = tmp_sc;
-        _last_succ[w] = tmp_ls;
-        u = w;
-      }
-      _pred[u_in] = in_arc;
-      _forward[u_in] = (u_in == _source[in_arc]);
-      _succ_num[u_in] = old_succ_num;
-
-      // Set limits for updating _last_succ form v_in and v_out
-      // towards the root
-      int up_limit_in = -1;
-      int up_limit_out = -1;
-      if (_last_succ[join] == v_in) {
-        up_limit_out = join;
-      } else {
-        up_limit_in = join;
+        // Handle the case when old_rev_thread equals to v_in
+        // (it also means that join and v_out coincide)
+        int thread_continue = old_rev_thread == v_in ?
+          _thread[old_last_succ] : _thread[v_in];
+
+        // Update _thread and _parent along the stem nodes (i.e. the nodes
+        // between u_in and u_out, whose parent have to be changed)
+        int stem = u_in;              // the current stem node
+        int par_stem = v_in;          // the new parent of stem
+        int next_stem;                // the next stem node
+        int last = _last_succ[u_in];  // the last successor of stem
+        int before, after = _thread[last];
+        _thread[v_in] = u_in;
+        _dirty_revs.clear();
+        _dirty_revs.push_back(v_in);
+        while (stem != u_out) {
+          // Insert the next stem node into the thread list
+          next_stem = _parent[stem];
+          _thread[last] = next_stem;
+          _dirty_revs.push_back(last);
+
+          // Remove the subtree of stem from the thread list
+          before = _rev_thread[stem];
+          _thread[before] = after;
+          _rev_thread[after] = before;
+
+          // Change the parent node and shift stem nodes
+          _parent[stem] = par_stem;
+          par_stem = stem;
+          stem = next_stem;
+
+          // Update last and after
+          last = _last_succ[stem] == _last_succ[par_stem] ?
+            _rev_thread[par_stem] : _last_succ[stem];
+          after = _thread[last];
+        }
+        _parent[u_out] = par_stem;
+        _thread[last] = thread_continue;
+        _rev_thread[thread_continue] = last;
+        _last_succ[u_out] = last;
+
+        // Remove the subtree of u_out from the thread list except for
+        // the case when old_rev_thread equals to v_in
+        if (old_rev_thread != v_in) {
+          _thread[old_rev_thread] = after;
+          _rev_thread[after] = old_rev_thread;
+        }
+
+        // Update _rev_thread using the new _thread values
+        for (int i = 0; i != int(_dirty_revs.size()); ++i) {
+          int u = _dirty_revs[i];
+          _rev_thread[_thread[u]] = u;
+        }
+
+        // Update _pred, _pred_dir, _last_succ and _succ_num for the
+        // stem nodes from u_out to u_in
+        int tmp_sc = 0, tmp_ls = _last_succ[u_out];
+        for (int u = u_out, p = _parent[u]; u != u_in; u = p, p = _parent[u]) {
+          _pred[u] = _pred[p];
+          _pred_dir[u] = -_pred_dir[p];
+          tmp_sc += _succ_num[u] - _succ_num[p];
+          _succ_num[u] = tmp_sc;
+          _last_succ[p] = tmp_ls;
+        }
+        _pred[u_in] = in_arc;
+        _pred_dir[u_in] = u_in == _source[in_arc] ? DIR_UP : DIR_DOWN;
+        _succ_num[u_in] = old_succ_num;
       }
 
       // Update _last_succ from v_in towards the root
-      for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
-           u = _parent[u]) {
-        _last_succ[u] = _last_succ[u_out];
-      }
+      int up_limit_out = _last_succ[join] == v_in ? join : -1;
+      int last_succ_out = _last_succ[u_out];
+      for (int u = v_in; u != -1 && _last_succ[u] == v_in; u = _parent[u]) {
+        _last_succ[u] = last_succ_out;
+      }
+
       // Update _last_succ from v_out towards the root
       if (join != old_rev_thread && v_in != old_rev_thread) {
-        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
+        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
              u = _parent[u]) {
           _last_succ[u] = old_rev_thread;
         }
-      } else {
-        for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
+      }
+      else if (last_succ_out != old_last_succ) {
+        for (int u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
              u = _parent[u]) {
-          _last_succ[u] = _last_succ[u_out];
+          _last_succ[u] = last_succ_out;
         }
       }
 
       // Update _succ_num from v_in to join
-      for (u = v_in; u != join; u = _parent[u]) {
+      for (int u = v_in; u != join; u = _parent[u]) {
         _succ_num[u] += old_succ_num;
       }
       // Update _succ_num from v_out to join
-      for (u = v_out; u != join; u = _parent[u]) {
+      for (int u = v_out; u != join; u = _parent[u]) {
         _succ_num[u] -= old_succ_num;
       }
     }
 
-    // Update potentials
+    // Update potentials in the subtree that has been moved
     void updatePotential() {
-      Cost sigma = _forward[u_in] ?
-        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
-        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
-      // Update potentials in the subtree, which has been moved
+      Cost sigma = _pi[v_in] - _pi[u_in] -
+                   _pred_dir[u_in] * _cost[in_arc];
       int end = _thread[_last_succ[u_in]];
       for (int u = u_in; u != end; u = _thread[u]) {