1.1 --- a/doc/groups.dox Fri Apr 03 18:59:15 2009 +0200
1.2 +++ b/doc/groups.dox Fri Apr 17 18:04:36 2009 +0200
1.3 @@ -317,15 +317,15 @@
1.4
1.5 The \e maximum \e flow \e problem is to find a flow of maximum value between
1.6 a single source and a single target. Formally, there is a \f$G=(V,A)\f$
1.7 -digraph, a \f$cap:A\rightarrow\mathbf{R}^+_0\f$ capacity function and
1.8 +digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
1.9 \f$s, t \in V\f$ source and target nodes.
1.10 -A maximum flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of the
1.11 +A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
1.12 following optimization problem.
1.13
1.14 -\f[ \max\sum_{a\in\delta_{out}(s)}f(a) - \sum_{a\in\delta_{in}(s)}f(a) \f]
1.15 -\f[ \sum_{a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)
1.16 - \qquad \forall v\in V\setminus\{s,t\} \f]
1.17 -\f[ 0 \leq f(a) \leq cap(a) \qquad \forall a\in A \f]
1.18 +\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
1.19 +\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
1.20 + \quad \forall u\in V\setminus\{s,t\} \f]
1.21 +\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
1.22
1.23 LEMON contains several algorithms for solving maximum flow problems:
1.24 - \ref EdmondsKarp Edmonds-Karp algorithm.
1.25 @@ -345,35 +345,103 @@
1.26
1.27 \brief Algorithms for finding minimum cost flows and circulations.
1.28
1.29 -This group describes the algorithms for finding minimum cost flows and
1.30 +This group contains the algorithms for finding minimum cost flows and
1.31 circulations.
1.32
1.33 The \e minimum \e cost \e flow \e problem is to find a feasible flow of
1.34 minimum total cost from a set of supply nodes to a set of demand nodes
1.35 -in a network with capacity constraints and arc costs.
1.36 +in a network with capacity constraints (lower and upper bounds)
1.37 +and arc costs.
1.38 Formally, let \f$G=(V,A)\f$ be a digraph,
1.39 \f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
1.40 -upper bounds for the flow values on the arcs,
1.41 +upper bounds for the flow values on the arcs, for which
1.42 +\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
1.43 \f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
1.44 -on the arcs, and
1.45 -\f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values
1.46 -of the nodes.
1.47 -A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of
1.48 -the following optimization problem.
1.49 +on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
1.50 +signed supply values of the nodes.
1.51 +If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
1.52 +supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
1.53 +\f$-sup(u)\f$ demand.
1.54 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
1.55 +of the following optimization problem.
1.56
1.57 -\f[ \min\sum_{a\in A} f(a) cost(a) \f]
1.58 -\f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a) =
1.59 - supply(v) \qquad \forall v\in V \f]
1.60 -\f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]
1.61 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
1.62 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
1.63 + sup(u) \quad \forall u\in V \f]
1.64 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
1.65
1.66 -LEMON contains several algorithms for solving minimum cost flow problems:
1.67 - - \ref CycleCanceling Cycle-canceling algorithms.
1.68 - - \ref CapacityScaling Successive shortest path algorithm with optional
1.69 +The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
1.70 +zero or negative in order to have a feasible solution (since the sum
1.71 +of the expressions on the left-hand side of the inequalities is zero).
1.72 +It means that the total demand must be greater or equal to the total
1.73 +supply and all the supplies have to be carried out from the supply nodes,
1.74 +but there could be demands that are not satisfied.
1.75 +If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
1.76 +constraints have to be satisfied with equality, i.e. all demands
1.77 +have to be satisfied and all supplies have to be used.
1.78 +
1.79 +If you need the opposite inequalities in the supply/demand constraints
1.80 +(i.e. the total demand is less than the total supply and all the demands
1.81 +have to be satisfied while there could be supplies that are not used),
1.82 +then you could easily transform the problem to the above form by reversing
1.83 +the direction of the arcs and taking the negative of the supply values
1.84 +(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
1.85 +However \ref NetworkSimplex algorithm also supports this form directly
1.86 +for the sake of convenience.
1.87 +
1.88 +A feasible solution for this problem can be found using \ref Circulation.
1.89 +
1.90 +Note that the above formulation is actually more general than the usual
1.91 +definition of the minimum cost flow problem, in which strict equalities
1.92 +are required in the supply/demand contraints, i.e.
1.93 +
1.94 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
1.95 + sup(u) \quad \forall u\in V. \f]
1.96 +
1.97 +However if the sum of the supply values is zero, then these two problems
1.98 +are equivalent. So if you need the equality form, you have to ensure this
1.99 +additional contraint for the algorithms.
1.100 +
1.101 +The dual solution of the minimum cost flow problem is represented by node
1.102 +potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
1.103 +An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
1.104 +is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
1.105 +node potentials the following \e complementary \e slackness optimality
1.106 +conditions hold.
1.107 +
1.108 + - For all \f$uv\in A\f$ arcs:
1.109 + - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
1.110 + - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
1.111 + - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
1.112 + - For all \f$u\in V\f$:
1.113 + - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
1.114 + then \f$\pi(u)=0\f$.
1.115 +
1.116 +Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
1.117 +\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
1.118 +\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
1.119 +
1.120 +All algorithms provide dual solution (node potentials) as well
1.121 +if an optimal flow is found.
1.122 +
1.123 +LEMON contains several algorithms for solving minimum cost flow problems.
1.124 + - \ref NetworkSimplex Primal Network Simplex algorithm with various
1.125 + pivot strategies.
1.126 + - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
1.127 + cost scaling.
1.128 + - \ref CapacityScaling Successive Shortest %Path algorithm with optional
1.129 capacity scaling.
1.130 - - \ref CostScaling Push-relabel and augment-relabel algorithms based on
1.131 - cost scaling.
1.132 - - \ref NetworkSimplex Primal network simplex algorithm with various
1.133 - pivot strategies.
1.134 + - \ref CancelAndTighten The Cancel and Tighten algorithm.
1.135 + - \ref CycleCanceling Cycle-Canceling algorithms.
1.136 +
1.137 +Most of these implementations support the general inequality form of the
1.138 +minimum cost flow problem, but CancelAndTighten and CycleCanceling
1.139 +only support the equality form due to the primal method they use.
1.140 +
1.141 +In general NetworkSimplex is the most efficient implementation,
1.142 +but in special cases other algorithms could be faster.
1.143 +For example, if the total supply and/or capacities are rather small,
1.144 +CapacityScaling is usually the fastest algorithm (without effective scaling).
1.145 */
1.146
1.147 /**
2.1 --- a/lemon/network_simplex.h Fri Apr 03 18:59:15 2009 +0200
2.2 +++ b/lemon/network_simplex.h Fri Apr 17 18:04:36 2009 +0200
2.3 @@ -30,6 +30,9 @@
2.4
2.5 #include <lemon/core.h>
2.6 #include <lemon/math.h>
2.7 +#include <lemon/maps.h>
2.8 +#include <lemon/circulation.h>
2.9 +#include <lemon/adaptors.h>
2.10
2.11 namespace lemon {
2.12
2.13 @@ -47,6 +50,8 @@
2.14 ///
2.15 /// In general this class is the fastest implementation available
2.16 /// in LEMON for the minimum cost flow problem.
2.17 + /// Moreover it supports both direction of the supply/demand inequality
2.18 + /// constraints. For more information see \ref ProblemType.
2.19 ///
2.20 /// \tparam GR The digraph type the algorithm runs on.
2.21 /// \tparam F The value type used for flow amounts, capacity bounds
2.22 @@ -58,7 +63,8 @@
2.23 /// be integer.
2.24 ///
2.25 /// \note %NetworkSimplex provides five different pivot rule
2.26 - /// implementations. For more information see \ref PivotRule.
2.27 + /// implementations, from which the most efficient one is used
2.28 + /// by default. For more information see \ref PivotRule.
2.29 template <typename GR, typename F = int, typename C = F>
2.30 class NetworkSimplex
2.31 {
2.32 @@ -68,10 +74,17 @@
2.33 typedef F Flow;
2.34 /// The cost type of the algorithm
2.35 typedef C Cost;
2.36 +#ifdef DOXYGEN
2.37 + /// The type of the flow map
2.38 + typedef GR::ArcMap<Flow> FlowMap;
2.39 + /// The type of the potential map
2.40 + typedef GR::NodeMap<Cost> PotentialMap;
2.41 +#else
2.42 /// The type of the flow map
2.43 typedef typename GR::template ArcMap<Flow> FlowMap;
2.44 /// The type of the potential map
2.45 typedef typename GR::template NodeMap<Cost> PotentialMap;
2.46 +#endif
2.47
2.48 public:
2.49
2.50 @@ -117,6 +130,58 @@
2.51 /// candidate list and extends this list in every iteration.
2.52 ALTERING_LIST
2.53 };
2.54 +
2.55 + /// \brief Enum type for selecting the problem type.
2.56 + ///
2.57 + /// Enum type for selecting the problem type, i.e. the direction of
2.58 + /// the inequalities in the supply/demand constraints of the
2.59 + /// \ref min_cost_flow "minimum cost flow problem".
2.60 + ///
2.61 + /// The default problem type is \c GEQ, since this form is supported
2.62 + /// by other minimum cost flow algorithms and the \ref Circulation
2.63 + /// algorithm as well.
2.64 + /// The \c LEQ problem type can be selected using the \ref problemType()
2.65 + /// function.
2.66 + ///
2.67 + /// Note that the equality form is a special case of both problem type.
2.68 + enum ProblemType {
2.69 +
2.70 + /// This option means that there are "<em>greater or equal</em>"
2.71 + /// constraints in the defintion, i.e. the exact formulation of the
2.72 + /// problem is the following.
2.73 + /**
2.74 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
2.75 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
2.76 + sup(u) \quad \forall u\in V \f]
2.77 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
2.78 + */
2.79 + /// It means that the total demand must be greater or equal to the
2.80 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
2.81 + /// negative) and all the supplies have to be carried out from
2.82 + /// the supply nodes, but there could be demands that are not
2.83 + /// satisfied.
2.84 + GEQ,
2.85 + /// It is just an alias for the \c GEQ option.
2.86 + CARRY_SUPPLIES = GEQ,
2.87 +
2.88 + /// This option means that there are "<em>less or equal</em>"
2.89 + /// constraints in the defintion, i.e. the exact formulation of the
2.90 + /// problem is the following.
2.91 + /**
2.92 + \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
2.93 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
2.94 + sup(u) \quad \forall u\in V \f]
2.95 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
2.96 + */
2.97 + /// It means that the total demand must be less or equal to the
2.98 + /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
2.99 + /// positive) and all the demands have to be satisfied, but there
2.100 + /// could be supplies that are not carried out from the supply
2.101 + /// nodes.
2.102 + LEQ,
2.103 + /// It is just an alias for the \c LEQ option.
2.104 + SATISFY_DEMANDS = LEQ
2.105 + };
2.106
2.107 private:
2.108
2.109 @@ -155,6 +220,7 @@
2.110 bool _pstsup;
2.111 Node _psource, _ptarget;
2.112 Flow _pstflow;
2.113 + ProblemType _ptype;
2.114
2.115 // Result maps
2.116 FlowMap *_flow_map;
2.117 @@ -586,13 +652,13 @@
2.118
2.119 /// \brief Constructor.
2.120 ///
2.121 - /// Constructor.
2.122 + /// The constructor of the class.
2.123 ///
2.124 /// \param graph The digraph the algorithm runs on.
2.125 NetworkSimplex(const GR& graph) :
2.126 _graph(graph),
2.127 _plower(NULL), _pupper(NULL), _pcost(NULL),
2.128 - _psupply(NULL), _pstsup(false),
2.129 + _psupply(NULL), _pstsup(false), _ptype(GEQ),
2.130 _flow_map(NULL), _potential_map(NULL),
2.131 _local_flow(false), _local_potential(false),
2.132 _node_id(graph)
2.133 @@ -611,6 +677,12 @@
2.134 if (_local_potential) delete _potential_map;
2.135 }
2.136
2.137 + /// \name Parameters
2.138 + /// The parameters of the algorithm can be specified using these
2.139 + /// functions.
2.140 +
2.141 + /// @{
2.142 +
2.143 /// \brief Set the lower bounds on the arcs.
2.144 ///
2.145 /// This function sets the lower bounds on the arcs.
2.146 @@ -760,6 +832,20 @@
2.147 _pstflow = k;
2.148 return *this;
2.149 }
2.150 +
2.151 + /// \brief Set the problem type.
2.152 + ///
2.153 + /// This function sets the problem type for the algorithm.
2.154 + /// If it is not used before calling \ref run(), the \ref GEQ problem
2.155 + /// type will be used.
2.156 + ///
2.157 + /// For more information see \ref ProblemType.
2.158 + ///
2.159 + /// \return <tt>(*this)</tt>
2.160 + NetworkSimplex& problemType(ProblemType problem_type) {
2.161 + _ptype = problem_type;
2.162 + return *this;
2.163 + }
2.164
2.165 /// \brief Set the flow map.
2.166 ///
2.167 @@ -795,6 +881,8 @@
2.168 _potential_map = ↦
2.169 return *this;
2.170 }
2.171 +
2.172 + /// @}
2.173
2.174 /// \name Execution Control
2.175 /// The algorithm can be executed using \ref run().
2.176 @@ -804,10 +892,11 @@
2.177 /// \brief Run the algorithm.
2.178 ///
2.179 /// This function runs the algorithm.
2.180 - /// The paramters can be specified using \ref lowerMap(),
2.181 + /// The paramters can be specified using functions \ref lowerMap(),
2.182 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
2.183 - /// \ref costMap(), \ref supplyMap() and \ref stSupply()
2.184 - /// functions. For example,
2.185 + /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
2.186 + /// \ref problemType(), \ref flowMap() and \ref potentialMap().
2.187 + /// For example,
2.188 /// \code
2.189 /// NetworkSimplex<ListDigraph> ns(graph);
2.190 /// ns.boundMaps(lower, upper).costMap(cost)
2.191 @@ -830,9 +919,10 @@
2.192 /// \brief Reset all the parameters that have been given before.
2.193 ///
2.194 /// This function resets all the paramaters that have been given
2.195 - /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
2.196 - /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
2.197 - /// \ref stSupply() functions before.
2.198 + /// before using functions \ref lowerMap(), \ref upperMap(),
2.199 + /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
2.200 + /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
2.201 + /// \ref flowMap() and \ref potentialMap().
2.202 ///
2.203 /// It is useful for multiple run() calls. If this function is not
2.204 /// used, all the parameters given before are kept for the next
2.205 @@ -869,6 +959,14 @@
2.206 _pcost = NULL;
2.207 _psupply = NULL;
2.208 _pstsup = false;
2.209 + _ptype = GEQ;
2.210 + if (_local_flow) delete _flow_map;
2.211 + if (_local_potential) delete _potential_map;
2.212 + _flow_map = NULL;
2.213 + _potential_map = NULL;
2.214 + _local_flow = false;
2.215 + _local_potential = false;
2.216 +
2.217 return *this;
2.218 }
2.219
2.220 @@ -1000,20 +1098,21 @@
2.221
2.222 // Initialize node related data
2.223 bool valid_supply = true;
2.224 + Flow sum_supply = 0;
2.225 if (!_pstsup && !_psupply) {
2.226 _pstsup = true;
2.227 _psource = _ptarget = NodeIt(_graph);
2.228 _pstflow = 0;
2.229 }
2.230 if (_psupply) {
2.231 - Flow sum = 0;
2.232 int i = 0;
2.233 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
2.234 _node_id[n] = i;
2.235 _supply[i] = (*_psupply)[n];
2.236 - sum += _supply[i];
2.237 + sum_supply += _supply[i];
2.238 }
2.239 - valid_supply = (sum == 0);
2.240 + valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
2.241 + (_ptype == LEQ && sum_supply >= 0);
2.242 } else {
2.243 int i = 0;
2.244 for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
2.245 @@ -1021,10 +1120,95 @@
2.246 _supply[i] = 0;
2.247 }
2.248 _supply[_node_id[_psource]] = _pstflow;
2.249 - _supply[_node_id[_ptarget]] = -_pstflow;
2.250 + _supply[_node_id[_ptarget]] = -_pstflow;
2.251 }
2.252 if (!valid_supply) return false;
2.253
2.254 + // Infinite capacity value
2.255 + Flow inf_cap =
2.256 + std::numeric_limits<Flow>::has_infinity ?
2.257 + std::numeric_limits<Flow>::infinity() :
2.258 + std::numeric_limits<Flow>::max();
2.259 +
2.260 + // Initialize artifical cost
2.261 + Cost art_cost;
2.262 + if (std::numeric_limits<Cost>::is_exact) {
2.263 + art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
2.264 + } else {
2.265 + art_cost = std::numeric_limits<Cost>::min();
2.266 + for (int i = 0; i != _arc_num; ++i) {
2.267 + if (_cost[i] > art_cost) art_cost = _cost[i];
2.268 + }
2.269 + art_cost = (art_cost + 1) * _node_num;
2.270 + }
2.271 +
2.272 + // Run Circulation to check if a feasible solution exists
2.273 + typedef ConstMap<Arc, Flow> ConstArcMap;
2.274 + FlowNodeMap *csup = NULL;
2.275 + bool local_csup = false;
2.276 + if (_psupply) {
2.277 + csup = _psupply;
2.278 + } else {
2.279 + csup = new FlowNodeMap(_graph, 0);
2.280 + (*csup)[_psource] = _pstflow;
2.281 + (*csup)[_ptarget] = -_pstflow;
2.282 + local_csup = true;
2.283 + }
2.284 + bool circ_result = false;
2.285 + if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
2.286 + // GEQ problem type
2.287 + if (_plower) {
2.288 + if (_pupper) {
2.289 + Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
2.290 + circ(_graph, *_plower, *_pupper, *csup);
2.291 + circ_result = circ.run();
2.292 + } else {
2.293 + Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
2.294 + circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
2.295 + circ_result = circ.run();
2.296 + }
2.297 + } else {
2.298 + if (_pupper) {
2.299 + Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
2.300 + circ(_graph, ConstArcMap(0), *_pupper, *csup);
2.301 + circ_result = circ.run();
2.302 + } else {
2.303 + Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
2.304 + circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
2.305 + circ_result = circ.run();
2.306 + }
2.307 + }
2.308 + } else {
2.309 + // LEQ problem type
2.310 + typedef ReverseDigraph<const GR> RevGraph;
2.311 + typedef NegMap<FlowNodeMap> NegNodeMap;
2.312 + RevGraph rgraph(_graph);
2.313 + NegNodeMap neg_csup(*csup);
2.314 + if (_plower) {
2.315 + if (_pupper) {
2.316 + Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
2.317 + circ(rgraph, *_plower, *_pupper, neg_csup);
2.318 + circ_result = circ.run();
2.319 + } else {
2.320 + Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
2.321 + circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
2.322 + circ_result = circ.run();
2.323 + }
2.324 + } else {
2.325 + if (_pupper) {
2.326 + Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
2.327 + circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
2.328 + circ_result = circ.run();
2.329 + } else {
2.330 + Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
2.331 + circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
2.332 + circ_result = circ.run();
2.333 + }
2.334 + }
2.335 + }
2.336 + if (local_csup) delete csup;
2.337 + if (!circ_result) return false;
2.338 +
2.339 // Set data for the artificial root node
2.340 _root = _node_num;
2.341 _parent[_root] = -1;
2.342 @@ -1033,8 +1217,12 @@
2.343 _rev_thread[0] = _root;
2.344 _succ_num[_root] = all_node_num;
2.345 _last_succ[_root] = _root - 1;
2.346 - _supply[_root] = 0;
2.347 - _pi[_root] = 0;
2.348 + _supply[_root] = -sum_supply;
2.349 + if (sum_supply < 0) {
2.350 + _pi[_root] = -art_cost;
2.351 + } else {
2.352 + _pi[_root] = art_cost;
2.353 + }
2.354
2.355 // Store the arcs in a mixed order
2.356 int k = std::max(int(sqrt(_arc_num)), 10);
2.357 @@ -1045,10 +1233,6 @@
2.358 }
2.359
2.360 // Initialize arc maps
2.361 - Flow inf_cap =
2.362 - std::numeric_limits<Flow>::has_infinity ?
2.363 - std::numeric_limits<Flow>::infinity() :
2.364 - std::numeric_limits<Flow>::max();
2.365 if (_pupper && _pcost) {
2.366 for (int i = 0; i != _arc_num; ++i) {
2.367 Arc e = _arc_ref[i];
2.368 @@ -1083,18 +1267,6 @@
2.369 }
2.370 }
2.371
2.372 - // Initialize artifical cost
2.373 - Cost art_cost;
2.374 - if (std::numeric_limits<Cost>::is_exact) {
2.375 - art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
2.376 - } else {
2.377 - art_cost = std::numeric_limits<Cost>::min();
2.378 - for (int i = 0; i != _arc_num; ++i) {
2.379 - if (_cost[i] > art_cost) art_cost = _cost[i];
2.380 - }
2.381 - art_cost = (art_cost + 1) * _node_num;
2.382 - }
2.383 -
2.384 // Remove non-zero lower bounds
2.385 if (_plower) {
2.386 for (int i = 0; i != _arc_num; ++i) {
2.387 @@ -1118,14 +1290,14 @@
2.388 _cost[e] = art_cost;
2.389 _cap[e] = inf_cap;
2.390 _state[e] = STATE_TREE;
2.391 - if (_supply[u] >= 0) {
2.392 + if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
2.393 _flow[e] = _supply[u];
2.394 _forward[u] = true;
2.395 - _pi[u] = -art_cost;
2.396 + _pi[u] = -art_cost + _pi[_root];
2.397 } else {
2.398 _flow[e] = -_supply[u];
2.399 _forward[u] = false;
2.400 - _pi[u] = art_cost;
2.401 + _pi[u] = art_cost + _pi[_root];
2.402 }
2.403 }
2.404
2.405 @@ -1382,11 +1554,6 @@
2.406 }
2.407 }
2.408
2.409 - // Check if the flow amount equals zero on all the artificial arcs
2.410 - for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
2.411 - if (_flow[e] > 0) return false;
2.412 - }
2.413 -
2.414 // Copy flow values to _flow_map
2.415 if (_plower) {
2.416 for (int i = 0; i != _arc_num; ++i) {
3.1 --- a/test/min_cost_flow_test.cc Fri Apr 03 18:59:15 2009 +0200
3.2 +++ b/test/min_cost_flow_test.cc Fri Apr 17 18:04:36 2009 +0200
3.3 @@ -33,19 +33,19 @@
3.4
3.5 char test_lgf[] =
3.6 "@nodes\n"
3.7 - "label sup1 sup2 sup3\n"
3.8 - " 1 20 27 0\n"
3.9 - " 2 -4 0 0\n"
3.10 - " 3 0 0 0\n"
3.11 - " 4 0 0 0\n"
3.12 - " 5 9 0 0\n"
3.13 - " 6 -6 0 0\n"
3.14 - " 7 0 0 0\n"
3.15 - " 8 0 0 0\n"
3.16 - " 9 3 0 0\n"
3.17 - " 10 -2 0 0\n"
3.18 - " 11 0 0 0\n"
3.19 - " 12 -20 -27 0\n"
3.20 + "label sup1 sup2 sup3 sup4 sup5\n"
3.21 + " 1 20 27 0 20 30\n"
3.22 + " 2 -4 0 0 -8 -3\n"
3.23 + " 3 0 0 0 0 0\n"
3.24 + " 4 0 0 0 0 0\n"
3.25 + " 5 9 0 0 6 11\n"
3.26 + " 6 -6 0 0 -5 -6\n"
3.27 + " 7 0 0 0 0 0\n"
3.28 + " 8 0 0 0 0 3\n"
3.29 + " 9 3 0 0 0 0\n"
3.30 + " 10 -2 0 0 -7 -2\n"
3.31 + " 11 0 0 0 -10 0\n"
3.32 + " 12 -20 -27 0 -30 -20\n"
3.33 "\n"
3.34 "@arcs\n"
3.35 " cost cap low1 low2\n"
3.36 @@ -76,6 +76,12 @@
3.37 "target 12\n";
3.38
3.39
3.40 +enum ProblemType {
3.41 + EQ,
3.42 + GEQ,
3.43 + LEQ
3.44 +};
3.45 +
3.46 // Check the interface of an MCF algorithm
3.47 template <typename GR, typename Flow, typename Cost>
3.48 class McfClassConcept
3.49 @@ -97,17 +103,19 @@
3.50 .costMap(cost)
3.51 .supplyMap(sup)
3.52 .stSupply(n, n, k)
3.53 + .flowMap(flow)
3.54 + .potentialMap(pot)
3.55 .run();
3.56 +
3.57 + const MCF& const_mcf = mcf;
3.58
3.59 - const typename MCF::FlowMap &fm = mcf.flowMap();
3.60 - const typename MCF::PotentialMap &pm = mcf.potentialMap();
3.61 + const typename MCF::FlowMap &fm = const_mcf.flowMap();
3.62 + const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
3.63
3.64 - v = mcf.totalCost();
3.65 - double x = mcf.template totalCost<double>();
3.66 - v = mcf.flow(a);
3.67 - v = mcf.potential(n);
3.68 - mcf.flowMap(flow);
3.69 - mcf.potentialMap(pot);
3.70 + v = const_mcf.totalCost();
3.71 + double x = const_mcf.template totalCost<double>();
3.72 + v = const_mcf.flow(a);
3.73 + v = const_mcf.potential(n);
3.74
3.75 ignore_unused_variable_warning(fm);
3.76 ignore_unused_variable_warning(pm);
3.77 @@ -142,7 +150,8 @@
3.78 template < typename GR, typename LM, typename UM,
3.79 typename SM, typename FM >
3.80 bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
3.81 - const SM& supply, const FM& flow )
3.82 + const SM& supply, const FM& flow,
3.83 + ProblemType type = EQ )
3.84 {
3.85 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
3.86
3.87 @@ -156,7 +165,10 @@
3.88 sum += flow[e];
3.89 for (InArcIt e(gr, n); e != INVALID; ++e)
3.90 sum -= flow[e];
3.91 - if (sum != supply[n]) return false;
3.92 + bool b = (type == EQ && sum == supply[n]) ||
3.93 + (type == GEQ && sum >= supply[n]) ||
3.94 + (type == LEQ && sum <= supply[n]);
3.95 + if (!b) return false;
3.96 }
3.97
3.98 return true;
3.99 @@ -165,9 +177,10 @@
3.100 // Check the feasibility of the given potentials (dual soluiton)
3.101 // using the "Complementary Slackness" optimality condition
3.102 template < typename GR, typename LM, typename UM,
3.103 - typename CM, typename FM, typename PM >
3.104 + typename CM, typename SM, typename FM, typename PM >
3.105 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
3.106 - const CM& cost, const FM& flow, const PM& pi )
3.107 + const CM& cost, const SM& supply, const FM& flow,
3.108 + const PM& pi )
3.109 {
3.110 TEMPLATE_DIGRAPH_TYPEDEFS(GR);
3.111
3.112 @@ -179,6 +192,16 @@
3.113 (red_cost > 0 && flow[e] == lower[e]) ||
3.114 (red_cost < 0 && flow[e] == upper[e]);
3.115 }
3.116 +
3.117 + for (NodeIt n(gr); opt && n != INVALID; ++n) {
3.118 + typename SM::Value sum = 0;
3.119 + for (OutArcIt e(gr, n); e != INVALID; ++e)
3.120 + sum += flow[e];
3.121 + for (InArcIt e(gr, n); e != INVALID; ++e)
3.122 + sum -= flow[e];
3.123 + opt = (sum == supply[n]) || (pi[n] == 0);
3.124 + }
3.125 +
3.126 return opt;
3.127 }
3.128
3.129 @@ -190,14 +213,15 @@
3.130 const GR& gr, const LM& lower, const UM& upper,
3.131 const CM& cost, const SM& supply,
3.132 bool result, typename CM::Value total,
3.133 - const std::string &test_id = "" )
3.134 + const std::string &test_id = "",
3.135 + ProblemType type = EQ )
3.136 {
3.137 check(mcf_result == result, "Wrong result " + test_id);
3.138 if (result) {
3.139 - check(checkFlow(gr, lower, upper, supply, mcf.flowMap()),
3.140 + check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
3.141 "The flow is not feasible " + test_id);
3.142 check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
3.143 - check(checkPotential(gr, lower, upper, cost, mcf.flowMap(),
3.144 + check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
3.145 mcf.potentialMap()),
3.146 "Wrong potentials " + test_id);
3.147 }
3.148 @@ -226,7 +250,7 @@
3.149 // Read the test digraph
3.150 Digraph gr;
3.151 Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
3.152 - Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr);
3.153 + Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
3.154 ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
3.155 Node v, w;
3.156
3.157 @@ -239,6 +263,8 @@
3.158 .nodeMap("sup1", s1)
3.159 .nodeMap("sup2", s2)
3.160 .nodeMap("sup3", s3)
3.161 + .nodeMap("sup4", s4)
3.162 + .nodeMap("sup5", s5)
3.163 .node("source", v)
3.164 .node("target", w)
3.165 .run();
3.166 @@ -247,6 +273,7 @@
3.167 {
3.168 NetworkSimplex<Digraph> mcf(gr);
3.169
3.170 + // Check the equality form
3.171 mcf.upperMap(u).costMap(c);
3.172 checkMcf(mcf, mcf.supplyMap(s1).run(),
3.173 gr, l1, u, c, s1, true, 5240, "#A1");
3.174 @@ -267,6 +294,28 @@
3.175 gr, l1, cu, cc, s3, true, 0, "#A7");
3.176 checkMcf(mcf, mcf.boundMaps(l2, u).run(),
3.177 gr, l2, u, cc, s3, false, 0, "#A8");
3.178 +
3.179 + // Check the GEQ form
3.180 + mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
3.181 + checkMcf(mcf, mcf.run(),
3.182 + gr, l1, u, c, s4, true, 3530, "#A9", GEQ);
3.183 + mcf.problemType(mcf.GEQ);
3.184 + checkMcf(mcf, mcf.lowerMap(l2).run(),
3.185 + gr, l2, u, c, s4, true, 4540, "#A10", GEQ);
3.186 + mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
3.187 + checkMcf(mcf, mcf.run(),
3.188 + gr, l2, u, c, s5, false, 0, "#A11", GEQ);
3.189 +
3.190 + // Check the LEQ form
3.191 + mcf.reset().problemType(mcf.LEQ);
3.192 + mcf.upperMap(u).costMap(c).supplyMap(s5);
3.193 + checkMcf(mcf, mcf.run(),
3.194 + gr, l1, u, c, s5, true, 5080, "#A12", LEQ);
3.195 + checkMcf(mcf, mcf.lowerMap(l2).run(),
3.196 + gr, l2, u, c, s5, true, 5930, "#A13", LEQ);
3.197 + mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
3.198 + checkMcf(mcf, mcf.run(),
3.199 + gr, l2, u, c, s4, false, 0, "#A14", LEQ);
3.200 }
3.201
3.202 // B. Test NetworkSimplex with each pivot rule