1.1 --- a/lemon/cycle_canceling.h	Thu Nov 15 07:05:29 2012 +0100
     1.2 +++ b/lemon/cycle_canceling.h	Thu Nov 15 07:17:48 2012 +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) {