Ticket #179: 179-3-implementation-83ce7ce39f21.patch

lemon/min_mean_cycle.h

@@ -74 +74 @@
     // The length of the arcs
     const LengthMap &_length;
 
-    // The total length of the found cycle
-    Value _cycle_length;
-    // The number of arcs on the found cycle
-    int _cycle_size;
-    // The found cycle
+    // Data for the found cycles
+    bool _curr_found, _best_found;
+    Value _curr_length, _best_length;
+    int _curr_size, _best_size;
+    Node _curr_node, _best_node;
+
     Path *_cycle_path;
+    bool _local_path;
 
-    bool _local_path;
-    bool _cycle_found;
-    Node _cycle_node;
+    // Internal data used by the algorithm
+    typename Digraph::template NodeMap<Arc> _policy;
+    typename Digraph::template NodeMap<bool> _reached;
+    typename Digraph::template NodeMap<int> _level;
+    typename Digraph::template NodeMap<double> _dist;
 
-    typename Digraph::template NodeMap<bool> _reached;
-    typename Digraph::template NodeMap<double> _dist;
-    typename Digraph::template NodeMap<Arc> _policy;
-
+    // Data for storing the strongly connected components
+    int _comp_num;
     typename Digraph::template NodeMap<int> _comp;
-    int _comp_num;
-
-    std::vector<Node> _nodes;
-    std::vector<Arc> _arcs;
+    std::vector<std::vector<Node> > _comp_nodes;
+    std::vector<Node>* _nodes;
+    typename Digraph::template NodeMap<std::vector<Arc> > _in_arcs;
+    
+    // Queue used for BFS search
+    std::vector<Node> _queue;
+    int _qfront, _qback;
+    
     Tolerance<double> _tol;
 
   public:
@@ -106 +112 @@
     /// \param length The lengths (costs) of the arcs.
     MinMeanCycle( const Digraph &digraph,
                   const LengthMap &length ) :
-      _gr(digraph), _length(length), _cycle_length(0), _cycle_size(-1),
-      _cycle_path(NULL), _local_path(false), _reached(digraph),
-      _dist(digraph), _policy(digraph), _comp(digraph)
+      _gr(digraph), _length(length), _cycle_path(NULL), _local_path(false),
+      _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph),
+      _comp(digraph), _in_arcs(digraph)
     {}
 
     /// Destructor.
@@ -172 +178 @@
     ///
     /// \return \c true if a directed cycle exists in the digraph.
     bool findMinMean() {
-      // Initialize
-      _tol.epsilon(1e-6);
-      if (!_cycle_path) {
-        _local_path = true;
-        _cycle_path = new Path;
-      }
-      _cycle_path->clear();
-      _cycle_found = false;
-
+      // Initialize and find strongly connected components
+      init();
+      findComponents();
+      
       // Find the minimum cycle mean in the components
-      _comp_num = stronglyConnectedComponents(_gr, _comp);
       for (int comp = 0; comp < _comp_num; ++comp) {
-        if (!initCurrentComponent(comp)) continue;
+        // Find the minimum mean cycle in the current component
+        if (!buildPolicyGraph(comp)) continue;
         while (true) {
-          if (!findPolicyCycles()) break;
-          contractPolicyGraph(comp);
+          findPolicyCycle();
           if (!computeNodeDistances()) break;
         }
+        // Update the best cycle (global minimum mean cycle)
+        if ( !_best_found || (_curr_found &&
+             _curr_length * _best_size < _best_length * _curr_size) ) {
+          _best_found = true;
+          _best_length = _curr_length;
+          _best_size = _curr_size;
+          _best_node = _curr_node;
+        }
       }
-      return _cycle_found;
+      return _best_found;
     }
 
     /// \brief Find a minimum mean directed cycle.
@@ -203 +211 @@
     ///
     /// \pre \ref findMinMean() must be called before using this function.
     bool findCycle() {
-      if (!_cycle_found) return false;
-      _cycle_path->addBack(_policy[_cycle_node]);
-      for ( Node v = _cycle_node;
-            (v = _gr.target(_policy[v])) != _cycle_node; ) {
+      if (!_best_found) return false;
+      _cycle_path->addBack(_policy[_best_node]);
+      for ( Node v = _best_node;
+            (v = _gr.target(_policy[v])) != _best_node; ) {
         _cycle_path->addBack(_policy[v]);
       }
       return true;
@@ -225 +233 @@
     ///
     /// This function returns the total length of the found cycle.
     ///
-    /// \pre \ref run() or \ref findCycle() must be called before
+    /// \pre \ref run() or \ref findMinMean() must be called before
     /// using this function.
     Value cycleLength() const {
-      return _cycle_length;
+      return _best_length;
     }
 
     /// \brief Return the number of arcs on the found cycle.
     ///
     /// This function returns the number of arcs on the found cycle.
     ///
-    /// \pre \ref run() or \ref findCycle() must be called before
+    /// \pre \ref run() or \ref findMinMean() must be called before
     /// using this function.
     int cycleArcNum() const {
-      return _cycle_size;
+      return _best_size;
     }
 
     /// \brief Return the mean length of the found cycle.
     ///
     /// This function returns the mean length of the found cycle.
     ///
-    /// \note <tt>mmc.cycleMean()</tt> is just a shortcut of the
+    /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
     /// following code.
     /// \code
-    ///   return double(mmc.cycleLength()) / mmc.cycleArcNum();
+    ///   return static_cast<double>(alg.cycleLength()) / alg.cycleArcNum();
     /// \endcode
     ///
     /// \pre \ref run() or \ref findMinMean() must be called before
     /// using this function.
     double cycleMean() const {
-      return double(_cycle_length) / _cycle_size;
+      return static_cast<double>(_best_length) / _best_size;
     }
 
     /// \brief Return the found cycle.
@@ -274 +282 @@
 
   private:
 
-    // Initialize the internal data structures for the current strongly
-    // connected component and create the policy graph.
-    // The policy graph can be represented by the _policy map because
-    // the out-degree of every node is 1.
-    bool initCurrentComponent(int comp) {
-      // Find the nodes of the current component
-      _nodes.clear();
-      for (NodeIt n(_gr); n != INVALID; ++n) {
-        if (_comp[n] == comp) _nodes.push_back(n);
+    // Initialize
+    void init() {
+      _tol.epsilon(1e-6);
+      if (!_cycle_path) {
+        _local_path = true;
+        _cycle_path = new Path;
       }
-      if (_nodes.size() <= 1) return false;
-      // Find the arcs of the current component
-      _arcs.clear();
-      for (ArcIt e(_gr); e != INVALID; ++e) {
-        if ( _comp[_gr.source(e)] == comp &&
-             _comp[_gr.target(e)] == comp )
-          _arcs.push_back(e);
+      _queue.resize(countNodes(_gr));
+      _best_found = false;
+      _best_length = 0;
+      _best_size = 1;
+      _cycle_path->clear();
+    }
+    
+    // Find strongly connected components and initialize _comp_nodes
+    // and _in_arcs
+    void findComponents() {
+      _comp_num = stronglyConnectedComponents(_gr, _comp);
+      _comp_nodes.resize(_comp_num);
+      if (_comp_num == 1) {
+        _comp_nodes[0].clear();
+        for (NodeIt n(_gr); n != INVALID; ++n) {
+          _comp_nodes[0].push_back(n);
+          _in_arcs[n].clear();
+          for (InArcIt a(_gr, n); a != INVALID; ++a) {
+            _in_arcs[n].push_back(a);
+          }
+        }
+      } else {
+        for (int i = 0; i < _comp_num; ++i)
+          _comp_nodes[i].clear();
+        for (NodeIt n(_gr); n != INVALID; ++n) {
+          int k = _comp[n];
+          _comp_nodes[k].push_back(n);
+          _in_arcs[n].clear();
+          for (InArcIt a(_gr, n); a != INVALID; ++a) {
+            if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a);
+          }
+        }
       }
-      // Initialize _reached, _dist, _policy maps
-      for (int i = 0; i < int(_nodes.size()); ++i) {
-        _reached[_nodes[i]] = false;
-        _policy[_nodes[i]] = INVALID;
+    }
+
+    // Build the policy graph in the given strongly connected component
+    // (the out-degree of every node is 1)
+    bool buildPolicyGraph(int comp) {
+      _nodes = &(_comp_nodes[comp]);
+      if (_nodes->size() < 1 ||
+          (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) {
+        return false;
       }
-      Node u; Arc e;
-      for (int j = 0; j < int(_arcs.size()); ++j) {
-        e = _arcs[j];
-        u = _gr.source(e);
-        if (!_reached[u] || _length[e] < _dist[u]) {
-          _dist[u] = _length[e];
-          _policy[u] = e;
-          _reached[u] = true;
+      for (int i = 0; i < int(_nodes->size()); ++i) {
+        _dist[(*_nodes)[i]] = std::numeric_limits<double>::max();
+      }
+      Node u, v;
+      Arc e;
+      for (int i = 0; i < int(_nodes->size()); ++i) {
+        v = (*_nodes)[i];
+        for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
+          e = _in_arcs[v][j];
+          u = _gr.source(e);
+          if (_length[e] < _dist[u]) {
+            _dist[u] = _length[e];
+            _policy[u] = e;
+          }
         }
       }
       return true;
     }
 
-    // Find all cycles in the policy graph.
-    // Set _cycle_found to true if a cycle is found and set
-    // _cycle_length, _cycle_size, _cycle_node to represent the minimum
-    // mean cycle in the policy graph.
-    bool findPolicyCycles() {
-      typename Digraph::template NodeMap<int> level(_gr, -1);
-      bool curr_cycle_found = false;
+    // Find the minimum mean cycle in the policy graph
+    void findPolicyCycle() {
+      for (int i = 0; i < int(_nodes->size()); ++i) {
+        _level[(*_nodes)[i]] = -1;
+      }
       Value clength;
       int csize;
-      int path_cnt = 0;
       Node u, v;
-      // Searching for cycles
-      for (int i = 0; i < int(_nodes.size()); ++i) {
-        if (level[_nodes[i]] < 0) {
-          u = _nodes[i];
-          level[u] = path_cnt;
-          while (level[u = _gr.target(_policy[u])] < 0)
-            level[u] = path_cnt;
-          if (level[u] == path_cnt) {
-            // A cycle is found
-            curr_cycle_found = true;
-            clength = _length[_policy[u]];
-            csize = 1;
-            for (v = u; (v = _gr.target(_policy[v])) != u; ) {
-              clength += _length[_policy[v]];
-              ++csize;
-            }
-            if ( !_cycle_found ||
-                 clength * _cycle_size < _cycle_length * csize ) {
-              _cycle_found = true;
-              _cycle_length = clength;
-              _cycle_size = csize;
-              _cycle_node = u;
-            }
+      _curr_found = false;
+      for (int i = 0; i < int(_nodes->size()); ++i) {
+        u = (*_nodes)[i];
+        if (_level[u] >= 0) continue;
+        for (; _level[u] < 0; u = _gr.target(_policy[u])) {
+          _level[u] = i;
+        }
+        if (_level[u] == i) {
+          // A cycle is found
+          clength = _length[_policy[u]];
+          csize = 1;
+          for (v = u; (v = _gr.target(_policy[v])) != u; ) {
+            clength += _length[_policy[v]];
+            ++csize;
           }
-          ++path_cnt;
-        }
-      }
-      return curr_cycle_found;
-    }
-
-    // Contract the policy graph to be connected by cutting all cycles
-    // except for the main cycle (i.e. the minimum mean cycle).
-    void contractPolicyGraph(int comp) {
-      // Find the component of the main cycle using reverse BFS search
-      typename Digraph::template NodeMap<int> found(_gr, false);
-      std::deque<Node> queue;
-      queue.push_back(_cycle_node);
-      found[_cycle_node] = true;
-      Node u, v;
-      while (!queue.empty()) {
-        v = queue.front(); queue.pop_front();
-        for (InArcIt e(_gr, v); e != INVALID; ++e) {
-          u = _gr.source(e);
-          if (_policy[u] == e && !found[u]) {
-            found[u] = true;
-            queue.push_back(u);
-          }
-        }
-      }
-      // Connect all other nodes to this component using reverse BFS search
-      queue.clear();
-      for (int i = 0; i < int(_nodes.size()); ++i)
-        if (found[_nodes[i]]) queue.push_back(_nodes[i]);
-      int found_cnt = queue.size();
-      while (found_cnt < int(_nodes.size())) {
-        v = queue.front(); queue.pop_front();
-        for (InArcIt e(_gr, v); e != INVALID; ++e) {
-          u = _gr.source(e);
-          if (_comp[u] == comp && !found[u]) {
-            found[u] = true;
-            ++found_cnt;
-            _policy[u] = e;
-            queue.push_back(u);
+          if ( !_curr_found ||
+               (clength * _curr_size < _curr_length * csize) ) {
+            _curr_found = true;
+            _curr_length = clength;
+            _curr_size = csize;
+            _curr_node = u;
           }
         }
       }
     }
 
-    // Compute node distances in the policy graph and update the
-    // policy graph if the node distances can be improved.
+    // Contract the policy graph and compute node distances
     bool computeNodeDistances() {
-      // Compute node distances using reverse BFS search
-      double cycle_mean = double(_cycle_length) / _cycle_size;
-      typename Digraph::template NodeMap<int> found(_gr, false);
-      std::deque<Node> queue;
-      queue.push_back(_cycle_node);
-      found[_cycle_node] = true;
-      _dist[_cycle_node] = 0;
+      // Find the component of the main cycle and compute node distances
+      // using reverse BFS
+      for (int i = 0; i < int(_nodes->size()); ++i) {
+        _reached[(*_nodes)[i]] = false;
+      }
+      double curr_mean = double(_curr_length) / _curr_size;
+      _qfront = _qback = 0;
+      _queue[0] = _curr_node;
+      _reached[_curr_node] = true;
+      _dist[_curr_node] = 0;
       Node u, v;
-      while (!queue.empty()) {
-        v = queue.front(); queue.pop_front();
-        for (InArcIt e(_gr, v); e != INVALID; ++e) {
+      Arc e;
+      while (_qfront <= _qback) {
+        v = _queue[_qfront++];
+        for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
+          e = _in_arcs[v][j];
           u = _gr.source(e);
-          if (_policy[u] == e && !found[u]) {
-            found[u] = true;
-            _dist[u] = _dist[v] + _length[e] - cycle_mean;
-            queue.push_back(u);
+          if (_policy[u] == e && !_reached[u]) {
+            _reached[u] = true;
+            _dist[u] = _dist[v] + _length[e] - curr_mean;
+            _queue[++_qback] = u;
           }
         }
       }
-      // Improving node distances
+
+      // Connect all other nodes to this component and compute node
+      // distances using reverse BFS
+      _qfront = 0;
+      while (_qback < int(_nodes->size())-1) {
+        v = _queue[_qfront++];
+        for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
+          e = _in_arcs[v][j];
+          u = _gr.source(e);
+          if (!_reached[u]) {
+            _reached[u] = true;
+            _policy[u] = e;
+            _dist[u] = _dist[v] + _length[e] - curr_mean;
+            _queue[++_qback] = u;
+          }
+        }
+      }
+
+      // Improve node distances
       bool improved = false;
-      for (int j = 0; j < int(_arcs.size()); ++j) {
-        Arc e = _arcs[j];
-        u = _gr.source(e); v = _gr.target(e);
-        double delta = _dist[v] + _length[e] - cycle_mean;
-        if (_tol.less(delta, _dist[u])) {
-          improved = true;
-          _dist[u] = delta;
-          _policy[u] = e;
+      for (int i = 0; i < int(_nodes->size()); ++i) {
+        v = (*_nodes)[i];
+        for (int j = 0; j < int(_in_arcs[v].size()); ++j) {
+          e = _in_arcs[v][j];
+          u = _gr.source(e);
+          double delta = _dist[v] + _length[e] - curr_mean;
+          if (_tol.less(delta, _dist[u])) {
+            _dist[u] = delta;
+            _policy[u] = e;
+            improved = true;
+          }
         }
       }
       return improved;