1.1 --- a/lemon/capacity_scaling.h Mon Jul 16 16:21:40 2018 +0200
1.2 +++ b/lemon/capacity_scaling.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 @@ -27,6 +27,7 @@
1.13 #include <vector>
1.14 #include <limits>
1.15 #include <lemon/core.h>
1.16 +#include <lemon/maps.h>
1.17 #include <lemon/bin_heap.h>
1.18
1.19 namespace lemon {
1.20 @@ -66,9 +67,16 @@
1.21 ///
1.22 /// \ref CapacityScaling implements the capacity scaling version
1.23 /// of the successive shortest path algorithm for finding a
1.24 - /// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,
1.25 - /// \ref edmondskarp72theoretical. It is an efficient dual
1.26 - /// solution method.
1.27 + /// \ref min_cost_flow "minimum cost flow" \cite amo93networkflows,
1.28 + /// \cite edmondskarp72theoretical. It is an efficient dual
1.29 + /// solution method, which runs in polynomial time
1.30 + /// \f$O(m\log U (n+m)\log n)\f$, where <i>U</i> denotes the maximum
1.31 + /// of node supply and arc capacity values.
1.32 + ///
1.33 + /// This algorithm is typically slower than \ref CostScaling and
1.34 + /// \ref NetworkSimplex, but in special cases, it can be more
1.35 + /// efficient than them.
1.36 + /// (For more information, see \ref min_cost_flow_algs "the module page".)
1.37 ///
1.38 /// Most of the parameters of the problem (except for the digraph)
1.39 /// can be given using separate functions, and the algorithm can be
1.40 @@ -86,10 +94,11 @@
1.41 /// In most cases, this parameter should not be set directly,
1.42 /// consider to use the named template parameters instead.
1.43 ///
1.44 - /// \warning Both number types must be signed and all input data must
1.45 - /// be integer.
1.46 - /// \warning This algorithm does not support negative costs for such
1.47 - /// arcs that have infinite upper bound.
1.48 + /// \warning Both \c V and \c C must be signed number types.
1.49 + /// \warning Capacity bounds and supply values must be integer, but
1.50 + /// arc costs can be arbitrary real numbers.
1.51 + /// \warning This algorithm does not support negative costs for
1.52 + /// arcs having infinite upper bound.
1.53 #ifdef DOXYGEN
1.54 template <typename GR, typename V, typename C, typename TR>
1.55 #else
1.56 @@ -110,7 +119,8 @@
1.57 /// The type of the heap used for internal Dijkstra computations
1.58 typedef typename TR::Heap Heap;
1.59
1.60 - /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
1.61 + /// \brief The \ref lemon::CapacityScalingDefaultTraits "traits class"
1.62 + /// of the algorithm
1.63 typedef TR Traits;
1.64
1.65 public:
1.66 @@ -154,7 +164,7 @@
1.67 int _root;
1.68
1.69 // Parameters of the problem
1.70 - bool _have_lower;
1.71 + bool _has_lower;
1.72 Value _sum_supply;
1.73
1.74 // Data structures for storing the digraph
1.75 @@ -347,10 +357,9 @@
1.76 /// \return <tt>(*this)</tt>
1.77 template <typename LowerMap>
1.78 CapacityScaling& lowerMap(const LowerMap& map) {
1.79 - _have_lower = true;
1.80 + _has_lower = true;
1.81 for (ArcIt a(_graph); a != INVALID; ++a) {
1.82 _lower[_arc_idf[a]] = map[a];
1.83 - _lower[_arc_idb[a]] = map[a];
1.84 }
1.85 return *this;
1.86 }
1.87 @@ -422,7 +431,7 @@
1.88 /// calling \ref run(), the supply of each node will be set to zero.
1.89 ///
1.90 /// Using this function has the same effect as using \ref supplyMap()
1.91 - /// with such a map in which \c k is assigned to \c s, \c -k is
1.92 + /// with a map in which \c k is assigned to \c s, \c -k is
1.93 /// assigned to \c t and all other nodes have zero supply value.
1.94 ///
1.95 /// \param s The source node.
1.96 @@ -534,7 +543,7 @@
1.97 _upper[j] = INF;
1.98 _cost[j] = _forward[j] ? 1 : -1;
1.99 }
1.100 - _have_lower = false;
1.101 + _has_lower = false;
1.102 return *this;
1.103 }
1.104
1.105 @@ -637,7 +646,7 @@
1.106 /// \brief Return the total cost of the found flow.
1.107 ///
1.108 /// This function returns the total cost of the found flow.
1.109 - /// Its complexity is O(e).
1.110 + /// Its complexity is O(m).
1.111 ///
1.112 /// \note The return type of the function can be specified as a
1.113 /// template parameter. For example,
1.114 @@ -675,7 +684,8 @@
1.115 return _res_cap[_arc_idb[a]];
1.116 }
1.117
1.118 - /// \brief Return the flow map (the primal solution).
1.119 + /// \brief Copy the flow values (the primal solution) into the
1.120 + /// given map.
1.121 ///
1.122 /// This function copies the flow value on each arc into the given
1.123 /// map. The \c Value type of the algorithm must be convertible to
1.124 @@ -699,7 +709,8 @@
1.125 return _pi[_node_id[n]];
1.126 }
1.127
1.128 - /// \brief Return the potential map (the dual solution).
1.129 + /// \brief Copy the potential values (the dual solution) into the
1.130 + /// given map.
1.131 ///
1.132 /// This function copies the potential (dual value) of each node
1.133 /// into the given map.
1.134 @@ -729,6 +740,11 @@
1.135 }
1.136 if (_sum_supply > 0) return INFEASIBLE;
1.137
1.138 + // Check lower and upper bounds
1.139 + LEMON_DEBUG(checkBoundMaps(),
1.140 + "Upper bounds must be greater or equal to the lower bounds");
1.141 +
1.142 +
1.143 // Initialize vectors
1.144 for (int i = 0; i != _root; ++i) {
1.145 _pi[i] = 0;
1.146 @@ -738,7 +754,7 @@
1.147 // Remove non-zero lower bounds
1.148 const Value MAX = std::numeric_limits<Value>::max();
1.149 int last_out;
1.150 - if (_have_lower) {
1.151 + if (_has_lower) {
1.152 for (int i = 0; i != _root; ++i) {
1.153 last_out = _first_out[i+1];
1.154 for (int j = _first_out[i]; j != last_out; ++j) {
1.155 @@ -823,6 +839,15 @@
1.156 return OPTIMAL;
1.157 }
1.158
1.159 + // Check if the upper bound is greater than or equal to the lower bound
1.160 + // on each forward arc.
1.161 + bool checkBoundMaps() {
1.162 + for (int j = 0; j != _res_arc_num; ++j) {
1.163 + if (_forward[j] && _upper[j] < _lower[j]) return false;
1.164 + }
1.165 + return true;
1.166 + }
1.167 +
1.168 ProblemType start() {
1.169 // Execute the algorithm
1.170 ProblemType pt;
1.171 @@ -832,10 +857,10 @@
1.172 pt = startWithoutScaling();
1.173
1.174 // Handle non-zero lower bounds
1.175 - if (_have_lower) {
1.176 + if (_has_lower) {
1.177 int limit = _first_out[_root];
1.178 for (int j = 0; j != limit; ++j) {
1.179 - if (!_forward[j]) _res_cap[j] += _lower[j];
1.180 + if (_forward[j]) _res_cap[_reverse[j]] += _lower[j];
1.181 }
1.182 }
1.183