alpar@956: /* -*- mode: C++; indent-tabs-mode: nil; -*- kpeter@805: * alpar@956: * This file is a part of LEMON, a generic C++ optimization library. kpeter@805: * alpar@956: * Copyright (C) 2003-2010 kpeter@805: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport kpeter@805: * (Egervary Research Group on Combinatorial Optimization, EGRES). kpeter@805: * kpeter@805: * Permission to use, modify and distribute this software is granted kpeter@805: * provided that this copyright notice appears in all copies. For kpeter@805: * precise terms see the accompanying LICENSE file. kpeter@805: * kpeter@805: * This software is provided "AS IS" with no warranty of any kind, kpeter@805: * express or implied, and with no claim as to its suitability for any kpeter@805: * purpose. kpeter@805: * kpeter@805: */ kpeter@805: kpeter@942: #ifndef LEMON_HOWARD_MMC_H kpeter@942: #define LEMON_HOWARD_MMC_H kpeter@805: kpeter@815: /// \ingroup min_mean_cycle kpeter@805: /// kpeter@805: /// \file kpeter@805: /// \brief Howard's algorithm for finding a minimum mean cycle. kpeter@805: kpeter@805: #include kpeter@810: #include kpeter@805: #include kpeter@805: #include kpeter@805: #include kpeter@805: #include kpeter@805: kpeter@805: namespace lemon { kpeter@805: kpeter@942: /// \brief Default traits class of HowardMmc class. kpeter@808: /// kpeter@942: /// Default traits class of HowardMmc class. kpeter@808: /// \tparam GR The type of the digraph. kpeter@942: /// \tparam CM The type of the cost map. kpeter@808: /// It must conform to the \ref concepts::ReadMap "ReadMap" concept. kpeter@808: #ifdef DOXYGEN kpeter@942: template kpeter@808: #else kpeter@942: template ::is_integer> kpeter@808: #endif kpeter@942: struct HowardMmcDefaultTraits kpeter@808: { kpeter@808: /// The type of the digraph kpeter@808: typedef GR Digraph; kpeter@942: /// The type of the cost map kpeter@942: typedef CM CostMap; kpeter@942: /// The type of the arc costs kpeter@942: typedef typename CostMap::Value Cost; kpeter@808: kpeter@942: /// \brief The large cost type used for internal computations kpeter@808: /// kpeter@942: /// The large cost type used for internal computations. kpeter@942: /// It is \c long \c long if the \c Cost type is integer, kpeter@808: /// otherwise it is \c double. kpeter@942: /// \c Cost must be convertible to \c LargeCost. kpeter@942: typedef double LargeCost; kpeter@808: kpeter@808: /// The tolerance type used for internal computations kpeter@942: typedef lemon::Tolerance Tolerance; kpeter@808: kpeter@808: /// \brief The path type of the found cycles kpeter@808: /// kpeter@808: /// The path type of the found cycles. kpeter@808: /// It must conform to the \ref lemon::concepts::Path "Path" concept kpeter@808: /// and it must have an \c addBack() function. kpeter@808: typedef lemon::Path Path; kpeter@808: }; kpeter@808: kpeter@942: // Default traits class for integer cost types kpeter@942: template kpeter@942: struct HowardMmcDefaultTraits kpeter@808: { kpeter@808: typedef GR Digraph; kpeter@942: typedef CM CostMap; kpeter@942: typedef typename CostMap::Value Cost; kpeter@808: #ifdef LEMON_HAVE_LONG_LONG kpeter@942: typedef long long LargeCost; kpeter@808: #else kpeter@942: typedef long LargeCost; kpeter@808: #endif kpeter@942: typedef lemon::Tolerance Tolerance; kpeter@808: typedef lemon::Path Path; kpeter@808: }; kpeter@808: kpeter@808: kpeter@815: /// \addtogroup min_mean_cycle kpeter@805: /// @{ kpeter@805: kpeter@805: /// \brief Implementation of Howard's algorithm for finding a minimum kpeter@805: /// mean cycle. kpeter@805: /// kpeter@811: /// This class implements Howard's policy iteration algorithm for finding kpeter@942: /// a directed cycle of minimum mean cost in a digraph kpeter@1164: /// \ref dasdan98minmeancycle, \ref dasdan04experimental. kpeter@815: /// This class provides the most efficient algorithm for the kpeter@815: /// minimum mean cycle problem, though the best known theoretical kpeter@815: /// bound on its running time is exponential. kpeter@805: /// kpeter@805: /// \tparam GR The type of the digraph the algorithm runs on. kpeter@942: /// \tparam CM The type of the cost map. The default kpeter@805: /// map type is \ref concepts::Digraph::ArcMap "GR::ArcMap". kpeter@891: /// \tparam TR The traits class that defines various types used by the kpeter@942: /// algorithm. By default, it is \ref HowardMmcDefaultTraits kpeter@942: /// "HowardMmcDefaultTraits". kpeter@891: /// In most cases, this parameter should not be set directly, kpeter@891: /// consider to use the named template parameters instead. kpeter@805: #ifdef DOXYGEN kpeter@942: template kpeter@805: #else kpeter@805: template < typename GR, kpeter@942: typename CM = typename GR::template ArcMap, kpeter@942: typename TR = HowardMmcDefaultTraits > kpeter@805: #endif kpeter@942: class HowardMmc kpeter@805: { kpeter@805: public: alpar@956: kpeter@808: /// The type of the digraph kpeter@808: typedef typename TR::Digraph Digraph; kpeter@942: /// The type of the cost map kpeter@942: typedef typename TR::CostMap CostMap; kpeter@942: /// The type of the arc costs kpeter@942: typedef typename TR::Cost Cost; kpeter@808: kpeter@942: /// \brief The large cost type kpeter@808: /// kpeter@942: /// The large cost type used for internal computations. kpeter@942: /// By default, it is \c long \c long if the \c Cost type is integer, kpeter@808: /// otherwise it is \c double. kpeter@942: typedef typename TR::LargeCost LargeCost; kpeter@808: kpeter@808: /// The tolerance type kpeter@808: typedef typename TR::Tolerance Tolerance; kpeter@808: kpeter@808: /// \brief The path type of the found cycles kpeter@808: /// kpeter@808: /// The path type of the found cycles. kpeter@942: /// Using the \ref HowardMmcDefaultTraits "default traits class", kpeter@808: /// it is \ref lemon::Path "Path". kpeter@808: typedef typename TR::Path Path; kpeter@808: kpeter@942: /// The \ref HowardMmcDefaultTraits "traits class" of the algorithm kpeter@808: typedef TR Traits; kpeter@805: kpeter@1178: /// \brief Constants for the causes of search termination. kpeter@1178: /// kpeter@1178: /// Enum type containing constants for the different causes of search kpeter@1178: /// termination. The \ref findCycleMean() function returns one of kpeter@1178: /// these values. kpeter@1178: enum TerminationCause { kpeter@1178: kpeter@1178: /// No directed cycle can be found in the digraph. kpeter@1178: NO_CYCLE = 0, kpeter@1178: kpeter@1178: /// Optimal solution (minimum cycle mean) is found. kpeter@1178: OPTIMAL = 1, kpeter@1178: kpeter@1178: /// The iteration count limit is reached. kpeter@1178: ITERATION_LIMIT kpeter@1178: }; kpeter@1178: kpeter@805: private: kpeter@805: kpeter@805: TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); alpar@956: kpeter@805: // The digraph the algorithm runs on kpeter@805: const Digraph &_gr; kpeter@942: // The cost of the arcs kpeter@942: const CostMap &_cost; kpeter@805: kpeter@807: // Data for the found cycles kpeter@807: bool _curr_found, _best_found; kpeter@942: LargeCost _curr_cost, _best_cost; kpeter@807: int _curr_size, _best_size; kpeter@807: Node _curr_node, _best_node; kpeter@807: kpeter@805: Path *_cycle_path; kpeter@807: bool _local_path; kpeter@805: kpeter@807: // Internal data used by the algorithm kpeter@807: typename Digraph::template NodeMap _policy; kpeter@807: typename Digraph::template NodeMap _reached; kpeter@807: typename Digraph::template NodeMap _level; kpeter@942: typename Digraph::template NodeMap _dist; kpeter@805: kpeter@807: // Data for storing the strongly connected components kpeter@807: int _comp_num; kpeter@805: typename Digraph::template NodeMap _comp; kpeter@807: std::vector > _comp_nodes; kpeter@807: std::vector* _nodes; kpeter@807: typename Digraph::template NodeMap > _in_arcs; alpar@956: kpeter@807: // Queue used for BFS search kpeter@807: std::vector _queue; kpeter@807: int _qfront, _qback; kpeter@808: kpeter@808: Tolerance _tolerance; alpar@956: kpeter@814: // Infinite constant kpeter@942: const LargeCost INF; kpeter@814: kpeter@808: public: alpar@956: kpeter@808: /// \name Named Template Parameters kpeter@808: /// @{ kpeter@808: kpeter@808: template kpeter@942: struct SetLargeCostTraits : public Traits { kpeter@942: typedef T LargeCost; kpeter@808: typedef lemon::Tolerance Tolerance; kpeter@808: }; kpeter@808: kpeter@808: /// \brief \ref named-templ-param "Named parameter" for setting kpeter@942: /// \c LargeCost type. kpeter@808: /// kpeter@942: /// \ref named-templ-param "Named parameter" for setting \c LargeCost kpeter@808: /// type. It is used for internal computations in the algorithm. kpeter@808: template kpeter@942: struct SetLargeCost kpeter@942: : public HowardMmc > { kpeter@942: typedef HowardMmc > Create; kpeter@808: }; kpeter@808: kpeter@808: template kpeter@808: struct SetPathTraits : public Traits { kpeter@808: typedef T Path; kpeter@808: }; kpeter@808: kpeter@808: /// \brief \ref named-templ-param "Named parameter" for setting kpeter@808: /// \c %Path type. kpeter@808: /// kpeter@808: /// \ref named-templ-param "Named parameter" for setting the \c %Path kpeter@808: /// type of the found cycles. kpeter@808: /// It must conform to the \ref lemon::concepts::Path "Path" concept kpeter@808: /// and it must have an \c addBack() function. kpeter@808: template kpeter@808: struct SetPath kpeter@942: : public HowardMmc > { kpeter@942: typedef HowardMmc > Create; kpeter@808: }; alpar@956: kpeter@808: /// @} kpeter@805: kpeter@941: protected: kpeter@941: kpeter@942: HowardMmc() {} kpeter@941: kpeter@805: public: kpeter@805: kpeter@805: /// \brief Constructor. kpeter@805: /// kpeter@805: /// The constructor of the class. kpeter@805: /// kpeter@805: /// \param digraph The digraph the algorithm runs on. kpeter@942: /// \param cost The costs of the arcs. kpeter@942: HowardMmc( const Digraph &digraph, kpeter@942: const CostMap &cost ) : kpeter@942: _gr(digraph), _cost(cost), _best_found(false), kpeter@942: _best_cost(0), _best_size(1), _cycle_path(NULL), _local_path(false), kpeter@807: _policy(digraph), _reached(digraph), _level(digraph), _dist(digraph), kpeter@814: _comp(digraph), _in_arcs(digraph), kpeter@942: INF(std::numeric_limits::has_infinity ? kpeter@942: std::numeric_limits::infinity() : kpeter@942: std::numeric_limits::max()) kpeter@805: {} kpeter@805: kpeter@805: /// Destructor. kpeter@942: ~HowardMmc() { kpeter@805: if (_local_path) delete _cycle_path; kpeter@805: } kpeter@805: kpeter@805: /// \brief Set the path structure for storing the found cycle. kpeter@805: /// kpeter@805: /// This function sets an external path structure for storing the kpeter@805: /// found cycle. kpeter@805: /// kpeter@805: /// If you don't call this function before calling \ref run() or kpeter@942: /// \ref findCycleMean(), it will allocate a local \ref Path "path" kpeter@805: /// structure. The destuctor deallocates this automatically kpeter@805: /// allocated object, of course. kpeter@805: /// kpeter@805: /// \note The algorithm calls only the \ref lemon::Path::addBack() kpeter@805: /// "addBack()" function of the given path structure. kpeter@805: /// kpeter@805: /// \return (*this) kpeter@942: HowardMmc& cycle(Path &path) { kpeter@805: if (_local_path) { kpeter@805: delete _cycle_path; kpeter@805: _local_path = false; kpeter@805: } kpeter@805: _cycle_path = &path; kpeter@805: return *this; kpeter@805: } kpeter@805: kpeter@816: /// \brief Set the tolerance used by the algorithm. kpeter@816: /// kpeter@816: /// This function sets the tolerance object used by the algorithm. kpeter@816: /// kpeter@816: /// \return (*this) kpeter@942: HowardMmc& tolerance(const Tolerance& tolerance) { kpeter@816: _tolerance = tolerance; kpeter@816: return *this; kpeter@816: } kpeter@816: kpeter@816: /// \brief Return a const reference to the tolerance. kpeter@816: /// kpeter@816: /// This function returns a const reference to the tolerance object kpeter@816: /// used by the algorithm. kpeter@816: const Tolerance& tolerance() const { kpeter@816: return _tolerance; kpeter@816: } kpeter@816: kpeter@805: /// \name Execution control kpeter@805: /// The simplest way to execute the algorithm is to call the \ref run() kpeter@805: /// function.\n kpeter@942: /// If you only need the minimum mean cost, you may call kpeter@942: /// \ref findCycleMean(). kpeter@805: kpeter@805: /// @{ kpeter@805: kpeter@805: /// \brief Run the algorithm. kpeter@805: /// kpeter@805: /// This function runs the algorithm. kpeter@806: /// It can be called more than once (e.g. if the underlying digraph kpeter@942: /// and/or the arc costs have been modified). kpeter@805: /// kpeter@805: /// \return \c true if a directed cycle exists in the digraph. kpeter@805: /// kpeter@806: /// \note mmc.run() is just a shortcut of the following code. kpeter@805: /// \code kpeter@942: /// return mmc.findCycleMean() && mmc.findCycle(); kpeter@805: /// \endcode kpeter@805: bool run() { kpeter@942: return findCycleMean() && findCycle(); kpeter@805: } kpeter@805: kpeter@1178: /// \brief Find the minimum cycle mean (or an upper bound). kpeter@805: /// kpeter@942: /// This function finds the minimum mean cost of the directed kpeter@1178: /// cycles in the digraph (or an upper bound for it). kpeter@805: /// kpeter@1178: /// By default, the function finds the exact minimum cycle mean, kpeter@1178: /// but an optional limit can also be specified for the number of kpeter@1178: /// iterations performed during the search process. kpeter@1178: /// The return value indicates if the optimal solution is found kpeter@1178: /// or the iteration limit is reached. In the latter case, an kpeter@1178: /// approximate solution is provided, which corresponds to a directed kpeter@1178: /// cycle whose mean cost is relatively small, but not necessarily kpeter@1178: /// minimal. kpeter@1178: /// kpeter@1178: /// \param limit The maximum allowed number of iterations during kpeter@1178: /// the search process. Its default value implies that the algorithm kpeter@1178: /// runs until it finds the exact optimal solution. kpeter@1178: /// kpeter@1178: /// \return The termination cause of the search process. kpeter@1178: /// For more information, see \ref TerminationCause. kpeter@1178: TerminationCause findCycleMean(int limit = std::numeric_limits::max()) { kpeter@807: // Initialize and find strongly connected components kpeter@807: init(); kpeter@807: findComponents(); alpar@956: kpeter@806: // Find the minimum cycle mean in the components kpeter@1178: int iter_count = 0; kpeter@1178: bool iter_limit_reached = false; kpeter@805: for (int comp = 0; comp < _comp_num; ++comp) { kpeter@807: // Find the minimum mean cycle in the current component kpeter@807: if (!buildPolicyGraph(comp)) continue; kpeter@805: while (true) { kpeter@1178: if (++iter_count > limit) { kpeter@1178: iter_limit_reached = true; kpeter@1178: break; kpeter@1178: } kpeter@807: findPolicyCycle(); kpeter@805: if (!computeNodeDistances()) break; kpeter@805: } kpeter@1178: kpeter@807: // Update the best cycle (global minimum mean cycle) kpeter@814: if ( _curr_found && (!_best_found || kpeter@942: _curr_cost * _best_size < _best_cost * _curr_size) ) { kpeter@807: _best_found = true; kpeter@942: _best_cost = _curr_cost; kpeter@807: _best_size = _curr_size; kpeter@807: _best_node = _curr_node; kpeter@807: } kpeter@1178: kpeter@1178: if (iter_limit_reached) break; kpeter@805: } kpeter@1178: kpeter@1178: if (iter_limit_reached) { kpeter@1178: return ITERATION_LIMIT; kpeter@1178: } else { kpeter@1178: return _best_found ? OPTIMAL : NO_CYCLE; kpeter@1178: } kpeter@805: } kpeter@805: kpeter@805: /// \brief Find a minimum mean directed cycle. kpeter@805: /// kpeter@942: /// This function finds a directed cycle of minimum mean cost kpeter@942: /// in the digraph using the data computed by findCycleMean(). kpeter@805: /// kpeter@805: /// \return \c true if a directed cycle exists in the digraph. kpeter@805: /// kpeter@942: /// \pre \ref findCycleMean() must be called before using this function. kpeter@805: bool findCycle() { kpeter@807: if (!_best_found) return false; kpeter@807: _cycle_path->addBack(_policy[_best_node]); kpeter@807: for ( Node v = _best_node; kpeter@807: (v = _gr.target(_policy[v])) != _best_node; ) { kpeter@805: _cycle_path->addBack(_policy[v]); kpeter@805: } kpeter@805: return true; kpeter@805: } kpeter@805: kpeter@805: /// @} kpeter@805: kpeter@805: /// \name Query Functions kpeter@806: /// The results of the algorithm can be obtained using these kpeter@805: /// functions.\n kpeter@805: /// The algorithm should be executed before using them. kpeter@805: kpeter@805: /// @{ kpeter@805: kpeter@942: /// \brief Return the total cost of the found cycle. kpeter@805: /// kpeter@942: /// This function returns the total cost of the found cycle. kpeter@805: /// kpeter@942: /// \pre \ref run() or \ref findCycleMean() must be called before kpeter@805: /// using this function. kpeter@942: Cost cycleCost() const { kpeter@942: return static_cast(_best_cost); kpeter@805: } kpeter@805: kpeter@805: /// \brief Return the number of arcs on the found cycle. kpeter@805: /// kpeter@805: /// This function returns the number of arcs on the found cycle. kpeter@805: /// kpeter@942: /// \pre \ref run() or \ref findCycleMean() must be called before kpeter@805: /// using this function. kpeter@942: int cycleSize() const { kpeter@807: return _best_size; kpeter@805: } kpeter@805: kpeter@942: /// \brief Return the mean cost of the found cycle. kpeter@805: /// kpeter@942: /// This function returns the mean cost of the found cycle. kpeter@805: /// kpeter@807: /// \note alg.cycleMean() is just a shortcut of the kpeter@805: /// following code. kpeter@805: /// \code kpeter@942: /// return static_cast(alg.cycleCost()) / alg.cycleSize(); kpeter@805: /// \endcode kpeter@805: /// kpeter@942: /// \pre \ref run() or \ref findCycleMean() must be called before kpeter@805: /// using this function. kpeter@805: double cycleMean() const { kpeter@942: return static_cast(_best_cost) / _best_size; kpeter@805: } kpeter@805: kpeter@805: /// \brief Return the found cycle. kpeter@805: /// kpeter@805: /// This function returns a const reference to the path structure kpeter@805: /// storing the found cycle. kpeter@805: /// kpeter@805: /// \pre \ref run() or \ref findCycle() must be called before using kpeter@805: /// this function. kpeter@805: const Path& cycle() const { kpeter@805: return *_cycle_path; kpeter@805: } kpeter@805: kpeter@805: ///@} kpeter@805: kpeter@805: private: kpeter@805: kpeter@807: // Initialize kpeter@807: void init() { kpeter@807: if (!_cycle_path) { kpeter@807: _local_path = true; kpeter@807: _cycle_path = new Path; kpeter@805: } kpeter@807: _queue.resize(countNodes(_gr)); kpeter@807: _best_found = false; kpeter@942: _best_cost = 0; kpeter@807: _best_size = 1; kpeter@807: _cycle_path->clear(); kpeter@807: } alpar@956: kpeter@807: // Find strongly connected components and initialize _comp_nodes kpeter@807: // and _in_arcs kpeter@807: void findComponents() { kpeter@807: _comp_num = stronglyConnectedComponents(_gr, _comp); kpeter@807: _comp_nodes.resize(_comp_num); kpeter@807: if (_comp_num == 1) { kpeter@807: _comp_nodes[0].clear(); kpeter@807: for (NodeIt n(_gr); n != INVALID; ++n) { kpeter@807: _comp_nodes[0].push_back(n); kpeter@807: _in_arcs[n].clear(); kpeter@807: for (InArcIt a(_gr, n); a != INVALID; ++a) { kpeter@807: _in_arcs[n].push_back(a); kpeter@807: } kpeter@807: } kpeter@807: } else { kpeter@807: for (int i = 0; i < _comp_num; ++i) kpeter@807: _comp_nodes[i].clear(); kpeter@807: for (NodeIt n(_gr); n != INVALID; ++n) { kpeter@807: int k = _comp[n]; kpeter@807: _comp_nodes[k].push_back(n); kpeter@807: _in_arcs[n].clear(); kpeter@807: for (InArcIt a(_gr, n); a != INVALID; ++a) { kpeter@807: if (_comp[_gr.source(a)] == k) _in_arcs[n].push_back(a); kpeter@807: } kpeter@807: } kpeter@805: } kpeter@807: } kpeter@807: kpeter@807: // Build the policy graph in the given strongly connected component kpeter@807: // (the out-degree of every node is 1) kpeter@807: bool buildPolicyGraph(int comp) { kpeter@807: _nodes = &(_comp_nodes[comp]); kpeter@807: if (_nodes->size() < 1 || kpeter@807: (_nodes->size() == 1 && _in_arcs[(*_nodes)[0]].size() == 0)) { kpeter@807: return false; kpeter@805: } kpeter@807: for (int i = 0; i < int(_nodes->size()); ++i) { kpeter@814: _dist[(*_nodes)[i]] = INF; kpeter@807: } kpeter@807: Node u, v; kpeter@807: Arc e; kpeter@807: for (int i = 0; i < int(_nodes->size()); ++i) { kpeter@807: v = (*_nodes)[i]; kpeter@807: for (int j = 0; j < int(_in_arcs[v].size()); ++j) { kpeter@807: e = _in_arcs[v][j]; kpeter@807: u = _gr.source(e); kpeter@942: if (_cost[e] < _dist[u]) { kpeter@942: _dist[u] = _cost[e]; kpeter@807: _policy[u] = e; kpeter@807: } kpeter@805: } kpeter@805: } kpeter@805: return true; kpeter@805: } kpeter@805: kpeter@807: // Find the minimum mean cycle in the policy graph kpeter@807: void findPolicyCycle() { kpeter@807: for (int i = 0; i < int(_nodes->size()); ++i) { kpeter@807: _level[(*_nodes)[i]] = -1; kpeter@807: } kpeter@942: LargeCost ccost; kpeter@805: int csize; kpeter@805: Node u, v; kpeter@807: _curr_found = false; kpeter@807: for (int i = 0; i < int(_nodes->size()); ++i) { kpeter@807: u = (*_nodes)[i]; kpeter@807: if (_level[u] >= 0) continue; kpeter@807: for (; _level[u] < 0; u = _gr.target(_policy[u])) { kpeter@807: _level[u] = i; kpeter@807: } kpeter@807: if (_level[u] == i) { kpeter@807: // A cycle is found kpeter@942: ccost = _cost[_policy[u]]; kpeter@807: csize = 1; kpeter@807: for (v = u; (v = _gr.target(_policy[v])) != u; ) { kpeter@942: ccost += _cost[_policy[v]]; kpeter@807: ++csize; kpeter@805: } kpeter@807: if ( !_curr_found || kpeter@942: (ccost * _curr_size < _curr_cost * csize) ) { kpeter@807: _curr_found = true; kpeter@942: _curr_cost = ccost; kpeter@807: _curr_size = csize; kpeter@807: _curr_node = u; kpeter@805: } kpeter@805: } kpeter@805: } kpeter@805: } kpeter@805: kpeter@807: // Contract the policy graph and compute node distances kpeter@805: bool computeNodeDistances() { kpeter@807: // Find the component of the main cycle and compute node distances kpeter@807: // using reverse BFS kpeter@807: for (int i = 0; i < int(_nodes->size()); ++i) { kpeter@807: _reached[(*_nodes)[i]] = false; kpeter@807: } kpeter@807: _qfront = _qback = 0; kpeter@807: _queue[0] = _curr_node; kpeter@807: _reached[_curr_node] = true; kpeter@807: _dist[_curr_node] = 0; kpeter@805: Node u, v; kpeter@807: Arc e; kpeter@807: while (_qfront <= _qback) { kpeter@807: v = _queue[_qfront++]; kpeter@807: for (int j = 0; j < int(_in_arcs[v].size()); ++j) { kpeter@807: e = _in_arcs[v][j]; kpeter@805: u = _gr.source(e); kpeter@807: if (_policy[u] == e && !_reached[u]) { kpeter@807: _reached[u] = true; kpeter@942: _dist[u] = _dist[v] + _cost[e] * _curr_size - _curr_cost; kpeter@807: _queue[++_qback] = u; kpeter@805: } kpeter@805: } kpeter@805: } kpeter@807: kpeter@807: // Connect all other nodes to this component and compute node kpeter@807: // distances using reverse BFS kpeter@807: _qfront = 0; kpeter@807: while (_qback < int(_nodes->size())-1) { kpeter@807: v = _queue[_qfront++]; kpeter@807: for (int j = 0; j < int(_in_arcs[v].size()); ++j) { kpeter@807: e = _in_arcs[v][j]; kpeter@807: u = _gr.source(e); kpeter@807: if (!_reached[u]) { kpeter@807: _reached[u] = true; kpeter@807: _policy[u] = e; kpeter@942: _dist[u] = _dist[v] + _cost[e] * _curr_size - _curr_cost; kpeter@807: _queue[++_qback] = u; kpeter@807: } kpeter@807: } kpeter@807: } kpeter@807: kpeter@807: // Improve node distances kpeter@805: bool improved = false; kpeter@807: for (int i = 0; i < int(_nodes->size()); ++i) { kpeter@807: v = (*_nodes)[i]; kpeter@807: for (int j = 0; j < int(_in_arcs[v].size()); ++j) { kpeter@807: e = _in_arcs[v][j]; kpeter@807: u = _gr.source(e); kpeter@942: LargeCost delta = _dist[v] + _cost[e] * _curr_size - _curr_cost; kpeter@808: if (_tolerance.less(delta, _dist[u])) { kpeter@807: _dist[u] = delta; kpeter@807: _policy[u] = e; kpeter@807: improved = true; kpeter@807: } kpeter@805: } kpeter@805: } kpeter@805: return improved; kpeter@805: } kpeter@805: kpeter@942: }; //class HowardMmc kpeter@805: kpeter@805: ///@} kpeter@805: kpeter@805: } //namespace lemon kpeter@805: kpeter@942: #endif //LEMON_HOWARD_MMC_H