1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/hartmann_orlin_mmc.h Sat Mar 13 22:01:38 2010 +0100
1.3 @@ -0,0 +1,650 @@
1.4 +/* -*- C++ -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library
1.7 + *
1.8 + * Copyright (C) 2003-2008
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#ifndef LEMON_HARTMANN_ORLIN_MMC_H
1.23 +#define LEMON_HARTMANN_ORLIN_MMC_H
1.24 +
1.25 +/// \ingroup min_mean_cycle
1.26 +///
1.27 +/// \file
1.28 +/// \brief Hartmann-Orlin's algorithm for finding a minimum mean cycle.
1.29 +
1.30 +#include <vector>
1.31 +#include <limits>
1.32 +#include <lemon/core.h>
1.33 +#include <lemon/path.h>
1.34 +#include <lemon/tolerance.h>
1.35 +#include <lemon/connectivity.h>
1.36 +
1.37 +namespace lemon {
1.38 +
1.39 + /// \brief Default traits class of HartmannOrlinMmc class.
1.40 + ///
1.41 + /// Default traits class of HartmannOrlinMmc class.
1.42 + /// \tparam GR The type of the digraph.
1.43 + /// \tparam CM The type of the cost map.
1.44 + /// It must conform to the \ref concepts::Rea_data "Rea_data" concept.
1.45 +#ifdef DOXYGEN
1.46 + template <typename GR, typename CM>
1.47 +#else
1.48 + template <typename GR, typename CM,
1.49 + bool integer = std::numeric_limits<typename CM::Value>::is_integer>
1.50 +#endif
1.51 + struct HartmannOrlinMmcDefaultTraits
1.52 + {
1.53 + /// The type of the digraph
1.54 + typedef GR Digraph;
1.55 + /// The type of the cost map
1.56 + typedef CM CostMap;
1.57 + /// The type of the arc costs
1.58 + typedef typename CostMap::Value Cost;
1.59 +
1.60 + /// \brief The large cost type used for internal computations
1.61 + ///
1.62 + /// The large cost type used for internal computations.
1.63 + /// It is \c long \c long if the \c Cost type is integer,
1.64 + /// otherwise it is \c double.
1.65 + /// \c Cost must be convertible to \c LargeCost.
1.66 + typedef double LargeCost;
1.67 +
1.68 + /// The tolerance type used for internal computations
1.69 + typedef lemon::Tolerance<LargeCost> Tolerance;
1.70 +
1.71 + /// \brief The path type of the found cycles
1.72 + ///
1.73 + /// The path type of the found cycles.
1.74 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
1.75 + /// and it must have an \c addFront() function.
1.76 + typedef lemon::Path<Digraph> Path;
1.77 + };
1.78 +
1.79 + // Default traits class for integer cost types
1.80 + template <typename GR, typename CM>
1.81 + struct HartmannOrlinMmcDefaultTraits<GR, CM, true>
1.82 + {
1.83 + typedef GR Digraph;
1.84 + typedef CM CostMap;
1.85 + typedef typename CostMap::Value Cost;
1.86 +#ifdef LEMON_HAVE_LONG_LONG
1.87 + typedef long long LargeCost;
1.88 +#else
1.89 + typedef long LargeCost;
1.90 +#endif
1.91 + typedef lemon::Tolerance<LargeCost> Tolerance;
1.92 + typedef lemon::Path<Digraph> Path;
1.93 + };
1.94 +
1.95 +
1.96 + /// \addtogroup min_mean_cycle
1.97 + /// @{
1.98 +
1.99 + /// \brief Implementation of the Hartmann-Orlin algorithm for finding
1.100 + /// a minimum mean cycle.
1.101 + ///
1.102 + /// This class implements the Hartmann-Orlin algorithm for finding
1.103 + /// a directed cycle of minimum mean cost in a digraph
1.104 + /// \ref amo93networkflows, \ref dasdan98minmeancycle.
1.105 + /// It is an improved version of \ref Karp "Karp"'s original algorithm,
1.106 + /// it applies an efficient early termination scheme.
1.107 + /// It runs in time O(ne) and uses space O(n<sup>2</sup>+e).
1.108 + ///
1.109 + /// \tparam GR The type of the digraph the algorithm runs on.
1.110 + /// \tparam CM The type of the cost map. The default
1.111 + /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
1.112 + /// \tparam TR The traits class that defines various types used by the
1.113 + /// algorithm. By default, it is \ref HartmannOrlinMmcDefaultTraits
1.114 + /// "HartmannOrlinMmcDefaultTraits<GR, CM>".
1.115 + /// In most cases, this parameter should not be set directly,
1.116 + /// consider to use the named template parameters instead.
1.117 +#ifdef DOXYGEN
1.118 + template <typename GR, typename CM, typename TR>
1.119 +#else
1.120 + template < typename GR,
1.121 + typename CM = typename GR::template ArcMap<int>,
1.122 + typename TR = HartmannOrlinMmcDefaultTraits<GR, CM> >
1.123 +#endif
1.124 + class HartmannOrlinMmc
1.125 + {
1.126 + public:
1.127 +
1.128 + /// The type of the digraph
1.129 + typedef typename TR::Digraph Digraph;
1.130 + /// The type of the cost map
1.131 + typedef typename TR::CostMap CostMap;
1.132 + /// The type of the arc costs
1.133 + typedef typename TR::Cost Cost;
1.134 +
1.135 + /// \brief The large cost type
1.136 + ///
1.137 + /// The large cost type used for internal computations.
1.138 + /// By default, it is \c long \c long if the \c Cost type is integer,
1.139 + /// otherwise it is \c double.
1.140 + typedef typename TR::LargeCost LargeCost;
1.141 +
1.142 + /// The tolerance type
1.143 + typedef typename TR::Tolerance Tolerance;
1.144 +
1.145 + /// \brief The path type of the found cycles
1.146 + ///
1.147 + /// The path type of the found cycles.
1.148 + /// Using the \ref HartmannOrlinMmcDefaultTraits "default traits class",
1.149 + /// it is \ref lemon::Path "Path<Digraph>".
1.150 + typedef typename TR::Path Path;
1.151 +
1.152 + /// The \ref HartmannOrlinMmcDefaultTraits "traits class" of the algorithm
1.153 + typedef TR Traits;
1.154 +
1.155 + private:
1.156 +
1.157 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
1.158 +
1.159 + // Data sturcture for path data
1.160 + struct PathData
1.161 + {
1.162 + LargeCost dist;
1.163 + Arc pred;
1.164 + PathData(LargeCost d, Arc p = INVALID) :
1.165 + dist(d), pred(p) {}
1.166 + };
1.167 +
1.168 + typedef typename Digraph::template NodeMap<std::vector<PathData> >
1.169 + PathDataNodeMap;
1.170 +
1.171 + private:
1.172 +
1.173 + // The digraph the algorithm runs on
1.174 + const Digraph &_gr;
1.175 + // The cost of the arcs
1.176 + const CostMap &_cost;
1.177 +
1.178 + // Data for storing the strongly connected components
1.179 + int _comp_num;
1.180 + typename Digraph::template NodeMap<int> _comp;
1.181 + std::vector<std::vector<Node> > _comp_nodes;
1.182 + std::vector<Node>* _nodes;
1.183 + typename Digraph::template NodeMap<std::vector<Arc> > _out_arcs;
1.184 +
1.185 + // Data for the found cycles
1.186 + bool _curr_found, _best_found;
1.187 + LargeCost _curr_cost, _best_cost;
1.188 + int _curr_size, _best_size;
1.189 + Node _curr_node, _best_node;
1.190 + int _curr_level, _best_level;
1.191 +
1.192 + Path *_cycle_path;
1.193 + bool _local_path;
1.194 +
1.195 + // Node map for storing path data
1.196 + PathDataNodeMap _data;
1.197 + // The processed nodes in the last round
1.198 + std::vector<Node> _process;
1.199 +
1.200 + Tolerance _tolerance;
1.201 +
1.202 + // Infinite constant
1.203 + const LargeCost INF;
1.204 +
1.205 + public:
1.206 +
1.207 + /// \name Named Template Parameters
1.208 + /// @{
1.209 +
1.210 + template <typename T>
1.211 + struct SetLargeCostTraits : public Traits {
1.212 + typedef T LargeCost;
1.213 + typedef lemon::Tolerance<T> Tolerance;
1.214 + };
1.215 +
1.216 + /// \brief \ref named-templ-param "Named parameter" for setting
1.217 + /// \c LargeCost type.
1.218 + ///
1.219 + /// \ref named-templ-param "Named parameter" for setting \c LargeCost
1.220 + /// type. It is used for internal computations in the algorithm.
1.221 + template <typename T>
1.222 + struct SetLargeCost
1.223 + : public HartmannOrlinMmc<GR, CM, SetLargeCostTraits<T> > {
1.224 + typedef HartmannOrlinMmc<GR, CM, SetLargeCostTraits<T> > Create;
1.225 + };
1.226 +
1.227 + template <typename T>
1.228 + struct SetPathTraits : public Traits {
1.229 + typedef T Path;
1.230 + };
1.231 +
1.232 + /// \brief \ref named-templ-param "Named parameter" for setting
1.233 + /// \c %Path type.
1.234 + ///
1.235 + /// \ref named-templ-param "Named parameter" for setting the \c %Path
1.236 + /// type of the found cycles.
1.237 + /// It must conform to the \ref lemon::concepts::Path "Path" concept
1.238 + /// and it must have an \c addFront() function.
1.239 + template <typename T>
1.240 + struct SetPath
1.241 + : public HartmannOrlinMmc<GR, CM, SetPathTraits<T> > {
1.242 + typedef HartmannOrlinMmc<GR, CM, SetPathTraits<T> > Create;
1.243 + };
1.244 +
1.245 + /// @}
1.246 +
1.247 + protected:
1.248 +
1.249 + HartmannOrlinMmc() {}
1.250 +
1.251 + public:
1.252 +
1.253 + /// \brief Constructor.
1.254 + ///
1.255 + /// The constructor of the class.
1.256 + ///
1.257 + /// \param digraph The digraph the algorithm runs on.
1.258 + /// \param cost The costs of the arcs.
1.259 + HartmannOrlinMmc( const Digraph &digraph,
1.260 + const CostMap &cost ) :
1.261 + _gr(digraph), _cost(cost), _comp(digraph), _out_arcs(digraph),
1.262 + _best_found(false), _best_cost(0), _best_size(1),
1.263 + _cycle_path(NULL), _local_path(false), _data(digraph),
1.264 + INF(std::numeric_limits<LargeCost>::has_infinity ?
1.265 + std::numeric_limits<LargeCost>::infinity() :
1.266 + std::numeric_limits<LargeCost>::max())
1.267 + {}
1.268 +
1.269 + /// Destructor.
1.270 + ~HartmannOrlinMmc() {
1.271 + if (_local_path) delete _cycle_path;
1.272 + }
1.273 +
1.274 + /// \brief Set the path structure for storing the found cycle.
1.275 + ///
1.276 + /// This function sets an external path structure for storing the
1.277 + /// found cycle.
1.278 + ///
1.279 + /// If you don't call this function before calling \ref run() or
1.280 + /// \ref findCycleMean(), it will allocate a local \ref Path "path"
1.281 + /// structure. The destuctor deallocates this automatically
1.282 + /// allocated object, of course.
1.283 + ///
1.284 + /// \note The algorithm calls only the \ref lemon::Path::addFront()
1.285 + /// "addFront()" function of the given path structure.
1.286 + ///
1.287 + /// \return <tt>(*this)</tt>
1.288 + HartmannOrlinMmc& cycle(Path &path) {
1.289 + if (_local_path) {
1.290 + delete _cycle_path;
1.291 + _local_path = false;
1.292 + }
1.293 + _cycle_path = &path;
1.294 + return *this;
1.295 + }
1.296 +
1.297 + /// \brief Set the tolerance used by the algorithm.
1.298 + ///
1.299 + /// This function sets the tolerance object used by the algorithm.
1.300 + ///
1.301 + /// \return <tt>(*this)</tt>
1.302 + HartmannOrlinMmc& tolerance(const Tolerance& tolerance) {
1.303 + _tolerance = tolerance;
1.304 + return *this;
1.305 + }
1.306 +
1.307 + /// \brief Return a const reference to the tolerance.
1.308 + ///
1.309 + /// This function returns a const reference to the tolerance object
1.310 + /// used by the algorithm.
1.311 + const Tolerance& tolerance() const {
1.312 + return _tolerance;
1.313 + }
1.314 +
1.315 + /// \name Execution control
1.316 + /// The simplest way to execute the algorithm is to call the \ref run()
1.317 + /// function.\n
1.318 + /// If you only need the minimum mean cost, you may call
1.319 + /// \ref findCycleMean().
1.320 +
1.321 + /// @{
1.322 +
1.323 + /// \brief Run the algorithm.
1.324 + ///
1.325 + /// This function runs the algorithm.
1.326 + /// It can be called more than once (e.g. if the underlying digraph
1.327 + /// and/or the arc costs have been modified).
1.328 + ///
1.329 + /// \return \c true if a directed cycle exists in the digraph.
1.330 + ///
1.331 + /// \note <tt>mmc.run()</tt> is just a shortcut of the following code.
1.332 + /// \code
1.333 + /// return mmc.findCycleMean() && mmc.findCycle();
1.334 + /// \endcode
1.335 + bool run() {
1.336 + return findCycleMean() && findCycle();
1.337 + }
1.338 +
1.339 + /// \brief Find the minimum cycle mean.
1.340 + ///
1.341 + /// This function finds the minimum mean cost of the directed
1.342 + /// cycles in the digraph.
1.343 + ///
1.344 + /// \return \c true if a directed cycle exists in the digraph.
1.345 + bool findCycleMean() {
1.346 + // Initialization and find strongly connected components
1.347 + init();
1.348 + findComponents();
1.349 +
1.350 + // Find the minimum cycle mean in the components
1.351 + for (int comp = 0; comp < _comp_num; ++comp) {
1.352 + if (!initComponent(comp)) continue;
1.353 + processRounds();
1.354 +
1.355 + // Update the best cycle (global minimum mean cycle)
1.356 + if ( _curr_found && (!_best_found ||
1.357 + _curr_cost * _best_size < _best_cost * _curr_size) ) {
1.358 + _best_found = true;
1.359 + _best_cost = _curr_cost;
1.360 + _best_size = _curr_size;
1.361 + _best_node = _curr_node;
1.362 + _best_level = _curr_level;
1.363 + }
1.364 + }
1.365 + return _best_found;
1.366 + }
1.367 +
1.368 + /// \brief Find a minimum mean directed cycle.
1.369 + ///
1.370 + /// This function finds a directed cycle of minimum mean cost
1.371 + /// in the digraph using the data computed by findCycleMean().
1.372 + ///
1.373 + /// \return \c true if a directed cycle exists in the digraph.
1.374 + ///
1.375 + /// \pre \ref findCycleMean() must be called before using this function.
1.376 + bool findCycle() {
1.377 + if (!_best_found) return false;
1.378 + IntNodeMap reached(_gr, -1);
1.379 + int r = _best_level + 1;
1.380 + Node u = _best_node;
1.381 + while (reached[u] < 0) {
1.382 + reached[u] = --r;
1.383 + u = _gr.source(_data[u][r].pred);
1.384 + }
1.385 + r = reached[u];
1.386 + Arc e = _data[u][r].pred;
1.387 + _cycle_path->addFront(e);
1.388 + _best_cost = _cost[e];
1.389 + _best_size = 1;
1.390 + Node v;
1.391 + while ((v = _gr.source(e)) != u) {
1.392 + e = _data[v][--r].pred;
1.393 + _cycle_path->addFront(e);
1.394 + _best_cost += _cost[e];
1.395 + ++_best_size;
1.396 + }
1.397 + return true;
1.398 + }
1.399 +
1.400 + /// @}
1.401 +
1.402 + /// \name Query Functions
1.403 + /// The results of the algorithm can be obtained using these
1.404 + /// functions.\n
1.405 + /// The algorithm should be executed before using them.
1.406 +
1.407 + /// @{
1.408 +
1.409 + /// \brief Return the total cost of the found cycle.
1.410 + ///
1.411 + /// This function returns the total cost of the found cycle.
1.412 + ///
1.413 + /// \pre \ref run() or \ref findCycleMean() must be called before
1.414 + /// using this function.
1.415 + Cost cycleCost() const {
1.416 + return static_cast<Cost>(_best_cost);
1.417 + }
1.418 +
1.419 + /// \brief Return the number of arcs on the found cycle.
1.420 + ///
1.421 + /// This function returns the number of arcs on the found cycle.
1.422 + ///
1.423 + /// \pre \ref run() or \ref findCycleMean() must be called before
1.424 + /// using this function.
1.425 + int cycleSize() const {
1.426 + return _best_size;
1.427 + }
1.428 +
1.429 + /// \brief Return the mean cost of the found cycle.
1.430 + ///
1.431 + /// This function returns the mean cost of the found cycle.
1.432 + ///
1.433 + /// \note <tt>alg.cycleMean()</tt> is just a shortcut of the
1.434 + /// following code.
1.435 + /// \code
1.436 + /// return static_cast<double>(alg.cycleCost()) / alg.cycleSize();
1.437 + /// \endcode
1.438 + ///
1.439 + /// \pre \ref run() or \ref findCycleMean() must be called before
1.440 + /// using this function.
1.441 + double cycleMean() const {
1.442 + return static_cast<double>(_best_cost) / _best_size;
1.443 + }
1.444 +
1.445 + /// \brief Return the found cycle.
1.446 + ///
1.447 + /// This function returns a const reference to the path structure
1.448 + /// storing the found cycle.
1.449 + ///
1.450 + /// \pre \ref run() or \ref findCycle() must be called before using
1.451 + /// this function.
1.452 + const Path& cycle() const {
1.453 + return *_cycle_path;
1.454 + }
1.455 +
1.456 + ///@}
1.457 +
1.458 + private:
1.459 +
1.460 + // Initialization
1.461 + void init() {
1.462 + if (!_cycle_path) {
1.463 + _local_path = true;
1.464 + _cycle_path = new Path;
1.465 + }
1.466 + _cycle_path->clear();
1.467 + _best_found = false;
1.468 + _best_cost = 0;
1.469 + _best_size = 1;
1.470 + _cycle_path->clear();
1.471 + for (NodeIt u(_gr); u != INVALID; ++u)
1.472 + _data[u].clear();
1.473 + }
1.474 +
1.475 + // Find strongly connected components and initialize _comp_nodes
1.476 + // and _out_arcs
1.477 + void findComponents() {
1.478 + _comp_num = stronglyConnectedComponents(_gr, _comp);
1.479 + _comp_nodes.resize(_comp_num);
1.480 + if (_comp_num == 1) {
1.481 + _comp_nodes[0].clear();
1.482 + for (NodeIt n(_gr); n != INVALID; ++n) {
1.483 + _comp_nodes[0].push_back(n);
1.484 + _out_arcs[n].clear();
1.485 + for (OutArcIt a(_gr, n); a != INVALID; ++a) {
1.486 + _out_arcs[n].push_back(a);
1.487 + }
1.488 + }
1.489 + } else {
1.490 + for (int i = 0; i < _comp_num; ++i)
1.491 + _comp_nodes[i].clear();
1.492 + for (NodeIt n(_gr); n != INVALID; ++n) {
1.493 + int k = _comp[n];
1.494 + _comp_nodes[k].push_back(n);
1.495 + _out_arcs[n].clear();
1.496 + for (OutArcIt a(_gr, n); a != INVALID; ++a) {
1.497 + if (_comp[_gr.target(a)] == k) _out_arcs[n].push_back(a);
1.498 + }
1.499 + }
1.500 + }
1.501 + }
1.502 +
1.503 + // Initialize path data for the current component
1.504 + bool initComponent(int comp) {
1.505 + _nodes = &(_comp_nodes[comp]);
1.506 + int n = _nodes->size();
1.507 + if (n < 1 || (n == 1 && _out_arcs[(*_nodes)[0]].size() == 0)) {
1.508 + return false;
1.509 + }
1.510 + for (int i = 0; i < n; ++i) {
1.511 + _data[(*_nodes)[i]].resize(n + 1, PathData(INF));
1.512 + }
1.513 + return true;
1.514 + }
1.515 +
1.516 + // Process all rounds of computing path data for the current component.
1.517 + // _data[v][k] is the cost of a shortest directed walk from the root
1.518 + // node to node v containing exactly k arcs.
1.519 + void processRounds() {
1.520 + Node start = (*_nodes)[0];
1.521 + _data[start][0] = PathData(0);
1.522 + _process.clear();
1.523 + _process.push_back(start);
1.524 +
1.525 + int k, n = _nodes->size();
1.526 + int next_check = 4;
1.527 + bool terminate = false;
1.528 + for (k = 1; k <= n && int(_process.size()) < n && !terminate; ++k) {
1.529 + processNextBuildRound(k);
1.530 + if (k == next_check || k == n) {
1.531 + terminate = checkTermination(k);
1.532 + next_check = next_check * 3 / 2;
1.533 + }
1.534 + }
1.535 + for ( ; k <= n && !terminate; ++k) {
1.536 + processNextFullRound(k);
1.537 + if (k == next_check || k == n) {
1.538 + terminate = checkTermination(k);
1.539 + next_check = next_check * 3 / 2;
1.540 + }
1.541 + }
1.542 + }
1.543 +
1.544 + // Process one round and rebuild _process
1.545 + void processNextBuildRound(int k) {
1.546 + std::vector<Node> next;
1.547 + Node u, v;
1.548 + Arc e;
1.549 + LargeCost d;
1.550 + for (int i = 0; i < int(_process.size()); ++i) {
1.551 + u = _process[i];
1.552 + for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
1.553 + e = _out_arcs[u][j];
1.554 + v = _gr.target(e);
1.555 + d = _data[u][k-1].dist + _cost[e];
1.556 + if (_tolerance.less(d, _data[v][k].dist)) {
1.557 + if (_data[v][k].dist == INF) next.push_back(v);
1.558 + _data[v][k] = PathData(d, e);
1.559 + }
1.560 + }
1.561 + }
1.562 + _process.swap(next);
1.563 + }
1.564 +
1.565 + // Process one round using _nodes instead of _process
1.566 + void processNextFullRound(int k) {
1.567 + Node u, v;
1.568 + Arc e;
1.569 + LargeCost d;
1.570 + for (int i = 0; i < int(_nodes->size()); ++i) {
1.571 + u = (*_nodes)[i];
1.572 + for (int j = 0; j < int(_out_arcs[u].size()); ++j) {
1.573 + e = _out_arcs[u][j];
1.574 + v = _gr.target(e);
1.575 + d = _data[u][k-1].dist + _cost[e];
1.576 + if (_tolerance.less(d, _data[v][k].dist)) {
1.577 + _data[v][k] = PathData(d, e);
1.578 + }
1.579 + }
1.580 + }
1.581 + }
1.582 +
1.583 + // Check early termination
1.584 + bool checkTermination(int k) {
1.585 + typedef std::pair<int, int> Pair;
1.586 + typename GR::template NodeMap<Pair> level(_gr, Pair(-1, 0));
1.587 + typename GR::template NodeMap<LargeCost> pi(_gr);
1.588 + int n = _nodes->size();
1.589 + LargeCost cost;
1.590 + int size;
1.591 + Node u;
1.592 +
1.593 + // Search for cycles that are already found
1.594 + _curr_found = false;
1.595 + for (int i = 0; i < n; ++i) {
1.596 + u = (*_nodes)[i];
1.597 + if (_data[u][k].dist == INF) continue;
1.598 + for (int j = k; j >= 0; --j) {
1.599 + if (level[u].first == i && level[u].second > 0) {
1.600 + // A cycle is found
1.601 + cost = _data[u][level[u].second].dist - _data[u][j].dist;
1.602 + size = level[u].second - j;
1.603 + if (!_curr_found || cost * _curr_size < _curr_cost * size) {
1.604 + _curr_cost = cost;
1.605 + _curr_size = size;
1.606 + _curr_node = u;
1.607 + _curr_level = level[u].second;
1.608 + _curr_found = true;
1.609 + }
1.610 + }
1.611 + level[u] = Pair(i, j);
1.612 + if (j != 0) {
1.613 + u = _gr.source(_data[u][j].pred);
1.614 + }
1.615 + }
1.616 + }
1.617 +
1.618 + // If at least one cycle is found, check the optimality condition
1.619 + LargeCost d;
1.620 + if (_curr_found && k < n) {
1.621 + // Find node potentials
1.622 + for (int i = 0; i < n; ++i) {
1.623 + u = (*_nodes)[i];
1.624 + pi[u] = INF;
1.625 + for (int j = 0; j <= k; ++j) {
1.626 + if (_data[u][j].dist < INF) {
1.627 + d = _data[u][j].dist * _curr_size - j * _curr_cost;
1.628 + if (_tolerance.less(d, pi[u])) pi[u] = d;
1.629 + }
1.630 + }
1.631 + }
1.632 +
1.633 + // Check the optimality condition for all arcs
1.634 + bool done = true;
1.635 + for (ArcIt a(_gr); a != INVALID; ++a) {
1.636 + if (_tolerance.less(_cost[a] * _curr_size - _curr_cost,
1.637 + pi[_gr.target(a)] - pi[_gr.source(a)]) ) {
1.638 + done = false;
1.639 + break;
1.640 + }
1.641 + }
1.642 + return done;
1.643 + }
1.644 + return (k == n);
1.645 + }
1.646 +
1.647 + }; //class HartmannOrlinMmc
1.648 +
1.649 + ///@}
1.650 +
1.651 +} //namespace lemon
1.652 +
1.653 +#endif //LEMON_HARTMANN_ORLIN_MMC_H