Merge #438 and #436
authorAlpar Juttner <alpar@cs.elte.hu>
Fri, 22 Feb 2013 14:12:48 +0100
changeset 1014bc726f4892c7
parent 1011 9a6c9d4ee77b
parent 1013 f6f6896a4724
child 1016 97975184f4aa
Merge #438 and #436
     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;