/* -*- 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's algorithm for finding a minimum mean (directed) cycle.

#include <vector>
#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 (directed) cycle.
  ///
  /// The \ref lemon::MinMeanCycle "MinMeanCycle" implements Karp's
  /// algorithm for finding a minimum mean (directed) cycle.
  ///
  /// \param Graph The directed graph type the algorithm runs on.
  /// \param LengthMap The type of the length (cost) map.
  ///
  /// \author Peter Kovacs

  template <typename Graph,
    typename LengthMap = typename Graph::template EdgeMap<int> >
  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::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;

  protected:

    /// \brief Data structure 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> >
      PathDataNodeMap;

  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.
    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();
      // 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 (int i = 0; i < n; ++i) {
	dmap[nodes[i]].resize(n + 1);
      }
    }

    /// \brief Processes all rounds of computing required path data for
    /// the current component.
    void processRounds() {
      dmap[nodes[0]][0] = PathData(true, 0);
      process.clear();
      // Processing the first round
      for (OutEdgeIt e(graph, nodes[0]); e != INVALID; ++e) {
	Node v = graph.target(e);
	if (comp[v] != comp_cnt || v == nodes[0]) continue;
	dmap[v][1] = PathData(true, length[e], e);
	process.push_back(v);
      }
      // Processing other rounds
      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 (int i = 0; i < process.size(); ++i) {
	for (OutEdgeIt e(graph, process[i]); e != INVALID; ++e) {
	  Node v = graph.target(e);
	  if (comp[v] != comp_cnt) continue;
	  if (!dmap[v][k].found) {
	    next.push_back(v);
	    dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e);
	  }
	  else if (dmap[process[i]][k-1].dist + length[e] < dmap[v][k].dist) {
	    dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e);
	  }
	}
      }
      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 (int i = 0; i < nodes.size(); ++i) {
	for (OutEdgeIt e(graph, nodes[i]); e != INVALID; ++e) {
	  Node v = graph.target(e);
	  if (comp[v] != comp_cnt) continue;
	  if ( !dmap[v][k].found ||
	       dmap[nodes[i]][k-1].dist + length[e] < dmap[v][k].dist ) {
	    dmap[v][k] = PathData(true, dmap[nodes[i]][k-1].dist + length[e], e);
	  }
	}
      }
    }

    /// \brief Finds the minimum cycle mean value in the current
    /// component.
    bool findCurrentMin(Length &min_length, int &min_size, Node &min_node) {
      bool found_min = false;
      for (int i = 0; i < nodes.size(); ++i) {
	int n = nodes.size();
	if (!dmap[nodes[i]][n].found) continue;
	Length len;
	int size;
	bool found_one = false;
	for (int k = 0; k < n; ++k) {
	  if (!dmap[nodes[i]][k].found) continue;
	  Length _len = dmap[nodes[i]][n].dist - dmap[nodes[i]][k].dist;
	  int _size = n - k;
	  if (!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 = nodes[i];
	}
      }
      return found_min;
    }

  public:

    /// \brief Runs the algorithm.
    ///
    /// Runs the algorithm.
    ///
    /// \return \c true if a cycle exists in the graph.
    ///
    /// \note Apart from the return value, <tt>m.run()</tt> is just a
    /// shortcut of the following code.
    /// \code
    ///   m.init();
    ///   m.findMinMean();
    ///   m.findCycle();
    /// \endcode
    bool run() {
      init();
      findMinMean();
      return findCycle();
    }

    /// \brief Initializes the internal data structures.
    ///
    /// Initializes the internal data structures.
    ///
    /// \sa reset()
    void init() {
      comp_num = stronglyConnectedComponents(graph, comp);
      if (!cycle_path) {
	local_path = true;
	cycle_path = new Path;
      }
    }

    /// \brief Finds the minimum cycle mean value in the graph.
    ///
    /// Computes all the required path data and finds the minimum cycle
    /// mean value in the graph.
    ///
    /// \return \c true if a cycle exists in the graph.
    ///
    /// \pre \ref init() must be called before using this function.
    bool findMinMean() {
      cycle_node = INVALID;
      for (comp_cnt = 0; comp_cnt < comp_num; ++comp_cnt) {
	initCurrent();
	processRounds();

	Length min_length;
	int min_size;
	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) ) {
	  cycle_length = min_length;
	  cycle_size = min_size;
	  cycle_node = min_node;
	}
      }
      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.
    ///
    /// \pre \ref init() and \ref findMinMean() must be called before
    /// using this function.
    bool findCycle() {
      if (cycle_node == INVALID) return false;
      cycle_length = 0;
      cycle_size = 0;
      IntNodeMap reached(graph, -1);
      int r = reached[cycle_node] = dmap[cycle_node].size() - 1;
      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;
      }
      return true;
    }

    /// \brief Resets the internal data structures.
    ///
    /// Resets the internal data structures so that \ref findMinMean()
    /// and \ref findCycle() can be called again (e.g. when the
    /// underlaying graph has been modified).
    ///
    /// \sa init()
    void reset() {
      for (NodeIt u(graph); u != INVALID; ++u)
	dmap[u].clear();
      cycle_node = INVALID;
      if (cycle_path) cycle_path->clear();
      comp_num = stronglyConnectedComponents(graph, comp);
    }

    /// \brief Returns the total length of the found cycle.
    ///
    /// Returns the total length of the found cycle.
    ///
    /// \pre \ref run() or \ref findCycle() must be called before
    /// using this function. If only \ref findMinMean() is called,
    /// the return value is not valid.
    Length cycleLength() const {
      return cycle_length;
    }

    /// \brief Returns the number of edges in the found cycle.
    ///
    /// Returns the number of edges in the found cycle.
    ///
    /// \pre \ref run() or \ref findCycle() must be called before
    /// using this function. If only \ref findMinMean() is called,
    /// the return value is not valid.
    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.
    ///
    /// \warning LengthMap::Value must be convertible to double.
    ///
    /// \note <tt>m.minMean()</tt> is just a shortcut of the following
    /// code.
    /// \code
    ///   return m.cycleLength() / double(m.cycleEdgeNum());
    /// \endcode
    /// However if only \ref findMinMean() is called before using this
    /// function, <tt>m.cycleLength()</tt> and <tt>m.cycleEdgeNum()</tt>
    /// 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 lemon::Path "Path"
    /// structure of the found cycle.
    ///
    /// Returns a const reference to the \ref lemon::Path "Path"
    /// structure of the found cycle.
    ///
    /// \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