/* -*- C++ -*-

  *

  * This file is a part of LEMON, a generic C++ optimization library

  *

  * Copyright (C) 2003-2007

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  *

  * Permission to use, modify and distribute this software is granted

  * provided that this copyright notice appears in all copies. For

  * precise terms see the accompanying LICENSE file.

  *

  * This software is provided "AS IS" with no warranty of any kind,

  * express or implied, and with no claim as to its suitability for any

  * purpose.

  *

  */

 #ifndef LEMON_MIN_MEAN_CYCLE_H

 #define LEMON_MIN_MEAN_CYCLE_H

 /// \ingroup min_cost_flow

 ///

 /// \file 

 /// \brief Karp algorithm for finding a minimum mean cycle.

 #include <lemon/graph_utils.h>

 #include <lemon/topology.h>

 #include <lemon/path.h>

 namespace lemon {

   /// \addtogroup min_cost_flow

   /// @{

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

   /// minimum mean cycle.

   ///

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

   /// algorithm for finding a minimum mean cycle.

   ///

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

   /// \param LengthMap The type of the length (cost) map.

   ///

   /// \author Peter Kovacs

 #ifdef DOXYGEN

   template <typename Graph, typename LengthMap>

 #else

   template <typename Graph,

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

 #endif

   class MinMeanCycle

   {

     typedef typename Graph::Node Node;

     typedef typename Graph::NodeIt NodeIt;

     typedef typename Graph::Edge Edge;

     typedef typename Graph::EdgeIt EdgeIt;

     typedef typename Graph::InEdgeIt InEdgeIt;

     typedef typename Graph::OutEdgeIt OutEdgeIt;

     typedef typename LengthMap::Value Length;

     typedef typename Graph::template NodeMap<int> IntNodeMap;

     typedef typename Graph::template NodeMap<Edge> PredNodeMap;

     typedef Path<Graph> Path;

     typedef std::vector<Node> NodeVector;

     typedef typename NodeVector::iterator NodeVectorIt;

   protected:

     /// \brief Data sturcture for path data.

     struct PathData {

       bool found;

       Length dist;

       Edge pred;

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

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

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

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

     };

   private:

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

   protected:

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

     /// 

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

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

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

     PathVectorNodeMap 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 The found cycle.

     Path *cycle_path;

     /// \brief 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), 

       cycle_length(0), cycle_size(-1), comp(_graph),

       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.

     void init() {

       comp_num = stronglyConnectedComponents(graph, comp);

       if (!cycle_path) {

 	local_path = true;

 	cycle_path = new Path;

       }

     }

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

     /// component.

     void initCurrent() {

       nodes.clear();

       // Finding the nodes of the current component

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

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

       }

       // Creating vectors for all nodes

       int n = nodes.size();

       for (NodeVectorIt vi = nodes.begin(); vi != nodes.end(); ++vi) {

 	dmap[*vi].resize(n + 1);

       }

     }

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

     void processRounds() {

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

       process.clear();

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

 	Node v = graph.target(oe);

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

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

 	process.push_back(v);

       }

       int n = nodes.size(), k;

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

 	processNextBuildRound(k);

       }

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

 	processNextFullRound(k);

       }

     }

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

     /// rebuilds \ref process vector.

     void processNextBuildRound(int k) {

       NodeVector next;

       for (NodeVectorIt ui = process.begin(); ui != process.end(); ++ui) {

 	for (OutEdgeIt oe(graph, *ui); oe != INVALID; ++oe) {

 	  Node v = graph.target(oe);

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

 	  if ( !dmap[v][k].found || (dmap[v][k].dist > dmap[*ui][k-1].dist + length[oe]) ) {

 	    if (!dmap[v][k].found) next.push_back(v); 

 	    dmap[v][k] = PathData(true, dmap[*ui][k-1].dist + length[oe], oe);

 	  }

 	}

       }

       process.swap(next);

     }

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

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

     void processNextFullRound(int k) {

       for (NodeVectorIt ui = nodes.begin(); ui != nodes.end(); ++ui) {

 	for (OutEdgeIt oe(graph, *ui); oe != INVALID; ++oe) {

 	  Node v = graph.target(oe);

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

 	  if ( !dmap[v][k].found || (dmap[v][k].dist > dmap[*ui][k-1].dist + length[oe]) ) {

 	    dmap[v][k] = PathData(true, dmap[*ui][k-1].dist + length[oe], oe);

 	  }

 	}

       }

     }

     /// \brief Finds the minimum cycle mean value according to the path

     /// data stored in \ref dmap.

     ///

     /// Finds the minimum cycle mean value according to the path data

     /// stored in \ref dmap.

     ///

     /// \retval min_length The total length of the found cycle.

     /// \retval min_size The number of edges in the found cycle.

     /// \retval min_node A node for obtaining a critical cycle.

     ///

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

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

       bool found_min = false;

       for (NodeVectorIt vi = nodes.begin(); vi != nodes.end(); ++vi) {

 	Length len;

 	int size;

 	bool found_one = false;

 	int n = nodes.size();

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

 	  Length _len = dmap[*vi][n].dist - dmap[*vi][k].dist;

 	  int _size = n - k; 

 	  if ( dmap[*vi][n].found && dmap[*vi][k].found && (!found_one || len * _size < _len * size) ) {

 	    found_one = true;

 	    len = _len;

 	    size = _size;

 	  }

 	}

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

 	  found_min = true;

 	  min_length = len;

 	  min_size = size;

 	  min_node = *vi;

 	}

       }

       return found_min;

     }

     /// \brief Finds a critical (minimum mean) cycle according to the

     /// path data stored in \ref dmap.

     void findCycle(const Node &min_n) {

       IntNodeMap reached(graph, -1);

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

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

       while (reached[u] < 0) {

 	reached[u] = --r;

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

       }

       r = reached[u];

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

       cycle_path->addFront(ce);

       Node v;

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

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

 	cycle_path->addFront(ce);

       }

     }

   public:

     /// \brief Runs the algorithm.

     ///

     /// Runs the algorithm.

     ///

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

     bool run() {

       init();

       // Searching for the minimum mean cycle in all components

       bool found_cycle = false;

       Node cycle_node;

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

 	initCurrent();

 	processRounds();

 	Length min_length;

 	int min_size;

 	Node min_node;

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

 	if ( found_min && (!found_cycle || min_length * cycle_size  < cycle_length * min_size) ) {

 	  found_cycle = true;

 	  cycle_length = min_length;

 	  cycle_size = min_size;

 	  cycle_node = min_node;

 	}

       }

       if (found_cycle) findCycle(cycle_node);

       return found_cycle;

     }

     /// \brief Returns the total length of the found critical cycle.

     ///

     /// Returns the total length of the found critical cycle.

     ///

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

     Length cycleLength() const { 

       return cycle_length;

     }

     /// \brief Returns the number of edges in the found critical 

     /// cycle.

     ///

     /// Returns the number of edges in the found critical cycle.

     ///

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

     int cycleEdgeNum() const { 

       return cycle_size;

     }

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

     ///

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

     ///

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

     ///

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

     double minMean() const { 

       return (double)cycle_length / (double)cycle_size;

     }

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

     /// structure of the found flow.

     ///

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

     /// structure of the found flow.

     ///

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

     ///

     /// \sa \ref cyclePath()

     const Path &cycle() const {

       return *cycle_path;

     }

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

     /// found cycle.

     /// 

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

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

     /// \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" 

     /// structure.

     /// 

     /// \return \c (*this)

     MinMeanCycle &cyclePath(Path& path) {

       if (local_path) {

 	delete cycle_path;

 	local_path = false;

       }

       cycle_path = &path;

       return *this;

     }

   }; //class MinMeanCycle

   ///@}

 } //namespace lemon

 #endif //LEMON_MIN_MEAN_CYCLE_H