/* -*- 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