# HG changeset patch
# User Alpar Juttner <alpar@cs.elte.hu>
# Date 1361538768 -3600
# Node ID bc726f4892c723d27e6bc8b515979f24aecc13b8
# Parent  9a6c9d4ee77b8b5f529e525e2539b9317faecdd0# Parent  f6f6896a4724480ece9a55cc7387d328ba8dabfb
Merge #438 and #436

diff -r 9a6c9d4ee77b -r bc726f4892c7 lemon/cycle_canceling.h
--- a/lemon/cycle_canceling.h	Wed Nov 28 12:08:47 2012 +0100
+++ b/lemon/cycle_canceling.h	Fri Feb 22 14:12:48 2013 +0100
@@ -35,6 +35,7 @@
 #include <lemon/circulation.h>
 #include <lemon/bellman_ford.h>
 #include <lemon/howard_mmc.h>
+#include <lemon/hartmann_orlin_mmc.h>
 
 namespace lemon {
 
@@ -927,19 +928,42 @@
 
     // Execute the "Minimum Mean Cycle Canceling" method
     void startMinMeanCycleCanceling() {
-      typedef SimplePath<StaticDigraph> SPath;
+      typedef Path<StaticDigraph> SPath;
       typedef typename SPath::ArcIt SPathArcIt;
       typedef typename HowardMmc<StaticDigraph, CostArcMap>
-        ::template SetPath<SPath>::Create MMC;
+        ::template SetPath<SPath>::Create HwMmc;
+      typedef typename HartmannOrlinMmc<StaticDigraph, CostArcMap>
+        ::template SetPath<SPath>::Create HoMmc;
+
+      const double HW_ITER_LIMIT_FACTOR = 1.0;
+      const int HW_ITER_LIMIT_MIN_VALUE = 5;
+
+      const int hw_iter_limit =
+          std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num),
+                   HW_ITER_LIMIT_MIN_VALUE);
 
       SPath cycle;
-      MMC mmc(_sgr, _cost_map);
-      mmc.cycle(cycle);
+      HwMmc hw_mmc(_sgr, _cost_map);
+      hw_mmc.cycle(cycle);
       buildResidualNetwork();
-      while (mmc.findCycleMean() && mmc.cycleCost() < 0) {
-        // Find the cycle
-        mmc.findCycle();
-
+      while (true) {
+        
+        typename HwMmc::TerminationCause hw_tc =
+            hw_mmc.findCycleMean(hw_iter_limit);
+        if (hw_tc == HwMmc::ITERATION_LIMIT) {
+          // Howard's algorithm reached the iteration limit, start a
+          // strongly polynomial algorithm instead
+          HoMmc ho_mmc(_sgr, _cost_map);
+          ho_mmc.cycle(cycle);
+          // Find a minimum mean cycle (Hartmann-Orlin algorithm)
+          if (!(ho_mmc.findCycleMean() && ho_mmc.cycleCost() < 0)) break;
+          ho_mmc.findCycle();
+        } else {
+          // Find a minimum mean cycle (Howard algorithm)
+          if (!(hw_tc == HwMmc::OPTIMAL && hw_mmc.cycleCost() < 0)) break;
+          hw_mmc.findCycle();
+        }
+        
         // Compute delta value
         Value delta = INF;
         for (SPathArcIt a(cycle); a != INVALID; ++a) {
@@ -964,6 +988,12 @@
       // Constants for the min mean cycle computations
       const double LIMIT_FACTOR = 1.0;
       const int MIN_LIMIT = 5;
+      const double HW_ITER_LIMIT_FACTOR = 1.0;
+      const int HW_ITER_LIMIT_MIN_VALUE = 5;
+
+      const int hw_iter_limit =
+          std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num),
+                   HW_ITER_LIMIT_MIN_VALUE);
 
       // Contruct auxiliary data vectors
       DoubleVector pi(_res_node_num, 0.0);
@@ -1137,17 +1167,30 @@
             }
           }
         } else {
-          typedef HowardMmc<StaticDigraph, CostArcMap> MMC;
+          typedef HowardMmc<StaticDigraph, CostArcMap> HwMmc;
+          typedef HartmannOrlinMmc<StaticDigraph, CostArcMap> HoMmc;
           typedef typename BellmanFord<StaticDigraph, CostArcMap>
             ::template SetDistMap<CostNodeMap>::Create BF;
 
           // Set epsilon to the minimum cycle mean
+          Cost cycle_cost = 0;
+          int cycle_size = 1;
           buildResidualNetwork();
-          MMC mmc(_sgr, _cost_map);
-          mmc.findCycleMean();
-          epsilon = -mmc.cycleMean();
-          Cost cycle_cost = mmc.cycleCost();
-          int cycle_size = mmc.cycleSize();
+          HwMmc hw_mmc(_sgr, _cost_map);
+          if (hw_mmc.findCycleMean(hw_iter_limit) == HwMmc::ITERATION_LIMIT) {
+            // Howard's algorithm reached the iteration limit, start a
+            // strongly polynomial algorithm instead
+            HoMmc ho_mmc(_sgr, _cost_map);
+            ho_mmc.findCycleMean();
+            epsilon = -ho_mmc.cycleMean();
+            cycle_cost = ho_mmc.cycleCost();
+            cycle_size = ho_mmc.cycleSize();
+          } else {
+            // Set epsilon
+            epsilon = -hw_mmc.cycleMean();
+            cycle_cost = hw_mmc.cycleCost();
+            cycle_size = hw_mmc.cycleSize();
+          }
 
           // Compute feasible potentials for the current epsilon
           for (int i = 0; i != int(_cost_vec.size()); ++i) {
diff -r 9a6c9d4ee77b -r bc726f4892c7 lemon/howard_mmc.h
--- a/lemon/howard_mmc.h	Wed Nov 28 12:08:47 2012 +0100
+++ b/lemon/howard_mmc.h	Fri Feb 22 14:12:48 2013 +0100
@@ -149,6 +149,23 @@
     /// The \ref HowardMmcDefaultTraits "traits class" of the algorithm
     typedef TR Traits;
 
+    /// \brief Constants for the causes of search termination.
+    ///
+    /// Enum type containing constants for the different causes of search
+    /// termination. The \ref findCycleMean() function returns one of
+    /// these values.
+    enum TerminationCause {
+      
+      /// No directed cycle can be found in the digraph.
+      NO_CYCLE = 0,
+    
+      /// Optimal solution (minimum cycle mean) is found.
+      OPTIMAL = 1,
+
+      /// The iteration count limit is reached.
+      ITERATION_LIMIT
+    };
+
   private:
 
     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
@@ -324,25 +341,46 @@
       return findCycleMean() && findCycle();
     }
 
-    /// \brief Find the minimum cycle mean.
+    /// \brief Find the minimum cycle mean (or an upper bound).
     ///
     /// This function finds the minimum mean cost of the directed
-    /// cycles in the digraph.
+    /// cycles in the digraph (or an upper bound for it).
     ///
-    /// \return \c true if a directed cycle exists in the digraph.
-    bool findCycleMean() {
+    /// By default, the function finds the exact minimum cycle mean,
+    /// but an optional limit can also be specified for the number of
+    /// iterations performed during the search process.
+    /// The return value indicates if the optimal solution is found
+    /// or the iteration limit is reached. In the latter case, an
+    /// approximate solution is provided, which corresponds to a directed
+    /// cycle whose mean cost is relatively small, but not necessarily
+    /// minimal.
+    ///
+    /// \param limit  The maximum allowed number of iterations during
+    /// the search process. Its default value implies that the algorithm 
+    /// runs until it finds the exact optimal solution.
+    ///
+    /// \return The termination cause of the search process.
+    /// For more information, see \ref TerminationCause. 
+    TerminationCause findCycleMean(int limit = std::numeric_limits<int>::max()) {
       // Initialize and find strongly connected components
       init();
       findComponents();
 
       // Find the minimum cycle mean in the components
+      int iter_count = 0;
+      bool iter_limit_reached = false;
       for (int comp = 0; comp < _comp_num; ++comp) {
         // Find the minimum mean cycle in the current component
         if (!buildPolicyGraph(comp)) continue;
         while (true) {
+          if (++iter_count > limit) {
+            iter_limit_reached = true;
+            break;
+          }
           findPolicyCycle();
           if (!computeNodeDistances()) break;
         }
+
         // Update the best cycle (global minimum mean cycle)
         if ( _curr_found && (!_best_found ||
              _curr_cost * _best_size < _best_cost * _curr_size) ) {
@@ -351,8 +389,15 @@
           _best_size = _curr_size;
           _best_node = _curr_node;
         }
+        
+        if (iter_limit_reached) break;
       }
-      return _best_found;
+
+      if (iter_limit_reached) {
+        return ITERATION_LIMIT;
+      } else {
+        return _best_found ? OPTIMAL : NO_CYCLE;
+      }
     }
 
     /// \brief Find a minimum mean directed cycle.
diff -r 9a6c9d4ee77b -r bc726f4892c7 test/min_mean_cycle_test.cc
--- a/test/min_mean_cycle_test.cc	Wed Nov 28 12:08:47 2012 +0100
+++ b/test/min_mean_cycle_test.cc	Fri Feb 22 14:12:48 2013 +0100
@@ -110,7 +110,7 @@
                  const SmartDigraph::ArcMap<int>& cm,
                  int cost, int size) {
   MMC alg(gr, lm);
-  alg.findCycleMean();
+  check(alg.findCycleMean(), "Wrong result");
   check(alg.cycleMean() == static_cast<double>(cost) / size,
         "Wrong cycle mean");
   alg.findCycle();
@@ -210,6 +210,13 @@
     checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l2, c2,  5, 2);
     checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l3, c3,  0, 1);
     checkMmcAlg<HowardMmc<GR, IntArcMap> >(gr, l4, c4, -1, 1);
+    
+    // Howard with iteration limit
+    HowardMmc<GR, IntArcMap> mmc(gr, l1);
+    check((mmc.findCycleMean(2) == HowardMmc<GR, IntArcMap>::ITERATION_LIMIT),
+      "Wrong termination cause");
+    check((mmc.findCycleMean(4) == HowardMmc<GR, IntArcMap>::OPTIMAL),
+      "Wrong termination cause");
   }
 
   return 0;