1.1 --- a/lemon/cycle_canceling.h Wed Nov 28 12:08:47 2012 +0100
1.2 +++ b/lemon/cycle_canceling.h Fri Feb 22 14:12:48 2013 +0100
1.3 @@ -35,6 +35,7 @@
1.4 #include <lemon/circulation.h>
1.5 #include <lemon/bellman_ford.h>
1.6 #include <lemon/howard_mmc.h>
1.7 +#include <lemon/hartmann_orlin_mmc.h>
1.8
1.9 namespace lemon {
1.10
1.11 @@ -927,19 +928,42 @@
1.12
1.13 // Execute the "Minimum Mean Cycle Canceling" method
1.14 void startMinMeanCycleCanceling() {
1.15 - typedef SimplePath<StaticDigraph> SPath;
1.16 + typedef Path<StaticDigraph> SPath;
1.17 typedef typename SPath::ArcIt SPathArcIt;
1.18 typedef typename HowardMmc<StaticDigraph, CostArcMap>
1.19 - ::template SetPath<SPath>::Create MMC;
1.20 + ::template SetPath<SPath>::Create HwMmc;
1.21 + typedef typename HartmannOrlinMmc<StaticDigraph, CostArcMap>
1.22 + ::template SetPath<SPath>::Create HoMmc;
1.23 +
1.24 + const double HW_ITER_LIMIT_FACTOR = 1.0;
1.25 + const int HW_ITER_LIMIT_MIN_VALUE = 5;
1.26 +
1.27 + const int hw_iter_limit =
1.28 + std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num),
1.29 + HW_ITER_LIMIT_MIN_VALUE);
1.30
1.31 SPath cycle;
1.32 - MMC mmc(_sgr, _cost_map);
1.33 - mmc.cycle(cycle);
1.34 + HwMmc hw_mmc(_sgr, _cost_map);
1.35 + hw_mmc.cycle(cycle);
1.36 buildResidualNetwork();
1.37 - while (mmc.findCycleMean() && mmc.cycleCost() < 0) {
1.38 - // Find the cycle
1.39 - mmc.findCycle();
1.40 -
1.41 + while (true) {
1.42 +
1.43 + typename HwMmc::TerminationCause hw_tc =
1.44 + hw_mmc.findCycleMean(hw_iter_limit);
1.45 + if (hw_tc == HwMmc::ITERATION_LIMIT) {
1.46 + // Howard's algorithm reached the iteration limit, start a
1.47 + // strongly polynomial algorithm instead
1.48 + HoMmc ho_mmc(_sgr, _cost_map);
1.49 + ho_mmc.cycle(cycle);
1.50 + // Find a minimum mean cycle (Hartmann-Orlin algorithm)
1.51 + if (!(ho_mmc.findCycleMean() && ho_mmc.cycleCost() < 0)) break;
1.52 + ho_mmc.findCycle();
1.53 + } else {
1.54 + // Find a minimum mean cycle (Howard algorithm)
1.55 + if (!(hw_tc == HwMmc::OPTIMAL && hw_mmc.cycleCost() < 0)) break;
1.56 + hw_mmc.findCycle();
1.57 + }
1.58 +
1.59 // Compute delta value
1.60 Value delta = INF;
1.61 for (SPathArcIt a(cycle); a != INVALID; ++a) {
1.62 @@ -964,6 +988,12 @@
1.63 // Constants for the min mean cycle computations
1.64 const double LIMIT_FACTOR = 1.0;
1.65 const int MIN_LIMIT = 5;
1.66 + const double HW_ITER_LIMIT_FACTOR = 1.0;
1.67 + const int HW_ITER_LIMIT_MIN_VALUE = 5;
1.68 +
1.69 + const int hw_iter_limit =
1.70 + std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num),
1.71 + HW_ITER_LIMIT_MIN_VALUE);
1.72
1.73 // Contruct auxiliary data vectors
1.74 DoubleVector pi(_res_node_num, 0.0);
1.75 @@ -1137,17 +1167,30 @@
1.76 }
1.77 }
1.78 } else {
1.79 - typedef HowardMmc<StaticDigraph, CostArcMap> MMC;
1.80 + typedef HowardMmc<StaticDigraph, CostArcMap> HwMmc;
1.81 + typedef HartmannOrlinMmc<StaticDigraph, CostArcMap> HoMmc;
1.82 typedef typename BellmanFord<StaticDigraph, CostArcMap>
1.83 ::template SetDistMap<CostNodeMap>::Create BF;
1.84
1.85 // Set epsilon to the minimum cycle mean
1.86 + Cost cycle_cost = 0;
1.87 + int cycle_size = 1;
1.88 buildResidualNetwork();
1.89 - MMC mmc(_sgr, _cost_map);
1.90 - mmc.findCycleMean();
1.91 - epsilon = -mmc.cycleMean();
1.92 - Cost cycle_cost = mmc.cycleCost();
1.93 - int cycle_size = mmc.cycleSize();
1.94 + HwMmc hw_mmc(_sgr, _cost_map);
1.95 + if (hw_mmc.findCycleMean(hw_iter_limit) == HwMmc::ITERATION_LIMIT) {
1.96 + // Howard's algorithm reached the iteration limit, start a
1.97 + // strongly polynomial algorithm instead
1.98 + HoMmc ho_mmc(_sgr, _cost_map);
1.99 + ho_mmc.findCycleMean();
1.100 + epsilon = -ho_mmc.cycleMean();
1.101 + cycle_cost = ho_mmc.cycleCost();
1.102 + cycle_size = ho_mmc.cycleSize();
1.103 + } else {
1.104 + // Set epsilon
1.105 + epsilon = -hw_mmc.cycleMean();
1.106 + cycle_cost = hw_mmc.cycleCost();
1.107 + cycle_size = hw_mmc.cycleSize();
1.108 + }
1.109
1.110 // Compute feasible potentials for the current epsilon
1.111 for (int i = 0; i != int(_cost_vec.size()); ++i) {
2.1 --- a/lemon/howard_mmc.h Wed Nov 28 12:08:47 2012 +0100
2.2 +++ b/lemon/howard_mmc.h Fri Feb 22 14:12:48 2013 +0100
2.3 @@ -149,6 +149,23 @@
2.4 /// The \ref HowardMmcDefaultTraits "traits class" of the algorithm
2.5 typedef TR Traits;
2.6
2.7 + /// \brief Constants for the causes of search termination.
2.8 + ///
2.9 + /// Enum type containing constants for the different causes of search
2.10 + /// termination. The \ref findCycleMean() function returns one of
2.11 + /// these values.
2.12 + enum TerminationCause {
2.13 +
2.14 + /// No directed cycle can be found in the digraph.
2.15 + NO_CYCLE = 0,
2.16 +
2.17 + /// Optimal solution (minimum cycle mean) is found.
2.18 + OPTIMAL = 1,
2.19 +
2.20 + /// The iteration count limit is reached.
2.21 + ITERATION_LIMIT
2.22 + };
2.23 +
2.24 private:
2.25
2.26 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
2.27 @@ -324,25 +341,46 @@
2.28 return findCycleMean() && findCycle();
2.29 }
2.30
2.31 - /// \brief Find the minimum cycle mean.
2.32 + /// \brief Find the minimum cycle mean (or an upper bound).
2.33 ///
2.34 /// This function finds the minimum mean cost of the directed
2.35 - /// cycles in the digraph.
2.36 + /// cycles in the digraph (or an upper bound for it).
2.37 ///
2.38 - /// \return \c true if a directed cycle exists in the digraph.
2.39 - bool findCycleMean() {
2.40 + /// By default, the function finds the exact minimum cycle mean,
2.41 + /// but an optional limit can also be specified for the number of
2.42 + /// iterations performed during the search process.
2.43 + /// The return value indicates if the optimal solution is found
2.44 + /// or the iteration limit is reached. In the latter case, an
2.45 + /// approximate solution is provided, which corresponds to a directed
2.46 + /// cycle whose mean cost is relatively small, but not necessarily
2.47 + /// minimal.
2.48 + ///
2.49 + /// \param limit The maximum allowed number of iterations during
2.50 + /// the search process. Its default value implies that the algorithm
2.51 + /// runs until it finds the exact optimal solution.
2.52 + ///
2.53 + /// \return The termination cause of the search process.
2.54 + /// For more information, see \ref TerminationCause.
2.55 + TerminationCause findCycleMean(int limit = std::numeric_limits<int>::max()) {
2.56 // Initialize and find strongly connected components
2.57 init();
2.58 findComponents();
2.59
2.60 // Find the minimum cycle mean in the components
2.61 + int iter_count = 0;
2.62 + bool iter_limit_reached = false;
2.63 for (int comp = 0; comp < _comp_num; ++comp) {
2.64 // Find the minimum mean cycle in the current component
2.65 if (!buildPolicyGraph(comp)) continue;
2.66 while (true) {
2.67 + if (++iter_count > limit) {
2.68 + iter_limit_reached = true;
2.69 + break;
2.70 + }
2.71 findPolicyCycle();
2.72 if (!computeNodeDistances()) break;
2.73 }
2.74 +
2.75 // Update the best cycle (global minimum mean cycle)
2.76 if ( _curr_found && (!_best_found ||
2.77 _curr_cost * _best_size < _best_cost * _curr_size) ) {
2.78 @@ -351,8 +389,15 @@
2.79 _best_size = _curr_size;
2.80 _best_node = _curr_node;
2.81 }
2.82 +
2.83 + if (iter_limit_reached) break;
2.84 }
2.85 - return _best_found;
2.86 +
2.87 + if (iter_limit_reached) {
2.88 + return ITERATION_LIMIT;
2.89 + } else {
2.90 + return _best_found ? OPTIMAL : NO_CYCLE;
2.91 + }
2.92 }
2.93
2.94 /// \brief Find a minimum mean directed cycle.
3.1 --- a/test/min_mean_cycle_test.cc Wed Nov 28 12:08:47 2012 +0100
3.2 +++ b/test/min_mean_cycle_test.cc Fri Feb 22 14:12:48 2013 +0100
3.3 @@ -110,7 +110,7 @@
3.4 const SmartDigraph::ArcMap<int>& cm,
3.5 int cost, int size) {
3.6 MMC alg(gr, lm);
3.7 - alg.findCycleMean();
3.8 + check(alg.findCycleMean(), "Wrong result");
3.9 check(alg.cycleMean() == static_cast<double>(cost) / size,
3.10 "Wrong cycle mean");
3.11 alg.findCycle();
3.12 @@ -210,6 +210,13 @@
3.13 checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l2, c2, 5, 2);
3.14 checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l3, c3, 0, 1);
3.15 checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l4, c4, -1, 1);
3.16 +
3.17 + // Howard with iteration limit
3.18 + HowardMmc<GR, IntArcMap> mmc(gr, l1);
3.19 + check((mmc.findCycleMean(2) == HowardMmc<GR, IntArcMap>::ITERATION_LIMIT),
3.20 + "Wrong termination cause");
3.21 + check((mmc.findCycleMean(4) == HowardMmc<GR, IntArcMap>::OPTIMAL),
3.22 + "Wrong termination cause");
3.23 }
3.24
3.25 return 0;