alpar@2409: /* -*- C++ -*-
alpar@2409:  *
alpar@2409:  * This file is a part of LEMON, a generic C++ optimization library
alpar@2409:  *
alpar@2409:  * Copyright (C) 2003-2007
alpar@2409:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@2409:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@2409:  *
alpar@2409:  * Permission to use, modify and distribute this software is granted
alpar@2409:  * provided that this copyright notice appears in all copies. For
alpar@2409:  * precise terms see the accompanying LICENSE file.
alpar@2409:  *
alpar@2409:  * This software is provided "AS IS" with no warranty of any kind,
alpar@2409:  * express or implied, and with no claim as to its suitability for any
alpar@2409:  * purpose.
alpar@2409:  *
alpar@2409:  */
alpar@2409: 
alpar@2409: #ifndef LEMON_MIN_MEAN_CYCLE_H
alpar@2409: #define LEMON_MIN_MEAN_CYCLE_H
alpar@2409: 
deba@2529: /// \ingroup shortest_path
alpar@2409: ///
deba@2437: /// \file
kpeter@2509: /// \brief Karp's algorithm for finding a minimum mean (directed) cycle.
alpar@2409: 
kpeter@2509: #include <vector>
alpar@2409: #include <lemon/graph_utils.h>
alpar@2409: #include <lemon/topology.h>
alpar@2409: #include <lemon/path.h>
alpar@2409: 
alpar@2409: namespace lemon {
alpar@2409: 
deba@2529:   /// \addtogroup shortest_path
alpar@2409:   /// @{
alpar@2409: 
deba@2413:   /// \brief Implementation of Karp's algorithm for finding a
deba@2413:   /// minimum mean (directed) cycle.
alpar@2409:   ///
deba@2413:   /// The \ref lemon::MinMeanCycle "MinMeanCycle" implements Karp's
deba@2413:   /// algorithm for finding a minimum mean (directed) cycle.
alpar@2409:   ///
alpar@2409:   /// \param Graph The directed graph type the algorithm runs on.
alpar@2409:   /// \param LengthMap The type of the length (cost) map.
alpar@2409:   ///
alpar@2409:   /// \author Peter Kovacs
alpar@2409: 
alpar@2409:   template <typename Graph,
alpar@2409:     typename LengthMap = typename Graph::template EdgeMap<int> >
alpar@2409:   class MinMeanCycle
alpar@2409:   {
alpar@2409:     typedef typename Graph::Node Node;
alpar@2409:     typedef typename Graph::NodeIt NodeIt;
alpar@2409:     typedef typename Graph::Edge Edge;
alpar@2409:     typedef typename Graph::EdgeIt EdgeIt;
alpar@2409:     typedef typename Graph::OutEdgeIt OutEdgeIt;
alpar@2409: 
alpar@2409:     typedef typename LengthMap::Value Length;
deba@2413: 
alpar@2409:     typedef typename Graph::template NodeMap<int> IntNodeMap;
alpar@2409:     typedef typename Graph::template NodeMap<Edge> PredNodeMap;
alpar@2409:     typedef Path<Graph> Path;
alpar@2409:     typedef std::vector<Node> NodeVector;
alpar@2409: 
alpar@2409:   protected:
deba@2413: 
kpeter@2509:     /// \brief Data structure for path data.
deba@2413:     struct PathData
deba@2413:     {
alpar@2409:       bool found;
alpar@2409:       Length dist;
alpar@2409:       Edge pred;
deba@2413:       PathData(bool _found = false, Length _dist = 0) :
alpar@2409: 	found(_found), dist(_dist), pred(INVALID) {}
deba@2413:       PathData(bool _found, Length _dist, Edge _pred) :
alpar@2409: 	found(_found), dist(_dist), pred(_pred) {}
alpar@2409:     };
deba@2413: 
alpar@2409:   private:
deba@2413: 
deba@2413:     typedef typename Graph::template NodeMap<std::vector<PathData> >
deba@2413:       PathDataNodeMap;
deba@2413: 
alpar@2409:   protected:
alpar@2409: 
alpar@2409:     /// \brief Node map for storing path data.
deba@2437:     ///
alpar@2409:     /// Node map for storing path data of all nodes in the current
deba@2437:     /// component. dmap[v][k] is the length of a shortest directed walk
alpar@2409:     /// to node v from the starting node containing exactly k edges.
deba@2413:     PathDataNodeMap dmap;
deba@2437: 
alpar@2409:     /// \brief The directed graph the algorithm runs on.
deba@2437:     const Graph &graph;
alpar@2409:     /// \brief The length of the edges.
deba@2437:     const LengthMap &length;
deba@2437: 
alpar@2409:     /// \brief The total length of the found cycle.
alpar@2409:     Length cycle_length;
alpar@2409:     /// \brief The number of edges in the found cycle.
alpar@2409:     int cycle_size;
deba@2437:     /// \brief A node for obtaining a minimum mean cycle.
deba@2413:     Node cycle_node;
deba@2413: 
alpar@2409:     /// \brief The found cycle.
alpar@2409:     Path *cycle_path;
deba@2437:     /// \brief The algorithm uses local \ref lemon::Path "Path"
deba@2413:     /// structure to store the found cycle.
alpar@2409:     bool local_path;
deba@2437: 
alpar@2409:     /// \brief Node map for identifying strongly connected components.
alpar@2409:     IntNodeMap comp;
alpar@2409:     /// \brief The number of strongly connected components.
alpar@2409:     int comp_num;
deba@2413:     /// \brief Counter for identifying the current component.
alpar@2409:     int comp_cnt;
alpar@2409:     /// \brief Nodes of the current component.
alpar@2409:     NodeVector nodes;
alpar@2409:     /// \brief The processed nodes in the last round.
alpar@2409:     NodeVector process;
deba@2437: 
alpar@2409:   public :
alpar@2409: 
alpar@2409:     /// \brief The constructor of the class.
alpar@2409:     ///
alpar@2409:     /// The constructor of the class.
alpar@2409:     ///
alpar@2409:     /// \param _graph The directed graph the algorithm runs on.
alpar@2409:     /// \param _length The length (cost) of the edges.
deba@2437:     MinMeanCycle( const Graph &_graph,
deba@2437: 		  const LengthMap &_length ) :
deba@2413:       graph(_graph), length(_length), dmap(_graph), comp(_graph),
deba@2413:       cycle_length(0), cycle_size(-1), cycle_node(INVALID),
alpar@2409:       cycle_path(NULL), local_path(false)
alpar@2409:     { }
alpar@2409: 
alpar@2409:     /// \brief The destructor of the class.
deba@2437:     ~MinMeanCycle() {
alpar@2409:       if (local_path) delete cycle_path;
alpar@2409:     }
alpar@2409: 
alpar@2409:   protected:
alpar@2409: 
alpar@2409:     /// \brief Initializes the internal data structures for the current
alpar@2409:     /// component.
alpar@2409:     void initCurrent() {
alpar@2409:       nodes.clear();
alpar@2409:       // Finding the nodes of the current component
alpar@2409:       for (NodeIt v(graph); v != INVALID; ++v) {
alpar@2409: 	if (comp[v] == comp_cnt) nodes.push_back(v);
alpar@2409:       }
alpar@2409:       // Creating vectors for all nodes
alpar@2409:       int n = nodes.size();
kpeter@2509:       for (int i = 0; i < n; ++i) {
deba@2437: 	dmap[nodes[i]].resize(n + 1);
alpar@2409:       }
alpar@2409:     }
deba@2437: 
deba@2437:     /// \brief Processes all rounds of computing required path data for
deba@2413:     /// the current component.
alpar@2409:     void processRounds() {
alpar@2409:       dmap[nodes[0]][0] = PathData(true, 0);
alpar@2409:       process.clear();
deba@2413:       // Processing the first round
deba@2413:       for (OutEdgeIt e(graph, nodes[0]); e != INVALID; ++e) {
deba@2413: 	Node v = graph.target(e);
alpar@2409: 	if (comp[v] != comp_cnt || v == nodes[0]) continue;
deba@2413: 	dmap[v][1] = PathData(true, length[e], e);
alpar@2409: 	process.push_back(v);
alpar@2409:       }
deba@2437:       // Processing other rounds
alpar@2409:       int n = nodes.size(), k;
deba@2437:       for (k = 2; k <= n && process.size() < n; ++k)
alpar@2409: 	processNextBuildRound(k);
deba@2437:       for ( ; k <= n; ++k)
alpar@2409: 	processNextFullRound(k);
alpar@2409:     }
deba@2437: 
alpar@2409:     /// \brief Processes one round of computing required path data and
alpar@2409:     /// rebuilds \ref process vector.
alpar@2409:     void processNextBuildRound(int k) {
alpar@2409:       NodeVector next;
deba@2437:       for (int i = 0; i < process.size(); ++i) {
deba@2437: 	for (OutEdgeIt e(graph, process[i]); e != INVALID; ++e) {
deba@2413: 	  Node v = graph.target(e);
alpar@2409: 	  if (comp[v] != comp_cnt) continue;
deba@2413: 	  if (!dmap[v][k].found) {
deba@2413: 	    next.push_back(v);
deba@2437: 	    dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e);
deba@2413: 	  }
deba@2437: 	  else if (dmap[process[i]][k-1].dist + length[e] < dmap[v][k].dist) {
deba@2437: 	    dmap[v][k] = PathData(true, dmap[process[i]][k-1].dist + length[e], e);
alpar@2409: 	  }
alpar@2409: 	}
alpar@2409:       }
alpar@2409:       process.swap(next);
alpar@2409:     }
alpar@2409: 
alpar@2409:     /// \brief Processes one round of computing required path data
alpar@2409:     /// using \ref nodes vector instead of \ref process vector.
alpar@2409:     void processNextFullRound(int k) {
deba@2437:       for (int i = 0; i < nodes.size(); ++i) {
deba@2437: 	for (OutEdgeIt e(graph, nodes[i]); e != INVALID; ++e) {
deba@2413: 	  Node v = graph.target(e);
alpar@2409: 	  if (comp[v] != comp_cnt) continue;
deba@2437: 	  if ( !dmap[v][k].found ||
deba@2437: 	       dmap[nodes[i]][k-1].dist + length[e] < dmap[v][k].dist ) {
deba@2437: 	    dmap[v][k] = PathData(true, dmap[nodes[i]][k-1].dist + length[e], e);
alpar@2409: 	  }
alpar@2409: 	}
alpar@2409:       }
alpar@2409:     }
deba@2437: 
deba@2437:     /// \brief Finds the minimum cycle mean value in the current
deba@2413:     /// component.
deba@2413:     bool findCurrentMin(Length &min_length, int &min_size, Node &min_node) {
alpar@2409:       bool found_min = false;
deba@2437:       for (int i = 0; i < nodes.size(); ++i) {
deba@2413: 	int n = nodes.size();
deba@2437: 	if (!dmap[nodes[i]][n].found) continue;
alpar@2409: 	Length len;
alpar@2409: 	int size;
alpar@2409: 	bool found_one = false;
alpar@2409: 	for (int k = 0; k < n; ++k) {
deba@2437: 	  if (!dmap[nodes[i]][k].found) continue;
deba@2437: 	  Length _len = dmap[nodes[i]][n].dist - dmap[nodes[i]][k].dist;
deba@2437: 	  int _size = n - k;
deba@2413: 	  if (!found_one || len * _size < _len * size) {
alpar@2409: 	    found_one = true;
alpar@2409: 	    len = _len;
alpar@2409: 	    size = _size;
alpar@2409: 	  }
alpar@2409: 	}
deba@2437: 	if ( found_one &&
deba@2413: 	     (!found_min || len * min_size < min_length * size) ) {
alpar@2409: 	  found_min = true;
alpar@2409: 	  min_length = len;
alpar@2409: 	  min_size = size;
deba@2437: 	  min_node = nodes[i];
alpar@2409: 	}
alpar@2409:       }
alpar@2409:       return found_min;
alpar@2409:     }
deba@2437: 
deba@2437:   public:
deba@2437: 
alpar@2409:     /// \brief Runs the algorithm.
alpar@2409:     ///
alpar@2409:     /// Runs the algorithm.
alpar@2409:     ///
alpar@2409:     /// \return \c true if a cycle exists in the graph.
deba@2413:     ///
kpeter@2517:     /// \note Apart from the return value, <tt>m.run()</tt> is just a
kpeter@2517:     /// shortcut of the following code.
deba@2413:     /// \code
deba@2413:     ///   m.init();
deba@2413:     ///   m.findMinMean();
deba@2413:     ///   m.findCycle();
deba@2413:     /// \endcode
alpar@2409:     bool run() {
alpar@2409:       init();
deba@2413:       findMinMean();
deba@2413:       return findCycle();
deba@2413:     }
deba@2437: 
deba@2413:     /// \brief Initializes the internal data structures.
kpeter@2517:     ///
kpeter@2517:     /// Initializes the internal data structures.
kpeter@2517:     ///
kpeter@2517:     /// \sa reset()
deba@2413:     void init() {
deba@2413:       comp_num = stronglyConnectedComponents(graph, comp);
deba@2413:       if (!cycle_path) {
deba@2413: 	local_path = true;
deba@2413: 	cycle_path = new Path;
deba@2413:       }
deba@2413:     }
deba@2413: 
deba@2413:     /// \brief Finds the minimum cycle mean value in the graph.
deba@2413:     ///
deba@2413:     /// Computes all the required path data and finds the minimum cycle
deba@2413:     /// mean value in the graph.
deba@2413:     ///
deba@2413:     /// \return \c true if a cycle exists in the graph.
deba@2437:     ///
deba@2413:     /// \pre \ref init() must be called before using this function.
deba@2413:     bool findMinMean() {
deba@2413:       cycle_node = INVALID;
alpar@2409:       for (comp_cnt = 0; comp_cnt < comp_num; ++comp_cnt) {
alpar@2409: 	initCurrent();
alpar@2409: 	processRounds();
deba@2437: 
alpar@2409: 	Length min_length;
alpar@2409: 	int min_size;
alpar@2409: 	Node min_node;
deba@2413: 	bool found_min = findCurrentMin(min_length, min_size, min_node);
deba@2437: 
deba@2437: 	if ( found_min && (cycle_node == INVALID ||
deba@2413: 	     min_length * cycle_size < cycle_length * min_size) ) {
alpar@2409: 	  cycle_length = min_length;
alpar@2409: 	  cycle_size = min_size;
alpar@2409: 	  cycle_node = min_node;
alpar@2409: 	}
alpar@2409:       }
deba@2413:       return (cycle_node != INVALID);
alpar@2409:     }
deba@2437: 
deba@2413:     /// \brief Finds a critical (minimum mean) cycle.
alpar@2409:     ///
deba@2413:     /// Finds a critical (minimum mean) cycle using the path data
deba@2413:     /// stored in \ref dmap.
deba@2413:     ///
deba@2413:     /// \return \c true if a cycle exists in the graph.
deba@2437:     ///
deba@2437:     /// \pre \ref init() and \ref findMinMean() must be called before
deba@2413:     /// using this function.
deba@2413:     bool findCycle() {
deba@2413:       if (cycle_node == INVALID) return false;
deba@2413:       cycle_length = 0;
deba@2413:       cycle_size = 0;
deba@2413:       IntNodeMap reached(graph, -1);
deba@2413:       int r = reached[cycle_node] = dmap[cycle_node].size() - 1;
deba@2413:       Node u = graph.source(dmap[cycle_node][r].pred);
deba@2413:       while (reached[u] < 0) {
deba@2413: 	reached[u] = --r;
deba@2413: 	u = graph.source(dmap[u][r].pred);
deba@2413:       }
deba@2413:       r = reached[u];
deba@2413:       Edge e = dmap[u][r].pred;
deba@2413:       cycle_path->addFront(e);
deba@2413:       cycle_length = cycle_length + length[e];
deba@2413:       ++cycle_size;
deba@2413:       Node v;
deba@2413:       while ((v = graph.source(e)) != u) {
deba@2413: 	e = dmap[v][--r].pred;
deba@2413: 	cycle_path->addFront(e);
deba@2413: 	cycle_length = cycle_length + length[e];
deba@2413: 	++cycle_size;
deba@2413:       }
deba@2413:       return true;
deba@2413:     }
deba@2413: 
deba@2413:     /// \brief Resets the internal data structures.
deba@2413:     ///
deba@2437:     /// Resets the internal data structures so that \ref findMinMean()
deba@2437:     /// and \ref findCycle() can be called again (e.g. when the
deba@2413:     /// underlaying graph has been modified).
kpeter@2517:     ///
kpeter@2517:     /// \sa init()
deba@2413:     void reset() {
deba@2413:       for (NodeIt u(graph); u != INVALID; ++u)
deba@2413: 	dmap[u].clear();
deba@2413:       cycle_node = INVALID;
deba@2413:       if (cycle_path) cycle_path->clear();
deba@2413:       comp_num = stronglyConnectedComponents(graph, comp);
deba@2413:     }
deba@2437: 
deba@2413:     /// \brief Returns the total length of the found cycle.
deba@2413:     ///
deba@2413:     /// Returns the total length of the found cycle.
alpar@2409:     ///
kpeter@2517:     /// \pre \ref run() or \ref findCycle() must be called before
kpeter@2517:     /// using this function. If only \ref findMinMean() is called,
kpeter@2517:     /// the return value is not valid.
deba@2437:     Length cycleLength() const {
alpar@2409:       return cycle_length;
alpar@2409:     }
deba@2437: 
deba@2413:     /// \brief Returns the number of edges in the found cycle.
alpar@2409:     ///
deba@2413:     /// Returns the number of edges in the found cycle.
alpar@2409:     ///
kpeter@2517:     /// \pre \ref run() or \ref findCycle() must be called before
kpeter@2517:     /// using this function. If only \ref findMinMean() is called,
kpeter@2517:     /// the return value is not valid.
deba@2437:     int cycleEdgeNum() const {
alpar@2409:       return cycle_size;
alpar@2409:     }
deba@2437: 
deba@2413:     /// \brief Returns the mean length of the found cycle.
alpar@2409:     ///
deba@2413:     /// Returns the mean length of the found cycle.
alpar@2409:     ///
kpeter@2517:     /// \pre \ref run() or \ref findMinMean() must be called before
kpeter@2517:     /// using this function.
alpar@2409:     ///
alpar@2409:     /// \warning LengthMap::Value must be convertible to double.
deba@2413:     ///
kpeter@2517:     /// \note <tt>m.minMean()</tt> is just a shortcut of the following
kpeter@2517:     /// code.
deba@2413:     /// \code
kpeter@2517:     ///   return m.cycleLength() / double(m.cycleEdgeNum());
deba@2413:     /// \endcode
kpeter@2517:     /// However if only \ref findMinMean() is called before using this
kpeter@2517:     /// function, <tt>m.cycleLength()</tt> and <tt>m.cycleEdgeNum()</tt>
kpeter@2517:     /// are not valid alone, only their ratio is valid.
deba@2437:     double minMean() const {
deba@2413:       return cycle_length / (double)cycle_size;
alpar@2409:     }
alpar@2409: 
alpar@2409:     /// \brief Returns a const reference to the \ref lemon::Path "Path"
deba@2413:     /// structure of the found cycle.
alpar@2409:     ///
alpar@2409:     /// Returns a const reference to the \ref lemon::Path "Path"
deba@2413:     /// structure of the found cycle.
alpar@2409:     ///
alpar@2409:     /// \pre \ref run() must be called before using this function.
alpar@2409:     ///
alpar@2409:     /// \sa \ref cyclePath()
deba@2437:     const Path& cycle() const {
alpar@2409:       return *cycle_path;
alpar@2409:     }
deba@2437: 
deba@2437:     /// \brief Sets the \ref lemon::Path "Path" structure storing the
alpar@2409:     /// found cycle.
deba@2437:     ///
deba@2437:     /// Sets the \ref lemon::Path "Path" structure storing the found
deba@2437:     /// cycle. If you don't use this function before calling
deba@2437:     /// \ref run(), it will allocate one. The destuctor deallocates
alpar@2409:     /// this automatically allocated map, of course.
alpar@2409:     ///
alpar@2409:     /// \note The algorithm calls only the \ref lemon::Path::addFront()
deba@2437:     /// "addFront()" method of the given \ref lemon::Path "Path"
alpar@2409:     /// structure.
deba@2437:     ///
alpar@2409:     /// \return \c (*this)
deba@2437:     MinMeanCycle& cyclePath(Path &path) {
alpar@2409:       if (local_path) {
alpar@2409: 	delete cycle_path;
alpar@2409: 	local_path = false;
alpar@2409:       }
alpar@2409:       cycle_path = &path;
alpar@2409:       return *this;
alpar@2409:     }
alpar@2409: 
alpar@2409:   }; //class MinMeanCycle
alpar@2409: 
alpar@2409:   ///@}
alpar@2409: 
alpar@2409: } //namespace lemon
alpar@2409: 
alpar@2409: #endif //LEMON_MIN_MEAN_CYCLE_H