1.1 --- a/lemon/network_simplex.h Fri Apr 03 18:59:15 2009 +0200
1.2 +++ b/lemon/network_simplex.h Fri Apr 17 18:04:36 2009 +0200
1.3 @@ -30,6 +30,9 @@
1.4
1.5 #include <lemon/core.h>
1.6 #include <lemon/math.h>
1.7 +#include <lemon/maps.h>
1.8 +#include <lemon/circulation.h>
1.9 +#include <lemon/adaptors.h>
1.10
1.11 namespace lemon {
1.12
1.13 @@ -47,6 +50,8 @@
1.14 ///
1.15 /// In general this class is the fastest implementation available
1.16 /// in LEMON for the minimum cost flow problem.
1.17 + /// Moreover it supports both direction of the supply/demand inequality
1.18 + /// constraints. For more information see \ref ProblemType.
1.19 ///
1.20 /// \tparam GR The digraph type the algorithm runs on.
1.21 /// \tparam F The value type used for flow amounts, capacity bounds
1.22 @@ -58,7 +63,8 @@
1.23 /// be integer.
1.24 ///
1.25 /// \note %NetworkSimplex provides five different pivot rule
1.26 - /// implementations. For more information see \ref PivotRule.
1.27 + /// implementations, from which the most efficient one is used
1.28 + /// by default. For more information see \ref PivotRule.
1.29 template <typename GR, typename F = int, typename C = F>
1.30 class NetworkSimplex
1.31 {
1.32 @@ -68,10 +74,17 @@
1.33 typedef F Flow;
1.34 /// The cost type of the algorithm
1.35 typedef C Cost;
1.36 +#ifdef DOXYGEN
1.37 + /// The type of the flow map
1.38 + typedef GR::ArcMap<Flow> FlowMap;
1.39 + /// The type of the potential map
1.40 + typedef GR::NodeMap<Cost> PotentialMap;
1.41 +#else
1.42 /// The type of the flow map
1.43 typedef typename GR::template ArcMap<Flow> FlowMap;
1.44 /// The type of the potential map
1.45 typedef typename GR::template NodeMap<Cost> PotentialMap;
1.46 +#endif
1.47
1.48 public:
1.49
1.50 @@ -117,6 +130,58 @@
1.51 /// candidate list and extends this list in every iteration.
1.52 ALTERING_LIST
1.53 };
1.54 +
1.55 + /// \brief Enum type for selecting the problem type.
1.56 + ///
1.57 + /// Enum type for selecting the problem type, i.e. the direction of
1.58 + /// the inequalities in the supply/demand constraints of the
1.59 + /// \ref min_cost_flow "minimum cost flow problem".
1.60 + ///
1.61 + /// The default problem type is \c GEQ, since this form is supported
1.62 + /// by other minimum cost flow algorithms and the \ref Circulation
1.63 + /// algorithm as well.
1.64 + /// The \c LEQ problem type can be selected using the \ref problemType()
1.65 + /// function.
1.66 + ///
1.67 + /// Note that the equality form is a special case of both problem type.
1.68 + enum ProblemType {
1.69 +
1.70 + /// This option means that there are "<em>greater or equal</em>"
1.71 + /// constraints in the defintion, i.e. the exact formulation of the
1.72 + /// problem is the following.
1.73 + /**
1.74 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.75 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
1.76 + sup(u) \quad \forall u\in V \f]
1.77 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.78 + */
1.79 + /// It means that the total demand must be greater or equal to the
1.80 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.81 + /// negative) and all the supplies have to be carried out from
1.82 + /// the supply nodes, but there could be demands that are not
1.83 + /// satisfied.
1.84 + GEQ,
1.85 + /// It is just an alias for the \c GEQ option.
1.86 + CARRY_SUPPLIES = GEQ,
1.87 +
1.88 + /// This option means that there are "<em>less or equal</em>"
1.89 + /// constraints in the defintion, i.e. the exact formulation of the
1.90 + /// problem is the following.
1.91 + /**
1.92 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.93 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
1.94 + sup(u) \quad \forall u\in V \f]
1.95 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.96 + */
1.97 + /// It means that the total demand must be less or equal to the
1.98 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
1.99 + /// positive) and all the demands have to be satisfied, but there
1.100 + /// could be supplies that are not carried out from the supply
1.101 + /// nodes.
1.102 + LEQ,
1.103 + /// It is just an alias for the \c LEQ option.
1.104 + SATISFY_DEMANDS = LEQ
1.105 + };
1.106
1.107 private:
1.108
1.109 @@ -155,6 +220,7 @@
1.110 bool _pstsup;
1.111 Node _psource, _ptarget;
1.112 Flow _pstflow;
1.113 + ProblemType _ptype;
1.114
1.115 // Result maps
1.116 FlowMap *_flow_map;
1.117 @@ -586,13 +652,13 @@
1.118
1.119 /// \brief Constructor.
1.120 ///
1.121 - /// Constructor.
1.122 + /// The constructor of the class.
1.123 ///
1.124 /// \param graph The digraph the algorithm runs on.
1.125 NetworkSimplex(const GR& graph) :
1.126 _graph(graph),
1.127 _plower(NULL), _pupper(NULL), _pcost(NULL),
1.128 - _psupply(NULL), _pstsup(false),
1.129 + _psupply(NULL), _pstsup(false), _ptype(GEQ),
1.130 _flow_map(NULL), _potential_map(NULL),
1.131 _local_flow(false), _local_potential(false),
1.132 _node_id(graph)
1.133 @@ -611,6 +677,12 @@
1.134 if (_local_potential) delete _potential_map;
1.135 }
1.136
1.137 + /// \name Parameters
1.138 + /// The parameters of the algorithm can be specified using these
1.139 + /// functions.
1.140 +
1.141 + /// @{
1.142 +
1.143 /// \brief Set the lower bounds on the arcs.
1.144 ///
1.145 /// This function sets the lower bounds on the arcs.
1.146 @@ -760,6 +832,20 @@
1.147 _pstflow = k;
1.148 return *this;
1.149 }
1.150 +
1.151 + /// \brief Set the problem type.
1.152 + ///
1.153 + /// This function sets the problem type for the algorithm.
1.154 + /// If it is not used before calling \ref run(), the \ref GEQ problem
1.155 + /// type will be used.
1.156 + ///
1.157 + /// For more information see \ref ProblemType.
1.158 + ///
1.159 + /// \return <tt>(*this)</tt>
1.160 + NetworkSimplex& problemType(ProblemType problem_type) {
1.161 + _ptype = problem_type;
1.162 + return *this;
1.163 + }
1.164
1.165 /// \brief Set the flow map.
1.166 ///
1.167 @@ -795,6 +881,8 @@
1.168 _potential_map = ↦
1.169 return *this;
1.170 }
1.171 +
1.172 + /// @}
1.173
1.174 /// \name Execution Control
1.175 /// The algorithm can be executed using \ref run().
1.176 @@ -804,10 +892,11 @@
1.177 /// \brief Run the algorithm.
1.178 ///
1.179 /// This function runs the algorithm.
1.180 - /// The paramters can be specified using \ref lowerMap(),
1.181 + /// The paramters can be specified using functions \ref lowerMap(),
1.182 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
1.183 - /// \ref costMap(), \ref supplyMap() and \ref stSupply()
1.184 - /// functions. For example,
1.185 + /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
1.186 + /// \ref problemType(), \ref flowMap() and \ref potentialMap().
1.187 + /// For example,
1.188 /// \code
1.189 /// NetworkSimplex<ListDigraph> ns(graph);
1.190 /// ns.boundMaps(lower, upper).costMap(cost)
1.191 @@ -830,9 +919,10 @@
1.192 /// \brief Reset all the parameters that have been given before.
1.193 ///
1.194 /// This function resets all the paramaters that have been given
1.195 - /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
1.196 - /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
1.197 - /// \ref stSupply() functions before.
1.198 + /// before using functions \ref lowerMap(), \ref upperMap(),
1.199 + /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
1.200 + /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
1.201 + /// \ref flowMap() and \ref potentialMap().
1.202 ///
1.203 /// It is useful for multiple run() calls. If this function is not
1.204 /// used, all the parameters given before are kept for the next
1.205 @@ -869,6 +959,14 @@
1.206 _pcost = NULL;
1.207 _psupply = NULL;
1.208 _pstsup = false;
1.209 + _ptype = GEQ;
1.210 + if (_local_flow) delete _flow_map;
1.211 + if (_local_potential) delete _potential_map;
1.212 + _flow_map = NULL;
1.213 + _potential_map = NULL;
1.214 + _local_flow = false;
1.215 + _local_potential = false;
1.216 +
1.217 return *this;
1.218 }
1.219
1.220 @@ -1000,20 +1098,21 @@
1.221
1.222 // Initialize node related data
1.223 bool valid_supply = true;
1.224 + Flow sum_supply = 0;
1.225 if (!_pstsup && !_psupply) {
1.226 _pstsup = true;
1.227 _psource = _ptarget = NodeIt(_graph);
1.228 _pstflow = 0;
1.229 }
1.230 if (_psupply) {
1.231 - Flow sum = 0;
1.232 int i = 0;
1.233 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.234 _node_id[n] = i;
1.235 _supply[i] = (*_psupply)[n];
1.236 - sum += _supply[i];
1.237 + sum_supply += _supply[i];
1.238 }
1.239 - valid_supply = (sum == 0);
1.240 + valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
1.241 + (_ptype == LEQ && sum_supply >= 0);
1.242 } else {
1.243 int i = 0;
1.244 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1.245 @@ -1021,10 +1120,95 @@
1.246 _supply[i] = 0;
1.247 }
1.248 _supply[_node_id[_psource]] = _pstflow;
1.249 - _supply[_node_id[_ptarget]] = -_pstflow;
1.250 + _supply[_node_id[_ptarget]] = -_pstflow;
1.251 }
1.252 if (!valid_supply) return false;
1.253
1.254 + // Infinite capacity value
1.255 + Flow inf_cap =
1.256 + std::numeric_limits<Flow>::has_infinity ?
1.257 + std::numeric_limits<Flow>::infinity() :
1.258 + std::numeric_limits<Flow>::max();
1.259 +
1.260 + // Initialize artifical cost
1.261 + Cost art_cost;
1.262 + if (std::numeric_limits<Cost>::is_exact) {
1.263 + art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1.264 + } else {
1.265 + art_cost = std::numeric_limits<Cost>::min();
1.266 + for (int i = 0; i != _arc_num; ++i) {
1.267 + if (_cost[i] > art_cost) art_cost = _cost[i];
1.268 + }
1.269 + art_cost = (art_cost + 1) * _node_num;
1.270 + }
1.271 +
1.272 + // Run Circulation to check if a feasible solution exists
1.273 + typedef ConstMap<Arc, Flow> ConstArcMap;
1.274 + FlowNodeMap *csup = NULL;
1.275 + bool local_csup = false;
1.276 + if (_psupply) {
1.277 + csup = _psupply;
1.278 + } else {
1.279 + csup = new FlowNodeMap(_graph, 0);
1.280 + (*csup)[_psource] = _pstflow;
1.281 + (*csup)[_ptarget] = -_pstflow;
1.282 + local_csup = true;
1.283 + }
1.284 + bool circ_result = false;
1.285 + if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
1.286 + // GEQ problem type
1.287 + if (_plower) {
1.288 + if (_pupper) {
1.289 + Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
1.290 + circ(_graph, *_plower, *_pupper, *csup);
1.291 + circ_result = circ.run();
1.292 + } else {
1.293 + Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
1.294 + circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
1.295 + circ_result = circ.run();
1.296 + }
1.297 + } else {
1.298 + if (_pupper) {
1.299 + Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
1.300 + circ(_graph, ConstArcMap(0), *_pupper, *csup);
1.301 + circ_result = circ.run();
1.302 + } else {
1.303 + Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
1.304 + circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
1.305 + circ_result = circ.run();
1.306 + }
1.307 + }
1.308 + } else {
1.309 + // LEQ problem type
1.310 + typedef ReverseDigraph<const GR> RevGraph;
1.311 + typedef NegMap<FlowNodeMap> NegNodeMap;
1.312 + RevGraph rgraph(_graph);
1.313 + NegNodeMap neg_csup(*csup);
1.314 + if (_plower) {
1.315 + if (_pupper) {
1.316 + Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
1.317 + circ(rgraph, *_plower, *_pupper, neg_csup);
1.318 + circ_result = circ.run();
1.319 + } else {
1.320 + Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
1.321 + circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
1.322 + circ_result = circ.run();
1.323 + }
1.324 + } else {
1.325 + if (_pupper) {
1.326 + Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
1.327 + circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
1.328 + circ_result = circ.run();
1.329 + } else {
1.330 + Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
1.331 + circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
1.332 + circ_result = circ.run();
1.333 + }
1.334 + }
1.335 + }
1.336 + if (local_csup) delete csup;
1.337 + if (!circ_result) return false;
1.338 +
1.339 // Set data for the artificial root node
1.340 _root = _node_num;
1.341 _parent[_root] = -1;
1.342 @@ -1033,8 +1217,12 @@
1.343 _rev_thread[0] = _root;
1.344 _succ_num[_root] = all_node_num;
1.345 _last_succ[_root] = _root - 1;
1.346 - _supply[_root] = 0;
1.347 - _pi[_root] = 0;
1.348 + _supply[_root] = -sum_supply;
1.349 + if (sum_supply < 0) {
1.350 + _pi[_root] = -art_cost;
1.351 + } else {
1.352 + _pi[_root] = art_cost;
1.353 + }
1.354
1.355 // Store the arcs in a mixed order
1.356 int k = std::max(int(sqrt(_arc_num)), 10);
1.357 @@ -1045,10 +1233,6 @@
1.358 }
1.359
1.360 // Initialize arc maps
1.361 - Flow inf_cap =
1.362 - std::numeric_limits<Flow>::has_infinity ?
1.363 - std::numeric_limits<Flow>::infinity() :
1.364 - std::numeric_limits<Flow>::max();
1.365 if (_pupper && _pcost) {
1.366 for (int i = 0; i != _arc_num; ++i) {
1.367 Arc e = _arc_ref[i];
1.368 @@ -1083,18 +1267,6 @@
1.369 }
1.370 }
1.371
1.372 - // Initialize artifical cost
1.373 - Cost art_cost;
1.374 - if (std::numeric_limits<Cost>::is_exact) {
1.375 - art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1.376 - } else {
1.377 - art_cost = std::numeric_limits<Cost>::min();
1.378 - for (int i = 0; i != _arc_num; ++i) {
1.379 - if (_cost[i] > art_cost) art_cost = _cost[i];
1.380 - }
1.381 - art_cost = (art_cost + 1) * _node_num;
1.382 - }
1.383 -
1.384 // Remove non-zero lower bounds
1.385 if (_plower) {
1.386 for (int i = 0; i != _arc_num; ++i) {
1.387 @@ -1118,14 +1290,14 @@
1.388 _cost[e] = art_cost;
1.389 _cap[e] = inf_cap;
1.390 _state[e] = STATE_TREE;
1.391 - if (_supply[u] >= 0) {
1.392 + if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
1.393 _flow[e] = _supply[u];
1.394 _forward[u] = true;
1.395 - _pi[u] = -art_cost;
1.396 + _pi[u] = -art_cost + _pi[_root];
1.397 } else {
1.398 _flow[e] = -_supply[u];
1.399 _forward[u] = false;
1.400 - _pi[u] = art_cost;
1.401 + _pi[u] = art_cost + _pi[_root];
1.402 }
1.403 }
1.404
1.405 @@ -1382,11 +1554,6 @@
1.406 }
1.407 }
1.408
1.409 - // Check if the flow amount equals zero on all the artificial arcs
1.410 - for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
1.411 - if (_flow[e] > 0) return false;
1.412 - }
1.413 -
1.414 // Copy flow values to _flow_map
1.415 if (_plower) {
1.416 for (int i = 0; i != _arc_num; ++i) {