1.1 --- a/lemon/cycle_canceling.h Mon Jul 16 16:21:40 2018 +0200
1.2 +++ b/lemon/cycle_canceling.h Wed Oct 17 19:14:07 2018 +0200
1.3 @@ -2,7 +2,7 @@
1.4 *
1.5 * This file is a part of LEMON, a generic C++ optimization library.
1.6 *
1.7 - * Copyright (C) 2003-2010
1.8 + * Copyright (C) 2003-2013
1.9 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 *
1.12 @@ -35,6 +35,7 @@
1.13 #include <lemon/circulation.h>
1.14 #include <lemon/bellman_ford.h>
1.15 #include <lemon/howard_mmc.h>
1.16 +#include <lemon/hartmann_orlin_mmc.h>
1.17
1.18 namespace lemon {
1.19
1.20 @@ -46,13 +47,14 @@
1.21 ///
1.22 /// \ref CycleCanceling implements three different cycle-canceling
1.23 /// algorithms for finding a \ref min_cost_flow "minimum cost flow"
1.24 - /// \ref amo93networkflows, \ref klein67primal,
1.25 - /// \ref goldberg89cyclecanceling.
1.26 - /// The most efficent one (both theoretically and practically)
1.27 - /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
1.28 - /// thus it is the default method.
1.29 - /// It is strongly polynomial, but in practice, it is typically much
1.30 - /// slower than the scaling algorithms and NetworkSimplex.
1.31 + /// \cite amo93networkflows, \cite klein67primal,
1.32 + /// \cite goldberg89cyclecanceling.
1.33 + /// The most efficent one is the \ref CANCEL_AND_TIGHTEN
1.34 + /// "Cancel-and-Tighten" algorithm, thus it is the default method.
1.35 + /// It runs in strongly polynomial time \f$O(n^2 m^2 \log n)\f$,
1.36 + /// but in practice, it is typically orders of magnitude slower than
1.37 + /// the scaling algorithms and \ref NetworkSimplex.
1.38 + /// (For more information, see \ref min_cost_flow_algs "the module page".)
1.39 ///
1.40 /// Most of the parameters of the problem (except for the digraph)
1.41 /// can be given using separate functions, and the algorithm can be
1.42 @@ -65,10 +67,11 @@
1.43 /// \tparam C The number type used for costs and potentials in the
1.44 /// algorithm. By default, it is the same as \c V.
1.45 ///
1.46 - /// \warning Both number types must be signed and all input data must
1.47 + /// \warning Both \c V and \c C must be signed number types.
1.48 + /// \warning All input data (capacities, supply values, and costs) must
1.49 /// be integer.
1.50 - /// \warning This algorithm does not support negative costs for such
1.51 - /// arcs that have infinite upper bound.
1.52 + /// \warning This algorithm does not support negative costs for
1.53 + /// arcs having infinite upper bound.
1.54 ///
1.55 /// \note For more information about the three available methods,
1.56 /// see \ref Method.
1.57 @@ -115,27 +118,28 @@
1.58 /// for the \ref run() function.
1.59 ///
1.60 /// \ref CycleCanceling provides three different cycle-canceling
1.61 - /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
1.62 - /// is used, which proved to be the most efficient and the most robust
1.63 - /// on various test inputs.
1.64 + /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel-and-Tighten"
1.65 + /// is used, which is by far the most efficient and the most robust.
1.66 /// However, the other methods can be selected using the \ref run()
1.67 /// function with the proper parameter.
1.68 enum Method {
1.69 /// A simple cycle-canceling method, which uses the
1.70 - /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
1.71 - /// number for detecting negative cycles in the residual network.
1.72 + /// \ref BellmanFord "Bellman-Ford" algorithm for detecting negative
1.73 + /// cycles in the residual network.
1.74 + /// The number of Bellman-Ford iterations is bounded by a successively
1.75 + /// increased limit.
1.76 SIMPLE_CYCLE_CANCELING,
1.77 /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
1.78 /// well-known strongly polynomial method
1.79 - /// \ref goldberg89cyclecanceling. It improves along a
1.80 + /// \cite goldberg89cyclecanceling. It improves along a
1.81 /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
1.82 - /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
1.83 + /// Its running time complexity is \f$O(n^2 m^3 \log n)\f$.
1.84 MINIMUM_MEAN_CYCLE_CANCELING,
1.85 - /// The "Cancel And Tighten" algorithm, which can be viewed as an
1.86 + /// The "Cancel-and-Tighten" algorithm, which can be viewed as an
1.87 /// improved version of the previous method
1.88 - /// \ref goldberg89cyclecanceling.
1.89 + /// \cite goldberg89cyclecanceling.
1.90 /// It is faster both in theory and in practice, its running time
1.91 - /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
1.92 + /// complexity is \f$O(n^2 m^2 \log n)\f$.
1.93 CANCEL_AND_TIGHTEN
1.94 };
1.95
1.96 @@ -191,7 +195,7 @@
1.97 int _root;
1.98
1.99 // Parameters of the problem
1.100 - bool _have_lower;
1.101 + bool _has_lower;
1.102 Value _sum_supply;
1.103
1.104 // Data structures for storing the digraph
1.105 @@ -274,10 +278,9 @@
1.106 /// \return <tt>(*this)</tt>
1.107 template <typename LowerMap>
1.108 CycleCanceling& lowerMap(const LowerMap& map) {
1.109 - _have_lower = true;
1.110 + _has_lower = true;
1.111 for (ArcIt a(_graph); a != INVALID; ++a) {
1.112 _lower[_arc_idf[a]] = map[a];
1.113 - _lower[_arc_idb[a]] = map[a];
1.114 }
1.115 return *this;
1.116 }
1.117 @@ -349,7 +352,7 @@
1.118 /// calling \ref run(), the supply of each node will be set to zero.
1.119 ///
1.120 /// Using this function has the same effect as using \ref supplyMap()
1.121 - /// with such a map in which \c k is assigned to \c s, \c -k is
1.122 + /// with a map in which \c k is assigned to \c s, \c -k is
1.123 /// assigned to \c t and all other nodes have zero supply value.
1.124 ///
1.125 /// \param s The source node.
1.126 @@ -467,7 +470,7 @@
1.127 _cost[j] = 0;
1.128 _cost[_reverse[j]] = 0;
1.129 }
1.130 - _have_lower = false;
1.131 + _has_lower = false;
1.132 return *this;
1.133 }
1.134
1.135 @@ -572,7 +575,7 @@
1.136 /// \brief Return the total cost of the found flow.
1.137 ///
1.138 /// This function returns the total cost of the found flow.
1.139 - /// Its complexity is O(e).
1.140 + /// Its complexity is O(m).
1.141 ///
1.142 /// \note The return type of the function can be specified as a
1.143 /// template parameter. For example,
1.144 @@ -610,7 +613,8 @@
1.145 return _res_cap[_arc_idb[a]];
1.146 }
1.147
1.148 - /// \brief Return the flow map (the primal solution).
1.149 + /// \brief Copy the flow values (the primal solution) into the
1.150 + /// given map.
1.151 ///
1.152 /// This function copies the flow value on each arc into the given
1.153 /// map. The \c Value type of the algorithm must be convertible to
1.154 @@ -634,7 +638,8 @@
1.155 return static_cast<Cost>(_pi[_node_id[n]]);
1.156 }
1.157
1.158 - /// \brief Return the potential map (the dual solution).
1.159 + /// \brief Copy the potential values (the dual solution) into the
1.160 + /// given map.
1.161 ///
1.162 /// This function copies the potential (dual value) of each node
1.163 /// into the given map.
1.164 @@ -664,6 +669,10 @@
1.165 }
1.166 if (_sum_supply > 0) return INFEASIBLE;
1.167
1.168 + // Check lower and upper bounds
1.169 + LEMON_DEBUG(checkBoundMaps(),
1.170 + "Upper bounds must be greater or equal to the lower bounds");
1.171 +
1.172
1.173 // Initialize vectors
1.174 for (int i = 0; i != _res_node_num; ++i) {
1.175 @@ -674,7 +683,7 @@
1.176 // Remove infinite upper bounds and check negative arcs
1.177 const Value MAX = std::numeric_limits<Value>::max();
1.178 int last_out;
1.179 - if (_have_lower) {
1.180 + if (_has_lower) {
1.181 for (int i = 0; i != _root; ++i) {
1.182 last_out = _first_out[i+1];
1.183 for (int j = _first_out[i]; j != last_out; ++j) {
1.184 @@ -717,7 +726,7 @@
1.185 for (NodeIt n(_graph); n != INVALID; ++n) {
1.186 sup[n] = _supply[_node_id[n]];
1.187 }
1.188 - if (_have_lower) {
1.189 + if (_has_lower) {
1.190 for (ArcIt a(_graph); a != INVALID; ++a) {
1.191 int j = _arc_idf[a];
1.192 Value c = _lower[j];
1.193 @@ -774,6 +783,15 @@
1.194 return OPTIMAL;
1.195 }
1.196
1.197 + // Check if the upper bound is greater than or equal to the lower bound
1.198 + // on each forward arc.
1.199 + bool checkBoundMaps() {
1.200 + for (int j = 0; j != _res_arc_num; ++j) {
1.201 + if (_forward[j] && _upper[j] < _lower[j]) return false;
1.202 + }
1.203 + return true;
1.204 + }
1.205 +
1.206 // Build a StaticDigraph structure containing the current
1.207 // residual network
1.208 void buildResidualNetwork() {
1.209 @@ -816,10 +834,10 @@
1.210 }
1.211
1.212 // Handle non-zero lower bounds
1.213 - if (_have_lower) {
1.214 + if (_has_lower) {
1.215 int limit = _first_out[_root];
1.216 for (int j = 0; j != limit; ++j) {
1.217 - if (!_forward[j]) _res_cap[j] += _lower[j];
1.218 + if (_forward[j]) _res_cap[_reverse[j]] += _lower[j];
1.219 }
1.220 }
1.221 }
1.222 @@ -922,18 +940,41 @@
1.223
1.224 // Execute the "Minimum Mean Cycle Canceling" method
1.225 void startMinMeanCycleCanceling() {
1.226 - typedef SimplePath<StaticDigraph> SPath;
1.227 + typedef Path<StaticDigraph> SPath;
1.228 typedef typename SPath::ArcIt SPathArcIt;
1.229 typedef typename HowardMmc<StaticDigraph, CostArcMap>
1.230 - ::template SetPath<SPath>::Create MMC;
1.231 + ::template SetPath<SPath>::Create HwMmc;
1.232 + typedef typename HartmannOrlinMmc<StaticDigraph, CostArcMap>
1.233 + ::template SetPath<SPath>::Create HoMmc;
1.234 +
1.235 + const double HW_ITER_LIMIT_FACTOR = 1.0;
1.236 + const int HW_ITER_LIMIT_MIN_VALUE = 5;
1.237 +
1.238 + const int hw_iter_limit =
1.239 + std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num),
1.240 + HW_ITER_LIMIT_MIN_VALUE);
1.241
1.242 SPath cycle;
1.243 - MMC mmc(_sgr, _cost_map);
1.244 - mmc.cycle(cycle);
1.245 + HwMmc hw_mmc(_sgr, _cost_map);
1.246 + hw_mmc.cycle(cycle);
1.247 buildResidualNetwork();
1.248 - while (mmc.findCycleMean() && mmc.cycleCost() < 0) {
1.249 - // Find the cycle
1.250 - mmc.findCycle();
1.251 + while (true) {
1.252 +
1.253 + typename HwMmc::TerminationCause hw_tc =
1.254 + hw_mmc.findCycleMean(hw_iter_limit);
1.255 + if (hw_tc == HwMmc::ITERATION_LIMIT) {
1.256 + // Howard's algorithm reached the iteration limit, start a
1.257 + // strongly polynomial algorithm instead
1.258 + HoMmc ho_mmc(_sgr, _cost_map);
1.259 + ho_mmc.cycle(cycle);
1.260 + // Find a minimum mean cycle (Hartmann-Orlin algorithm)
1.261 + if (!(ho_mmc.findCycleMean() && ho_mmc.cycleCost() < 0)) break;
1.262 + ho_mmc.findCycle();
1.263 + } else {
1.264 + // Find a minimum mean cycle (Howard algorithm)
1.265 + if (!(hw_tc == HwMmc::OPTIMAL && hw_mmc.cycleCost() < 0)) break;
1.266 + hw_mmc.findCycle();
1.267 + }
1.268
1.269 // Compute delta value
1.270 Value delta = INF;
1.271 @@ -954,11 +995,17 @@
1.272 }
1.273 }
1.274
1.275 - // Execute the "Cancel And Tighten" method
1.276 + // Execute the "Cancel-and-Tighten" method
1.277 void startCancelAndTighten() {
1.278 // Constants for the min mean cycle computations
1.279 const double LIMIT_FACTOR = 1.0;
1.280 const int MIN_LIMIT = 5;
1.281 + const double HW_ITER_LIMIT_FACTOR = 1.0;
1.282 + const int HW_ITER_LIMIT_MIN_VALUE = 5;
1.283 +
1.284 + const int hw_iter_limit =
1.285 + std::max(static_cast<int>(HW_ITER_LIMIT_FACTOR * _node_num),
1.286 + HW_ITER_LIMIT_MIN_VALUE);
1.287
1.288 // Contruct auxiliary data vectors
1.289 DoubleVector pi(_res_node_num, 0.0);
1.290 @@ -1132,17 +1179,30 @@
1.291 }
1.292 }
1.293 } else {
1.294 - typedef HowardMmc<StaticDigraph, CostArcMap> MMC;
1.295 + typedef HowardMmc<StaticDigraph, CostArcMap> HwMmc;
1.296 + typedef HartmannOrlinMmc<StaticDigraph, CostArcMap> HoMmc;
1.297 typedef typename BellmanFord<StaticDigraph, CostArcMap>
1.298 ::template SetDistMap<CostNodeMap>::Create BF;
1.299
1.300 // Set epsilon to the minimum cycle mean
1.301 + Cost cycle_cost = 0;
1.302 + int cycle_size = 1;
1.303 buildResidualNetwork();
1.304 - MMC mmc(_sgr, _cost_map);
1.305 - mmc.findCycleMean();
1.306 - epsilon = -mmc.cycleMean();
1.307 - Cost cycle_cost = mmc.cycleCost();
1.308 - int cycle_size = mmc.cycleSize();
1.309 + HwMmc hw_mmc(_sgr, _cost_map);
1.310 + if (hw_mmc.findCycleMean(hw_iter_limit) == HwMmc::ITERATION_LIMIT) {
1.311 + // Howard's algorithm reached the iteration limit, start a
1.312 + // strongly polynomial algorithm instead
1.313 + HoMmc ho_mmc(_sgr, _cost_map);
1.314 + ho_mmc.findCycleMean();
1.315 + epsilon = -ho_mmc.cycleMean();
1.316 + cycle_cost = ho_mmc.cycleCost();
1.317 + cycle_size = ho_mmc.cycleSize();
1.318 + } else {
1.319 + // Set epsilon
1.320 + epsilon = -hw_mmc.cycleMean();
1.321 + cycle_cost = hw_mmc.cycleCost();
1.322 + cycle_size = hw_mmc.cycleSize();
1.323 + }
1.324
1.325 // Compute feasible potentials for the current epsilon
1.326 for (int i = 0; i != int(_cost_vec.size()); ++i) {