source:lemon-0.x/lemon/min_mean_cycle.h@2555:a84e52e99f57

Last change on this file since 2555:a84e52e99f57 was 2555:a84e52e99f57, checked in by Peter Kovacs, 15 years ago

Reimplemented MinMeanCycle? to be much more efficient.
The new version implements Howard's algorithm instead of Karp's algorithm and
it is at least 10-20 times faster on all the 40-50 random graphs we have tested.

File size: 13.1 KB
Line
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_MIN_MEAN_CYCLE_H
20#define LEMON_MIN_MEAN_CYCLE_H
21
22/// \ingroup shortest_path
23///
24/// \file
25/// \brief Howard's algorithm for finding a minimum mean directed cycle.
26
27#include <vector>
28#include <lemon/graph_utils.h>
29#include <lemon/topology.h>
30#include <lemon/path.h>
31#include <lemon/tolerance.h>
32
33namespace lemon {
34
36  /// @{
37
38  /// \brief Implementation of Howard's algorithm for finding a minimum
39  /// mean directed cycle.
40  ///
41  /// \ref MinMeanCycle implements Howard's algorithm for finding a
42  /// minimum mean directed cycle.
43  ///
44  /// \param Graph The directed graph type the algorithm runs on.
45  /// \param LengthMap The type of the length (cost) map.
46  ///
47  /// \warning \c LengthMap::Value must be convertible to \c double.
48  ///
49  /// \author Peter Kovacs
50
51  template < typename Graph,
52             typename LengthMap = typename Graph::template EdgeMap<int> >
53  class MinMeanCycle
54  {
55    GRAPH_TYPEDEFS(typename Graph);
56
57    typedef typename LengthMap::Value Length;
58    typedef Path<Graph> Path;
59
60  protected:
61
62    typename Graph::template NodeMap<bool> _reached;
63    typename Graph::template NodeMap<double> _dist;
64    typename Graph::template NodeMap<Edge> _policy;
65
66    /// The directed graph the algorithm runs on.
67    const Graph &_graph;
68    /// The length of the edges.
69    const LengthMap &_length;
70
71    /// The total length of the found cycle.
72    Length _cycle_length;
73    /// The number of edges on the found cycle.
74    int _cycle_size;
75    /// The found cycle.
76    Path *_cycle_path;
77
78    bool _local_path;
79    bool _cycle_found;
80    Node _cycle_node;
81
82    typename Graph::template NodeMap<int> _component;
83    int _component_num;
84
85    std::vector<Node> _nodes;
86    std::vector<Edge> _edges;
87    Tolerance<double> _tolerance;
88
89  public:
90
91    /// \brief The constructor of the class.
92    ///
93    /// The constructor of the class.
94    ///
95    /// \param graph The directed graph the algorithm runs on.
96    /// \param length The length (cost) of the edges.
97    MinMeanCycle( const Graph &graph,
98                  const LengthMap &length ) :
99      _graph(graph), _length(length), _dist(graph), _reached(graph),
100      _policy(graph), _component(graph), _cycle_length(0),
101      _cycle_size(-1), _cycle_path(NULL), _local_path(false)
102    { }
103
104    /// The destructor of the class.
105    ~MinMeanCycle() {
106      if (_local_path) delete _cycle_path;
107    }
108
109  protected:
110
111    // Initializes the internal data structures for the current strongly
112    // connected component and creating the policy graph.
113    // The policy graph can be represented by the _policy map because
114    // the out degree of every node is 1.
115    bool initCurrentComponent(int comp) {
116      // Finding the nodes of the current component
117      _nodes.clear();
118      for (NodeIt n(_graph); n != INVALID; ++n) {
119        if (_component[n] == comp) _nodes.push_back(n);
120      }
121      if (_nodes.size() <= 1) return false;
122      // Finding the edges of the current component
123      _edges.clear();
124      for (EdgeIt e(_graph); e != INVALID; ++e) {
125        if ( _component[_graph.source(e)] == comp &&
126             _component[_graph.target(e)] == comp )
127          _edges.push_back(e);
128      }
129      // Initializing _reached, _dist, _policy maps
130      for (int i = 0; i < _nodes.size(); ++i) {
131        _reached[_nodes[i]] = false;
132        _policy[_nodes[i]] = INVALID;
133      }
134      Node u; Edge e;
135      for (int j = 0; j < _edges.size(); ++j) {
136        e = _edges[j];
137        u = _graph.source(e);
138        if (!_reached[u] || _length[e] < _dist[u]) {
139          _dist[u] = _length[e];
140          _policy[u] = e;
141          _reached[u] = true;
142        }
143      }
144      return true;
145    }
146
147    // Finds all cycles in the policy graph.
148    // Sets _cycle_found to true if a cycle is found and sets
149    // _cycle_length, _cycle_size, _cycle_node to represent the minimum
150    // mean cycle in the policy graph.
151    bool findPolicyCycles(int comp) {
152      typename Graph::template NodeMap<int> level(_graph, -1);
153      _cycle_found = false;
154      Length clength;
155      int csize;
156      int path_cnt = 0;
157      Node u, v;
158      // Searching for cycles
159      for (int i = 0; i < _nodes.size(); ++i) {
160        if (level[_nodes[i]] < 0) {
161          u = _nodes[i];
162          level[u] = path_cnt;
163          while (level[u = _graph.target(_policy[u])] < 0)
164            level[u] = path_cnt;
165          if (level[u] == path_cnt) {
166            // A cycle is found
167            clength = _length[_policy[u]];
168            csize = 1;
169            for (v = u; (v = _graph.target(_policy[v])) != u; ) {
170              clength += _length[_policy[v]];
171              ++csize;
172            }
173            if ( !_cycle_found ||
174                 clength * _cycle_size < _cycle_length * csize ) {
175              _cycle_found = true;
176              _cycle_length = clength;
177              _cycle_size = csize;
178              _cycle_node = u;
179            }
180          }
181          ++path_cnt;
182        }
183      }
184      return _cycle_found;
185    }
186
187    // Contracts the policy graph to be connected by cutting all cycles
188    // except for the main cycle (i.e. the minimum mean cycle).
189    void contractPolicyGraph(int comp) {
190      // Finding the component of the main cycle using
191      // reverse BFS search
192      typename Graph::template NodeMap<int> found(_graph, false);
193      std::deque<Node> queue;
194      queue.push_back(_cycle_node);
195      found[_cycle_node] = true;
196      Node u, v;
197      while (!queue.empty()) {
198        v = queue.front(); queue.pop_front();
199        for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
200          u = _graph.source(e);
201          if (_component[u] == comp && !found[u] && _policy[u] == e) {
202            found[u] = true;
203            queue.push_back(u);
204          }
205        }
206      }
207      // Connecting all other nodes to this component using
208      // reverse BFS search
209      queue.clear();
210      for (int i = 0; i < _nodes.size(); ++i)
211        if (found[_nodes[i]]) queue.push_back(_nodes[i]);
212      int found_cnt = queue.size();
213      while (found_cnt < _nodes.size() && !queue.empty()) {
214        v = queue.front(); queue.pop_front();
215        for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
216          u = _graph.source(e);
217          if (_component[u] == comp && !found[u]) {
218            found[u] = true;
219            ++found_cnt;
220            _policy[u] = e;
221            queue.push_back(u);
222          }
223        }
224      }
225    }
226
227    // Computes node distances in the policy graph and updates the
228    // policy graph if the node distances can be improved.
229    bool computeNodeDistances(int comp) {
230      // Computing node distances using reverse BFS search
231      double cycle_mean = (double)_cycle_length / _cycle_size;
232      typename Graph::template NodeMap<int> found(_graph, false);
233      std::deque<Node> queue;
234      queue.push_back(_cycle_node);
235      found[_cycle_node] = true;
236      Node u, v;
237      while (!queue.empty()) {
238        v = queue.front(); queue.pop_front();
239        for (InEdgeIt e(_graph, v); e != INVALID; ++e) {
240          u = _graph.source(e);
241          if (_component[u] == comp && !found[u] && _policy[u] == e) {
242            found[u] = true;
243            _dist[u] = _dist[v] + _length[e] - cycle_mean;
244            queue.push_back(u);
245          }
246        }
247      }
248      // Improving node distances
249      bool improved = false;
250      for (int j = 0; j < _edges.size(); ++j) {
251        Edge e = _edges[j];
252        u = _graph.source(e); v = _graph.target(e);
253        double delta = _dist[v] + _length[e] - cycle_mean;
254        if (_tolerance.less(delta, _dist[u])) {
255          improved = true;
256          _dist[u] = delta;
257          _policy[u] = e;
258        }
259      }
260      return improved;
261    }
262
263  public:
264
265    /// \brief Runs the algorithm.
266    ///
267    /// Runs the algorithm.
268    ///
269    /// \return Returns \c true if a directed cycle exists in the graph.
270    ///
271    /// \note Apart from the return value, <tt>mmc.run()</tt> is just a
272    /// shortcut of the following code.
273    /// \code
274    ///   mmc.init();
275    ///   mmc.findMinMean();
276    ///   mmc.findCycle();
277    /// \endcode
278    bool run() {
279      init();
280      return findMinMean() && findCycle();
281    }
282
283    /// \brief Initializes the internal data structures.
284    ///
285    /// Initializes the internal data structures.
286    ///
287    /// \sa reset()
288    void init() {
289      _tolerance.epsilon(1e-8);
290      if (!_cycle_path) {
291        _local_path = true;
292        _cycle_path = new Path;
293      }
294      _cycle_found = false;
295      _component_num = stronglyConnectedComponents(_graph, _component);
296    }
297
298    /// \brief Resets the internal data structures.
299    ///
300    /// Resets the internal data structures so that \ref findMinMean()
301    /// and \ref findCycle() can be called again (e.g. when the
302    /// underlaying graph has been modified).
303    ///
304    /// \sa init()
305    void reset() {
306      if (_cycle_path) _cycle_path->clear();
307      _cycle_found = false;
308      _component_num = stronglyConnectedComponents(_graph, _component);
309    }
310
311    /// \brief Finds the minimum cycle mean length in the graph.
312    ///
313    /// Computes all the required data and finds the minimum cycle mean
314    /// length in the graph.
315    ///
316    /// \return Returns \c true if a directed cycle exists in the graph.
317    ///
318    /// \pre \ref init() must be called before using this function.
319    bool findMinMean() {
320      // Finding the minimum mean cycle in the components
321      for (int comp = 0; comp < _component_num; ++comp) {
322        if (!initCurrentComponent(comp)) continue;
323        while (true) {
324          if (!findPolicyCycles(comp)) return false;
325          contractPolicyGraph(comp);
326          if (!computeNodeDistances(comp)) return true;
327        }
328      }
329    }
330
331    /// \brief Finds a critical (minimum mean) directed cycle.
332    ///
333    /// Finds a critical (minimum mean) directed cycle using the data
334    /// computed in the \ref findMinMean() function.
335    ///
336    /// \return Returns \c true if a directed cycle exists in the graph.
337    ///
338    /// \pre \ref init() and \ref findMinMean() must be called before
339    /// using this function.
340    bool findCycle() {
341      if (!_cycle_found) return false;
343      for ( Node v = _cycle_node;
344            (v = _graph.target(_policy[v])) != _cycle_node; ) {
346      }
347      return true;
348    }
349
350    /// \brief Returns the total length of the found cycle.
351    ///
352    /// Returns the total length of the found cycle.
353    ///
354    /// \pre \ref run() or \ref findMinMean() must be called before
355    /// using this function.
356    Length cycleLength() const {
357      return _cycle_length;
358    }
359
360    /// \brief Returns the number of edges on the found cycle.
361    ///
362    /// Returns the number of edges on the found cycle.
363    ///
364    /// \pre \ref run() or \ref findMinMean() must be called before
365    /// using this function.
366    int cycleEdgeNum() const {
367      return _cycle_size;
368    }
369
370    /// \brief Returns the mean length of the found cycle.
371    ///
372    /// Returns the mean length of the found cycle.
373    ///
374    /// \pre \ref run() or \ref findMinMean() must be called before
375    /// using this function.
376    ///
377    /// \note <tt>mmc.cycleMean()</tt> is just a shortcut of the
378    /// following code.
379    /// \code
380    ///   return double(mmc.cycleLength()) / mmc.cycleEdgeNum();
381    /// \endcode
382    double cycleMean() const {
383      return (double)_cycle_length / _cycle_size;
384    }
385
386    /// \brief Returns a const reference to the \ref Path "path"
387    /// structure storing the found cycle.
388    ///
389    /// Returns a const reference to the \ref Path "path"
390    /// structure storing the found cycle.
391    ///
392    /// \pre \ref run() or \ref findCycle() must be called before using
393    /// this function.
394    ///
395    /// \sa cyclePath()
396    const Path& cycle() const {
397      return *_cycle_path;
398    }
399
400    /// \brief Sets the \ref Path "path" structure for storing the found
401    /// cycle.
402    ///
403    /// Sets an external \ref Path "path" structure for storing the
404    /// found cycle.
405    ///
406    /// If you don't call this function before calling \ref run() or
407    /// \ref init(), it will allocate a local \ref Path "path"
408    /// structure.
409    /// The destuctor deallocates this automatically allocated map,
410    /// of course.
411    ///
412    /// \note The algorithm calls only the \ref lemon::Path::addBack()
413    /// "addBack()" function of the given \ref Path "path" structure.
414    ///
415    /// \return <tt>(*this)</tt>
416    ///
417    /// \sa cycle()
418    MinMeanCycle& cyclePath(Path &path) {
419      if (_local_path) {
420        delete _cycle_path;
421        _local_path = false;
422      }
423      _cycle_path = &path;
424      return *this;
425    }
426
427  }; //class MinMeanCycle
428
429  ///@}
430
431} //namespace lemon
432
433#endif //LEMON_MIN_MEAN_CYCLE_H
Note: See TracBrowser for help on using the repository browser.