/// \brief Karp's algorithm for finding a minimum mean (directed) cycle.

 /// \brief Howard's algorithm for finding a minimum mean directed cycle.

   /// \brief Implementation of Karp's algorithm for finding a

   /// \brief Implementation of Howard's algorithm for finding a minimum

   /// minimum mean (directed) cycle.

   /// mean directed cycle.

   /// The \ref lemon::MinMeanCycle "MinMeanCycle" implements Karp's

   /// \ref MinMeanCycle implements Howard's algorithm for finding a

   /// algorithm for finding a minimum mean (directed) cycle.

   /// minimum mean directed cycle.

   template <typename Graph,

   template < typename Graph,

     typename LengthMap = typename Graph::template EdgeMap<int> >

              typename LengthMap = typename Graph::template EdgeMap<int> >

     typedef typename Graph::Node Node;

     GRAPH_TYPEDEFS(typename Graph);

     typedef typename Graph::NodeIt NodeIt;

     typedef typename Graph::Edge Edge;

     typedef typename Graph::EdgeIt EdgeIt;

     typedef typename Graph::OutEdgeIt OutEdgeIt;

     /// \brief Data structure for path data.

     typename Graph::template NodeMap<bool> _reached;

     struct PathData

     typename Graph::template NodeMap<double> _dist;

     {

     typename Graph::template NodeMap<Edge> _policy;

       bool found;

       Length dist;

     /// The directed graph the algorithm runs on.

       Edge pred;

     const Graph &_graph;

       PathData(bool _found = false, Length _dist = 0) :

     /// The length of the edges.

 	found(_found), dist(_dist), pred(INVALID) {}

     const LengthMap &_length;

       PathData(bool _found, Length _dist, Edge _pred) :

 	found(_found), dist(_dist), pred(_pred) {}

     /// The total length of the found cycle.

     };

     Length _cycle_length;

     /// The number of edges on the found cycle.

   private:

     int _cycle_size;

     /// The found cycle.

     typedef typename Graph::template NodeMap<std::vector<PathData> >

     Path *_cycle_path;

       PathDataNodeMap;

     bool _local_path;

     bool _cycle_found;

     Node _cycle_node;

     typename Graph::template NodeMap<int> _component;

     int _component_num;

     std::vector<Node> _nodes;

     std::vector<Edge> _edges;

     Tolerance<double> _tolerance;

   public:

     /// \brief The constructor of the class.

     ///

     /// The constructor of the class.

     ///

     /// \param graph The directed graph the algorithm runs on.

     /// \param length The length (cost) of the edges.

     MinMeanCycle( const Graph &graph,

                   const LengthMap &length ) :

       _graph(graph), _length(length), _dist(graph), _reached(graph),

       _policy(graph), _component(graph), _cycle_length(0),

       _cycle_size(-1), _cycle_path(NULL), _local_path(false)

     { }

     /// The destructor of the class.

     ~MinMeanCycle() {

       if (_local_path) delete _cycle_path;

     }

     /// \brief Node map for storing path data.

     // Initializes the internal data structures for the current strongly

     ///

     // connected component and creating the policy graph.

     /// Node map for storing path data of all nodes in the current

     // The policy graph can be represented by the _policy map because

     /// component. dmap[v][k] is the length of a shortest directed walk

     // the out degree of every node is 1.

     /// to node v from the starting node containing exactly k edges.

     bool initCurrentComponent(int comp) {

     PathDataNodeMap dmap;

     /// \brief The directed graph the algorithm runs on.

     const Graph &graph;

     /// \brief The length of the edges.

     const LengthMap &length;

     /// \brief The total length of the found cycle.

     Length cycle_length;

     /// \brief The number of edges in the found cycle.

     int cycle_size;

     /// \brief A node for obtaining a minimum mean cycle.

     Node cycle_node;

     /// \brief The found cycle.

     Path *cycle_path;

     /// \brief The algorithm uses local \ref lemon::Path "Path"

     /// structure to store the found cycle.

     bool local_path;

     /// \brief Node map for identifying strongly connected components.

     IntNodeMap comp;

     /// \brief The number of strongly connected components.

     int comp_num;

     /// \brief Counter for identifying the current component.

     int comp_cnt;

     /// \brief Nodes of the current component.

     NodeVector nodes;

     /// \brief The processed nodes in the last round.

     NodeVector process;

   public :

     /// \brief The constructor of the class.

     ///

     /// The constructor of the class.

     ///

     /// \param _graph The directed graph the algorithm runs on.

     /// \param _length The length (cost) of the edges.

     MinMeanCycle( const Graph &_graph,

 		  const LengthMap &_length ) :

       graph(_graph), length(_length), dmap(_graph), comp(_graph),

       cycle_length(0), cycle_size(-1), cycle_node(INVALID),

       cycle_path(NULL), local_path(false)

     { }

     /// \brief The destructor of the class.

     ~MinMeanCycle() {

       if (local_path) delete cycle_path;

     }

   protected:

     /// \brief Initializes the internal data structures for the current

     /// component.

     void initCurrent() {

       nodes.clear();

       for (NodeIt v(graph); v != INVALID; ++v) {

       _nodes.clear();

 	if (comp[v] == comp_cnt) nodes.push_back(v);

       for (NodeIt n(_graph); n != INVALID; ++n) {

       }

         if (_component[n] == comp) _nodes.push_back(n);

       // Creating vectors for all nodes

       }

       int n = nodes.size();

       if (_nodes.size() <= 1) return false;

       for (int i = 0; i < n; ++i) {

       // Finding the edges of the current component

 	dmap[nodes[i]].resize(n + 1);

       _edges.clear();

       }

       for (EdgeIt e(_graph); e != INVALID; ++e) {

     }

         if ( _component[_graph.source(e)] == comp &&

              _component[_graph.target(e)] == comp )

     /// \brief Processes all rounds of computing required path data for

           _edges.push_back(e);

     /// the current component.

       }

     void processRounds() {

       // Initializing _reached, _dist, _policy maps

       dmap[nodes[0]][0] = PathData(true, 0);

       for (int i = 0; i < _nodes.size(); ++i) {

       process.clear();

         _reached[_nodes[i]] = false;

       // Processing the first round

         _policy[_nodes[i]] = INVALID;

       for (OutEdgeIt e(graph, nodes[0]); e != INVALID; ++e) {

       }

 	Node v = graph.target(e);

       Node u; Edge e;

 	if (comp[v] != comp_cnt || v == nodes[0]) continue;

       for (int j = 0; j < _edges.size(); ++j) {

 	dmap[v][1] = PathData(true, length[e], e);

         e = _edges[j];

 	process.push_back(v);

         u = _graph.source(e);

       }

         if (!_reached[u] || _length[e] < _dist[u]) {

       // Processing other rounds

           _dist[u] = _length[e];

       int n = nodes.size(), k;

           _policy[u] = e;

       for (k = 2; k <= n && process.size() < n; ++k)

           _reached[u] = true;

 	processNextBuildRound(k);

         }

       for ( ; k <= n; ++k)

       }

 	processNextFullRound(k);

       return true;

     }

     }

     /// \brief Processes one round of computing required path data and

     // Finds all cycles in the policy graph.

     /// rebuilds \ref process vector.

     // Sets _cycle_found to true if a cycle is found and sets

     void processNextBuildRound(int k) {

     // _cycle_length, _cycle_size, _cycle_node to represent the minimum

       NodeVector next;

     // mean cycle in the policy graph.

       for (int i = 0; i < process.size(); ++i) {

     bool findPolicyCycles(int comp) {

 	for (OutEdgeIt e(graph, process[i]); e != INVALID; ++e) {

       typename Graph::template NodeMap<int> level(_graph, -1);

 	  Node v = graph.target(e);

       _cycle_found = false;

 	  if (comp[v] != comp_cnt) continue;

       Length clength;

 	  if (!dmap[v][k].found) {

       int csize;

 	    next.push_back(v);

       int path_cnt = 0;

 	    dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e);

       Node u, v;

 	  }

       // Searching for cycles

 	  else if (dmap[process[i]][k-1].dist + length[e] < dmap[v][k].dist) {

       for (int i = 0; i < _nodes.size(); ++i) {

 	    dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e);

         if (level[_nodes[i]] < 0) {

 	  }

           u = _nodes[i];

 	}

           level[u] = path_cnt;

       }

           while (level[u = _graph.target(_policy[u])] < 0)

       process.swap(next);

             level[u] = path_cnt;

     }

           if (level[u] == path_cnt) {

             // A cycle is found

     /// \brief Processes one round of computing required path data

             clength = _length[_policy[u]];

     /// using \ref nodes vector instead of \ref process vector.

             csize = 1;

     void processNextFullRound(int k) {

             for (v = u; (v = _graph.target(_policy[v])) != u; ) {

       for (int i = 0; i < nodes.size(); ++i) {

               clength += _length[_policy[v]];

 	for (OutEdgeIt e(graph, nodes[i]); e != INVALID; ++e) {

               ++csize;

 	  Node v = graph.target(e);

             }

 	  if (comp[v] != comp_cnt) continue;

             if ( !_cycle_found ||

 	  if ( !dmap[v][k].found ||

                  clength * _cycle_size < _cycle_length * csize ) {

 	       dmap[nodes[i]][k-1].dist + length[e] < dmap[v][k].dist ) {

               _cycle_found = true;

 	    dmap[v][k] = PathData(true, dmap[nodes[i]][k-1].dist + length[e], e);

               _cycle_length = clength;

 	  }

               _cycle_size = csize;

 	}

               _cycle_node = u;

       }

             }

     }

           }

           ++path_cnt;

     /// \brief Finds the minimum cycle mean value in the current

         }

     /// component.

       }

     bool findCurrentMin(Length &min_length, int &min_size, Node &min_node) {

       return _cycle_found;

       bool found_min = false;

     }

       for (int i = 0; i < nodes.size(); ++i) {

 	int n = nodes.size();

     // Contracts the policy graph to be connected by cutting all cycles

 	if (!dmap[nodes[i]][n].found) continue;

     // except for the main cycle (i.e. the minimum mean cycle).

 	Length len;

     void contractPolicyGraph(int comp) {

 	int size;

       // Finding the component of the main cycle using

 	bool found_one = false;

       // reverse BFS search

 	for (int k = 0; k < n; ++k) {

       typename Graph::template NodeMap<int> found(_graph, false);

 	  if (!dmap[nodes[i]][k].found) continue;

       std::deque<Node> queue;

 	  Length _len = dmap[nodes[i]][n].dist - dmap[nodes[i]][k].dist;

       queue.push_back(_cycle_node);

 	  int _size = n - k;

       found[_cycle_node] = true;

 	  if (!found_one || len * _size < _len * size) {

       Node u, v;

 	    found_one = true;

       while (!queue.empty()) {

 	    len = _len;

         v = queue.front(); queue.pop_front();

 	    size = _size;

         for (InEdgeIt e(_graph, v); e != INVALID; ++e) {

 	  }

           u = _graph.source(e);

 	}

           if (_component[u] == comp && !found[u] && _policy[u] == e) {

 	if ( found_one &&

             found[u] = true;

 	     (!found_min || len * min_size < min_length * size) ) {

             queue.push_back(u);

 	  found_min = true;

           }

 	  min_length = len;

         }

 	  min_size = size;

       }

 	  min_node = nodes[i];

       // Connecting all other nodes to this component using

 	}

       // reverse BFS search

       }

       queue.clear();

       return found_min;

       for (int i = 0; i < _nodes.size(); ++i)

         if (found[_nodes[i]]) queue.push_back(_nodes[i]);

       int found_cnt = queue.size();

       while (found_cnt < _nodes.size() && !queue.empty()) {

         v = queue.front(); queue.pop_front();

         for (InEdgeIt e(_graph, v); e != INVALID; ++e) {

           u = _graph.source(e);

           if (_component[u] == comp && !found[u]) {

             found[u] = true;

             ++found_cnt;

             _policy[u] = e;

             queue.push_back(u);

           }

         }

       }

     }

     // Computes node distances in the policy graph and updates the

     // policy graph if the node distances can be improved.

     bool computeNodeDistances(int comp) {

       // Computing node distances using reverse BFS search

       double cycle_mean = (double)_cycle_length / _cycle_size;

       typename Graph::template NodeMap<int> found(_graph, 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 (InEdgeIt e(_graph, v); e != INVALID; ++e) {

           u = _graph.source(e);

           if (_component[u] == comp && !found[u] && _policy[u] == e) {

             found[u] = true;

             _dist[u] = _dist[v] + _length[e] - cycle_mean;

             queue.push_back(u);

           }

         }

       }

       // Improving node distances

       bool improved = false;

       for (int j = 0; j < _edges.size(); ++j) {

         Edge e = _edges[j];

         u = _graph.source(e); v = _graph.target(e);

         double delta = _dist[v] + _length[e] - cycle_mean;

         if (_tolerance.less(delta, _dist[u])) {

           improved = true;

           _dist[u] = delta;

           _policy[u] = e;

         }

       }

       return improved;

     /// \return \c true if a cycle exists in the graph.

     /// \return Returns \c true if a directed cycle exists in the graph.

     ///

     ///

     /// \note Apart from the return value, <tt>m.run()</tt> is just a

     /// \note Apart from the return value, <tt>mmc.run()</tt> is just a

     ///   m.init();

     ///   mmc.init();

     ///   m.findMinMean();

     ///   mmc.findMinMean();

     ///   m.findCycle();

     ///   mmc.findCycle();

       findMinMean();

       return findMinMean() && findCycle();

       return findCycle();

       comp_num = stronglyConnectedComponents(graph, comp);

       _tolerance.epsilon(1e-8);

       if (!cycle_path) {

       if (!_cycle_path) {

 	local_path = true;

         _local_path = true;

 	cycle_path = new Path;

         _cycle_path = new Path;

       }

       }

     }

       _cycle_found = false;

       _component_num = stronglyConnectedComponents(_graph, _component);

     /// \brief Finds the minimum cycle mean value in the graph.

     }

     ///

     /// Computes all the required path data and finds the minimum cycle

     /// \brief Resets the internal data structures.

     /// mean value in the graph.

     ///

     ///

     /// Resets the internal data structures so that \ref findMinMean()

     /// \return \c true if a cycle exists in the graph.

     /// and \ref findCycle() can be called again (e.g. when the

     /// underlaying graph has been modified).

     ///

     /// \sa init()

     void reset() {

       if (_cycle_path) _cycle_path->clear();

       _cycle_found = false;

       _component_num = stronglyConnectedComponents(_graph, _component);

     }

     /// \brief Finds the minimum cycle mean length in the graph.

     ///

     /// Computes all the required data and finds the minimum cycle mean

     /// length in the graph.

     ///

     /// \return Returns \c true if a directed cycle exists in the graph.

       cycle_node = INVALID;

       // Finding the minimum mean cycle in the components

       for (comp_cnt = 0; comp_cnt < comp_num; ++comp_cnt) {

       for (int comp = 0; comp < _component_num; ++comp) {

 	initCurrent();

         if (!initCurrentComponent(comp)) continue;

 	processRounds();

         while (true) {

           if (!findPolicyCycles(comp)) return false;

 	Length min_length;

           contractPolicyGraph(comp);

 	int min_size;

           if (!computeNodeDistances(comp)) return true;

 	Node min_node;

         }

 	bool found_min = findCurrentMin(min_length, min_size, min_node);

       }

     }

 	if ( found_min && (cycle_node == INVALID ||

 	     min_length * cycle_size < cycle_length * min_size) ) {

     /// \brief Finds a critical (minimum mean) directed cycle.

 	  cycle_length = min_length;

     ///

 	  cycle_size = min_size;

     /// Finds a critical (minimum mean) directed cycle using the data

 	  cycle_node = min_node;

     /// computed in the \ref findMinMean() function.

 	}

     ///

       }

     /// \return Returns \c true if a directed cycle exists in the graph.

       return (cycle_node != INVALID);

     }

     /// \brief Finds a critical (minimum mean) cycle.

     ///

     /// Finds a critical (minimum mean) cycle using the path data

     /// stored in \ref dmap.

     ///

     /// \return \c true if a cycle exists in the graph.

       if (cycle_node == INVALID) return false;

       if (!_cycle_found) return false;

       cycle_length = 0;

       _cycle_path->addBack(_policy[_cycle_node]);

       cycle_size = 0;

       for ( Node v = _cycle_node;

       IntNodeMap reached(graph, -1);

             (v = _graph.target(_policy[v])) != _cycle_node; ) {

       int r = reached[cycle_node] = dmap[cycle_node].size() - 1;

         _cycle_path->addBack(_policy[v]);

       Node u = graph.source(dmap[cycle_node][r].pred);

       while (reached[u] < 0) {

 	reached[u] = --r;

 	u = graph.source(dmap[u][r].pred);

       }

       r = reached[u];

       Edge e = dmap[u][r].pred;

       cycle_path->addFront(e);

       cycle_length = cycle_length + length[e];

       ++cycle_size;

       Node v;

       while ((v = graph.source(e)) != u) {

 	e = dmap[v][--r].pred;

 	cycle_path->addFront(e);

 	cycle_length = cycle_length + length[e];

 	++cycle_size;

     ///

     Length cycleLength() const {

     /// \warning LengthMap::Value must be convertible to double.

       return _cycle_length;

     ///

     }

     /// \note <tt>m.minMean()</tt> is just a shortcut of the following

     /// code.

     /// \brief Returns the number of edges on the found cycle.

     ///

     /// Returns the number of edges on the found cycle.

     ///

     /// \pre \ref run() or \ref findMinMean() must be called before

     /// using this function.

     int cycleEdgeNum() const {

       return _cycle_size;

     }

     /// \brief Returns the mean length of the found cycle.

     ///

     /// Returns the mean length of the found cycle.

     ///

     /// \pre \ref run() or \ref findMinMean() must be called before

     /// using this function.

     ///

     /// \note <tt>mmc.cycleMean()</tt> is just a shortcut of the

     /// following code.

     ///   return m.cycleLength() / double(m.cycleEdgeNum());

     ///   return double(mmc.cycleLength()) / mmc.cycleEdgeNum();

     /// However if only \ref findMinMean() is called before using this

     double cycleMean() const {

     /// function, <tt>m.cycleLength()</tt> and <tt>m.cycleEdgeNum()</tt>

       return (double)_cycle_length / _cycle_size;

     /// are not valid alone, only their ratio is valid.

     }

     double minMean() const {

       return cycle_length / (double)cycle_size;

     /// \brief Returns a const reference to the \ref Path "path"

     }

     /// structure storing the found cycle.

     ///

     /// \brief Returns a const reference to the \ref lemon::Path "Path"

     /// Returns a const reference to the \ref Path "path"

     /// structure of the found cycle.

     /// structure storing the found cycle.

     ///

     ///

     /// Returns a const reference to the \ref lemon::Path "Path"

     /// \pre \ref run() or \ref findCycle() must be called before using

     /// structure of the found cycle.

     /// this function.

     ///

     ///

     /// \pre \ref run() must be called before using this function.

     /// \sa cyclePath()

     ///

     /// \sa \ref cyclePath()

       return *cycle_path;

       return *_cycle_path;

     }

     }

     /// \brief Sets the \ref lemon::Path "Path" structure storing the

     /// \brief Sets the \ref Path "path" structure for storing the found

     /// cycle.

     ///

     /// Sets an external \ref Path "path" structure for storing the

     /// Sets the \ref lemon::Path "Path" structure storing the found

     /// If you don't call this function before calling \ref run() or

     /// cycle. If you don't use this function before calling

     /// \ref init(), it will allocate a local \ref Path "path"

     /// \ref run(), it will allocate one. The destuctor deallocates

     /// this automatically allocated map, of course.

     ///

     /// \note The algorithm calls only the \ref lemon::Path::addFront()

     /// "addFront()" method of the given \ref lemon::Path "Path"

     ///

     /// The destuctor deallocates this automatically allocated map,

     /// \return \c (*this)

     /// of course.

     ///

     /// \note The algorithm calls only the \ref lemon::Path::addBack()

     /// "addBack()" function of the given \ref Path "path" structure.

     ///

     /// \return <tt>(*this)</tt>

     ///

     /// \sa cycle()

       if (local_path) {

       if (_local_path) {

 	delete cycle_path;

         delete _cycle_path;

 	local_path = false;

         _local_path = false;

       }

       }

       cycle_path = &path;

       _cycle_path = &path;