1.1 --- a/lemon/Makefile.am Tue Mar 13 15:42:06 2007 +0000
1.2 +++ b/lemon/Makefile.am Tue Mar 13 16:32:35 2007 +0000
1.3 @@ -89,6 +89,7 @@
1.4 lemon/matrix_maps.h \
1.5 lemon/max_matching.h \
1.6 lemon/min_cost_arborescence.h \
1.7 + lemon/min_mean_cycle.h \
1.8 lemon/mip_glpk.h \
1.9 lemon/mip_cplex.h \
1.10 lemon/nagamochi_ibaraki.h \
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/min_mean_cycle.h Tue Mar 13 16:32:35 2007 +0000
2.3 @@ -0,0 +1,380 @@
2.4 +/* -*- C++ -*-
2.5 + *
2.6 + * This file is a part of LEMON, a generic C++ optimization library
2.7 + *
2.8 + * Copyright (C) 2003-2007
2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
2.11 + *
2.12 + * Permission to use, modify and distribute this software is granted
2.13 + * provided that this copyright notice appears in all copies. For
2.14 + * precise terms see the accompanying LICENSE file.
2.15 + *
2.16 + * This software is provided "AS IS" with no warranty of any kind,
2.17 + * express or implied, and with no claim as to its suitability for any
2.18 + * purpose.
2.19 + *
2.20 + */
2.21 +
2.22 +#ifndef LEMON_MIN_MEAN_CYCLE_H
2.23 +#define LEMON_MIN_MEAN_CYCLE_H
2.24 +
2.25 +/// \ingroup min_cost_flow
2.26 +///
2.27 +/// \file
2.28 +/// \brief Karp algorithm for finding a minimum mean cycle.
2.29 +
2.30 +#include <lemon/graph_utils.h>
2.31 +#include <lemon/topology.h>
2.32 +#include <lemon/path.h>
2.33 +
2.34 +namespace lemon {
2.35 +
2.36 + /// \addtogroup min_cost_flow
2.37 + /// @{
2.38 +
2.39 + /// \brief Implementation of Karp's algorithm for finding a
2.40 + /// minimum mean cycle.
2.41 + ///
2.42 + /// The \ref lemon::MinMeanCycle "MinMeanCycle" implements Karp's
2.43 + /// algorithm for finding a minimum mean cycle.
2.44 + ///
2.45 + /// \param Graph The directed graph type the algorithm runs on.
2.46 + /// \param LengthMap The type of the length (cost) map.
2.47 + ///
2.48 + /// \author Peter Kovacs
2.49 +
2.50 +#ifdef DOXYGEN
2.51 + template <typename Graph, typename LengthMap>
2.52 +#else
2.53 + template <typename Graph,
2.54 + typename LengthMap = typename Graph::template EdgeMap<int> >
2.55 +#endif
2.56 +
2.57 + class MinMeanCycle
2.58 + {
2.59 + typedef typename Graph::Node Node;
2.60 + typedef typename Graph::NodeIt NodeIt;
2.61 + typedef typename Graph::Edge Edge;
2.62 + typedef typename Graph::EdgeIt EdgeIt;
2.63 + typedef typename Graph::InEdgeIt InEdgeIt;
2.64 + typedef typename Graph::OutEdgeIt OutEdgeIt;
2.65 +
2.66 + typedef typename LengthMap::Value Length;
2.67 +
2.68 + typedef typename Graph::template NodeMap<int> IntNodeMap;
2.69 + typedef typename Graph::template NodeMap<Edge> PredNodeMap;
2.70 + typedef Path<Graph> Path;
2.71 + typedef std::vector<Node> NodeVector;
2.72 + typedef typename NodeVector::iterator NodeVectorIt;
2.73 +
2.74 + protected:
2.75 +
2.76 + /// \brief Data sturcture for path data.
2.77 + struct PathData {
2.78 + bool found;
2.79 + Length dist;
2.80 + Edge pred;
2.81 + PathData(bool _found = false, Length _dist = 0) :
2.82 + found(_found), dist(_dist), pred(INVALID) {}
2.83 + PathData(bool _found, Length _dist, Edge _pred) :
2.84 + found(_found), dist(_dist), pred(_pred) {}
2.85 + };
2.86 +
2.87 + private:
2.88 +
2.89 + typedef typename Graph::template NodeMap<std::vector<PathData> > PathVectorNodeMap;
2.90 +
2.91 + protected:
2.92 +
2.93 + /// \brief Node map for storing path data.
2.94 + ///
2.95 + /// Node map for storing path data of all nodes in the current
2.96 + /// component. dmap[v][k] is the length of a shortest directed walk
2.97 + /// to node v from the starting node containing exactly k edges.
2.98 + PathVectorNodeMap dmap;
2.99 +
2.100 + /// \brief The directed graph the algorithm runs on.
2.101 + const Graph& graph;
2.102 + /// \brief The length of the edges.
2.103 + const LengthMap& length;
2.104 +
2.105 + /// \brief The total length of the found cycle.
2.106 + Length cycle_length;
2.107 + /// \brief The number of edges in the found cycle.
2.108 + int cycle_size;
2.109 + /// \brief The found cycle.
2.110 + Path *cycle_path;
2.111 + /// \brief The found cycle.
2.112 + bool local_path;
2.113 +
2.114 + /// \brief Node map for identifying strongly connected components.
2.115 + IntNodeMap comp;
2.116 + /// \brief The number of strongly connected components.
2.117 + int comp_num;
2.118 + /// \brief %Counter for identifying the current component.
2.119 + int comp_cnt;
2.120 + /// \brief Nodes of the current component.
2.121 + NodeVector nodes;
2.122 + /// \brief The processed nodes in the last round.
2.123 + NodeVector process;
2.124 +
2.125 + public :
2.126 +
2.127 + /// \brief The constructor of the class.
2.128 + ///
2.129 + /// The constructor of the class.
2.130 + ///
2.131 + /// \param _graph The directed graph the algorithm runs on.
2.132 + /// \param _length The length (cost) of the edges.
2.133 + MinMeanCycle( const Graph& _graph,
2.134 + const LengthMap& _length ) :
2.135 + graph(_graph), length(_length), dmap(_graph),
2.136 + cycle_length(0), cycle_size(-1), comp(_graph),
2.137 + cycle_path(NULL), local_path(false)
2.138 + { }
2.139 +
2.140 + /// \brief The destructor of the class.
2.141 + ~MinMeanCycle() {
2.142 + if (local_path) delete cycle_path;
2.143 + }
2.144 +
2.145 + protected:
2.146 +
2.147 + /// \brief Initializes the internal data structures.
2.148 + void init() {
2.149 + comp_num = stronglyConnectedComponents(graph, comp);
2.150 + if (!cycle_path) {
2.151 + local_path = true;
2.152 + cycle_path = new Path;
2.153 + }
2.154 + }
2.155 +
2.156 + /// \brief Initializes the internal data structures for the current
2.157 + /// component.
2.158 + void initCurrent() {
2.159 + nodes.clear();
2.160 + // Finding the nodes of the current component
2.161 + for (NodeIt v(graph); v != INVALID; ++v) {
2.162 + if (comp[v] == comp_cnt) nodes.push_back(v);
2.163 + }
2.164 +
2.165 + // Creating vectors for all nodes
2.166 + int n = nodes.size();
2.167 + for (NodeVectorIt vi = nodes.begin(); vi != nodes.end(); ++vi) {
2.168 + dmap[*vi].resize(n + 1);
2.169 + }
2.170 + }
2.171 +
2.172 + /// \brief Processes all rounds of computing required path data.
2.173 + void processRounds() {
2.174 + dmap[nodes[0]][0] = PathData(true, 0);
2.175 + process.clear();
2.176 + for (OutEdgeIt oe(graph, nodes[0]); oe != INVALID; ++oe) {
2.177 + Node v = graph.target(oe);
2.178 + if (comp[v] != comp_cnt || v == nodes[0]) continue;
2.179 + dmap[v][1] = PathData(true, length[oe], oe);
2.180 + process.push_back(v);
2.181 + }
2.182 + int n = nodes.size(), k;
2.183 + for (k = 2; k <= n && process.size() < n; ++k) {
2.184 + processNextBuildRound(k);
2.185 + }
2.186 + for ( ; k <= n; ++k) {
2.187 + processNextFullRound(k);
2.188 + }
2.189 + }
2.190 +
2.191 + /// \brief Processes one round of computing required path data and
2.192 + /// rebuilds \ref process vector.
2.193 + void processNextBuildRound(int k) {
2.194 + NodeVector next;
2.195 + for (NodeVectorIt ui = process.begin(); ui != process.end(); ++ui) {
2.196 + for (OutEdgeIt oe(graph, *ui); oe != INVALID; ++oe) {
2.197 + Node v = graph.target(oe);
2.198 + if (comp[v] != comp_cnt) continue;
2.199 + if ( !dmap[v][k].found || (dmap[v][k].dist > dmap[*ui][k-1].dist + length[oe]) ) {
2.200 + if (!dmap[v][k].found) next.push_back(v);
2.201 + dmap[v][k] = PathData(true, dmap[*ui][k-1].dist + length[oe], oe);
2.202 + }
2.203 + }
2.204 + }
2.205 + process.swap(next);
2.206 + }
2.207 +
2.208 + /// \brief Processes one round of computing required path data
2.209 + /// using \ref nodes vector instead of \ref process vector.
2.210 + void processNextFullRound(int k) {
2.211 + for (NodeVectorIt ui = nodes.begin(); ui != nodes.end(); ++ui) {
2.212 + for (OutEdgeIt oe(graph, *ui); oe != INVALID; ++oe) {
2.213 + Node v = graph.target(oe);
2.214 + if (comp[v] != comp_cnt) continue;
2.215 + if ( !dmap[v][k].found || (dmap[v][k].dist > dmap[*ui][k-1].dist + length[oe]) ) {
2.216 + dmap[v][k] = PathData(true, dmap[*ui][k-1].dist + length[oe], oe);
2.217 + }
2.218 + }
2.219 + }
2.220 + }
2.221 +
2.222 + /// \brief Finds the minimum cycle mean value according to the path
2.223 + /// data stored in \ref dmap.
2.224 + ///
2.225 + /// Finds the minimum cycle mean value according to the path data
2.226 + /// stored in \ref dmap.
2.227 + ///
2.228 + /// \retval min_length The total length of the found cycle.
2.229 + /// \retval min_size The number of edges in the found cycle.
2.230 + /// \retval min_node A node for obtaining a critical cycle.
2.231 + ///
2.232 + /// \return \c true if a cycle exists in the graph.
2.233 + bool findMinCycleMean(Length &min_length, int &min_size, Node &min_node) {
2.234 + bool found_min = false;
2.235 + for (NodeVectorIt vi = nodes.begin(); vi != nodes.end(); ++vi) {
2.236 + Length len;
2.237 + int size;
2.238 + bool found_one = false;
2.239 + int n = nodes.size();
2.240 + for (int k = 0; k < n; ++k) {
2.241 + Length _len = dmap[*vi][n].dist - dmap[*vi][k].dist;
2.242 + int _size = n - k;
2.243 + if ( dmap[*vi][n].found && dmap[*vi][k].found && (!found_one || len * _size < _len * size) ) {
2.244 + found_one = true;
2.245 + len = _len;
2.246 + size = _size;
2.247 + }
2.248 + }
2.249 + if (found_one && (!found_min || len * min_size < min_length * size)) {
2.250 + found_min = true;
2.251 + min_length = len;
2.252 + min_size = size;
2.253 + min_node = *vi;
2.254 + }
2.255 + }
2.256 + return found_min;
2.257 + }
2.258 +
2.259 + /// \brief Finds a critical (minimum mean) cycle according to the
2.260 + /// path data stored in \ref dmap.
2.261 + void findCycle(const Node &min_n) {
2.262 + IntNodeMap reached(graph, -1);
2.263 + int r = reached[min_n] = dmap[min_n].size() - 1;
2.264 + Node u = graph.source(dmap[min_n][r].pred);
2.265 + while (reached[u] < 0) {
2.266 + reached[u] = --r;
2.267 + u = graph.source(dmap[u][r].pred);
2.268 + }
2.269 + r = reached[u];
2.270 + Edge ce = dmap[u][r].pred;
2.271 + cycle_path->addFront(ce);
2.272 + Node v;
2.273 + while ((v = graph.source(ce)) != u) {
2.274 + ce = dmap[v][--r].pred;
2.275 + cycle_path->addFront(ce);
2.276 + }
2.277 + }
2.278 +
2.279 + public:
2.280 +
2.281 + /// \brief Runs the algorithm.
2.282 + ///
2.283 + /// Runs the algorithm.
2.284 + ///
2.285 + /// \return \c true if a cycle exists in the graph.
2.286 + bool run() {
2.287 + init();
2.288 + // Searching for the minimum mean cycle in all components
2.289 + bool found_cycle = false;
2.290 + Node cycle_node;
2.291 + for (comp_cnt = 0; comp_cnt < comp_num; ++comp_cnt) {
2.292 +
2.293 + initCurrent();
2.294 + processRounds();
2.295 +
2.296 + Length min_length;
2.297 + int min_size;
2.298 + Node min_node;
2.299 + bool found_min = findMinCycleMean(min_length, min_size, min_node);
2.300 +
2.301 + if ( found_min && (!found_cycle || min_length * cycle_size < cycle_length * min_size) ) {
2.302 + found_cycle = true;
2.303 + cycle_length = min_length;
2.304 + cycle_size = min_size;
2.305 + cycle_node = min_node;
2.306 + }
2.307 + }
2.308 + if (found_cycle) findCycle(cycle_node);
2.309 + return found_cycle;
2.310 + }
2.311 +
2.312 + /// \brief Returns the total length of the found critical cycle.
2.313 + ///
2.314 + /// Returns the total length of the found critical cycle.
2.315 + ///
2.316 + /// \pre \ref run() must be called before using this function.
2.317 + Length cycleLength() const {
2.318 + return cycle_length;
2.319 + }
2.320 +
2.321 + /// \brief Returns the number of edges in the found critical
2.322 + /// cycle.
2.323 + ///
2.324 + /// Returns the number of edges in the found critical cycle.
2.325 + ///
2.326 + /// \pre \ref run() must be called before using this function.
2.327 + int cycleEdgeNum() const {
2.328 + return cycle_size;
2.329 + }
2.330 +
2.331 + /// \brief Returns the mean length of the found critical cycle.
2.332 + ///
2.333 + /// Returns the mean length of the found critical cycle.
2.334 + ///
2.335 + /// \pre \ref run() must be called before using this function.
2.336 + ///
2.337 + /// \warning LengthMap::Value must be convertible to double.
2.338 + double minMean() const {
2.339 + return (double)cycle_length / (double)cycle_size;
2.340 + }
2.341 +
2.342 + /// \brief Returns a const reference to the \ref lemon::Path "Path"
2.343 + /// structure of the found flow.
2.344 + ///
2.345 + /// Returns a const reference to the \ref lemon::Path "Path"
2.346 + /// structure of the found flow.
2.347 + ///
2.348 + /// \pre \ref run() must be called before using this function.
2.349 + ///
2.350 + /// \sa \ref cyclePath()
2.351 + const Path &cycle() const {
2.352 + return *cycle_path;
2.353 + }
2.354 +
2.355 + /// \brief Sets the \ref lemon::Path "Path" structure storing the
2.356 + /// found cycle.
2.357 + ///
2.358 + /// Sets the \ref lemon::Path "Path" structure storing the found
2.359 + /// cycle. If you don't use this function before calling
2.360 + /// \ref run(), it will allocate one. The destuctor deallocates
2.361 + /// this automatically allocated map, of course.
2.362 + ///
2.363 + /// \note The algorithm calls only the \ref lemon::Path::addFront()
2.364 + /// "addFront()" method of the given \ref lemon::Path "Path"
2.365 + /// structure.
2.366 + ///
2.367 + /// \return \c (*this)
2.368 + MinMeanCycle &cyclePath(Path& path) {
2.369 + if (local_path) {
2.370 + delete cycle_path;
2.371 + local_path = false;
2.372 + }
2.373 + cycle_path = &path;
2.374 + return *this;
2.375 + }
2.376 +
2.377 + }; //class MinMeanCycle
2.378 +
2.379 + ///@}
2.380 +
2.381 +} //namespace lemon
2.382 +
2.383 +#endif //LEMON_MIN_MEAN_CYCLE_H