0
4
0
10
10
202
243
... | ... |
@@ -352,17 +352,17 @@ |
352 | 352 |
minimum total cost from a set of supply nodes to a set of demand nodes |
353 | 353 |
in a network with capacity constraints (lower and upper bounds) |
354 | 354 |
and arc costs. |
355 |
Formally, let \f$G=(V,A)\f$ be a digraph, |
|
356 |
\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and |
|
355 |
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$, |
|
356 |
\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and |
|
357 | 357 |
upper bounds for the flow values on the arcs, for which |
358 |
\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$. |
|
359 |
\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow |
|
360 |
|
|
358 |
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, |
|
359 |
\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow |
|
360 |
on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the |
|
361 | 361 |
signed supply values of the nodes. |
362 | 362 |
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ |
363 | 363 |
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with |
364 | 364 |
\f$-sup(u)\f$ demand. |
365 |
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z} |
|
365 |
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution |
|
366 | 366 |
of the following optimization problem. |
367 | 367 |
|
368 | 368 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
... | ... |
@@ -404,7 +404,7 @@ |
404 | 404 |
|
405 | 405 |
The dual solution of the minimum cost flow problem is represented by node |
406 | 406 |
potentials \f$\pi: V\rightarrow\mathbf{Z}\f$. |
407 |
An \f$f: A\rightarrow\mathbf{Z} |
|
407 |
An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem |
|
408 | 408 |
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$ |
409 | 409 |
node potentials the following \e complementary \e slackness optimality |
410 | 410 |
conditions hold. |
... | ... |
@@ -413,15 +413,15 @@ |
413 | 413 |
- if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; |
414 | 414 |
- if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; |
415 | 415 |
- if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. |
416 |
- For all \f$u\in V\f$: |
|
416 |
- For all \f$u\in V\f$ nodes: |
|
417 | 417 |
- if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$, |
418 | 418 |
then \f$\pi(u)=0\f$. |
419 | 419 |
|
420 | 420 |
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc |
421 |
\f$uv\in A\f$ with respect to the |
|
421 |
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e. |
|
422 | 422 |
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f] |
423 | 423 |
|
424 |
All algorithms provide dual solution (node potentials) as well |
|
424 |
All algorithms provide dual solution (node potentials) as well, |
|
425 | 425 |
if an optimal flow is found. |
426 | 426 |
|
427 | 427 |
LEMON contains several algorithms for solving minimum cost flow problems. |
... | ... |
@@ -30,9 +30,6 @@ |
30 | 30 |
|
31 | 31 |
#include <lemon/core.h> |
32 | 32 |
#include <lemon/math.h> |
33 |
#include <lemon/maps.h> |
|
34 |
#include <lemon/circulation.h> |
|
35 |
#include <lemon/adaptors.h> |
|
36 | 33 |
|
37 | 34 |
namespace lemon { |
38 | 35 |
|
... | ... |
@@ -50,8 +47,13 @@ |
50 | 47 |
/// |
51 | 48 |
/// In general this class is the fastest implementation available |
52 | 49 |
/// in LEMON for the minimum cost flow problem. |
53 |
/// Moreover it supports both direction of the supply/demand inequality |
|
54 |
/// constraints. For more information see \ref ProblemType. |
|
50 |
/// Moreover it supports both directions of the supply/demand inequality |
|
51 |
/// constraints. For more information see \ref SupplyType. |
|
52 |
/// |
|
53 |
/// Most of the parameters of the problem (except for the digraph) |
|
54 |
/// can be given using separate functions, and the algorithm can be |
|
55 |
/// executed using the \ref run() function. If some parameters are not |
|
56 |
/// specified, then default values will be used. |
|
55 | 57 |
/// |
56 | 58 |
/// \tparam GR The digraph type the algorithm runs on. |
57 | 59 |
/// \tparam F The value type used for flow amounts, capacity bounds |
... | ... |
@@ -88,11 +90,80 @@ |
88 | 90 |
|
89 | 91 |
public: |
90 | 92 |
|
91 |
/// \brief |
|
93 |
/// \brief Problem type constants for the \c run() function. |
|
92 | 94 |
/// |
93 |
/// Enum type |
|
95 |
/// Enum type containing the problem type constants that can be |
|
96 |
/// returned by the \ref run() function of the algorithm. |
|
97 |
enum ProblemType { |
|
98 |
/// The problem has no feasible solution (flow). |
|
99 |
INFEASIBLE, |
|
100 |
/// The problem has optimal solution (i.e. it is feasible and |
|
101 |
/// bounded), and the algorithm has found optimal flow and node |
|
102 |
/// potentials (primal and dual solutions). |
|
103 |
OPTIMAL, |
|
104 |
/// The objective function of the problem is unbounded, i.e. |
|
105 |
/// there is a directed cycle having negative total cost and |
|
106 |
/// infinite upper bound. |
|
107 |
UNBOUNDED |
|
108 |
}; |
|
109 |
|
|
110 |
/// \brief Constants for selecting the type of the supply constraints. |
|
111 |
/// |
|
112 |
/// Enum type containing constants for selecting the supply type, |
|
113 |
/// i.e. the direction of the inequalities in the supply/demand |
|
114 |
/// constraints of the \ref min_cost_flow "minimum cost flow problem". |
|
115 |
/// |
|
116 |
/// The default supply type is \c GEQ, since this form is supported |
|
117 |
/// by other minimum cost flow algorithms and the \ref Circulation |
|
118 |
/// algorithm, as well. |
|
119 |
/// The \c LEQ problem type can be selected using the \ref supplyType() |
|
94 | 120 |
/// function. |
95 | 121 |
/// |
122 |
/// Note that the equality form is a special case of both supply types. |
|
123 |
enum SupplyType { |
|
124 |
|
|
125 |
/// This option means that there are <em>"greater or equal"</em> |
|
126 |
/// supply/demand constraints in the definition, i.e. the exact |
|
127 |
/// formulation of the problem is the following. |
|
128 |
/** |
|
129 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
130 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq |
|
131 |
sup(u) \quad \forall u\in V \f] |
|
132 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
133 |
*/ |
|
134 |
/// It means that the total demand must be greater or equal to the |
|
135 |
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or |
|
136 |
/// negative) and all the supplies have to be carried out from |
|
137 |
/// the supply nodes, but there could be demands that are not |
|
138 |
/// satisfied. |
|
139 |
GEQ, |
|
140 |
/// It is just an alias for the \c GEQ option. |
|
141 |
CARRY_SUPPLIES = GEQ, |
|
142 |
|
|
143 |
/// This option means that there are <em>"less or equal"</em> |
|
144 |
/// supply/demand constraints in the definition, i.e. the exact |
|
145 |
/// formulation of the problem is the following. |
|
146 |
/** |
|
147 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
148 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq |
|
149 |
sup(u) \quad \forall u\in V \f] |
|
150 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
151 |
*/ |
|
152 |
/// It means that the total demand must be less or equal to the |
|
153 |
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or |
|
154 |
/// positive) and all the demands have to be satisfied, but there |
|
155 |
/// could be supplies that are not carried out from the supply |
|
156 |
/// nodes. |
|
157 |
LEQ, |
|
158 |
/// It is just an alias for the \c LEQ option. |
|
159 |
SATISFY_DEMANDS = LEQ |
|
160 |
}; |
|
161 |
|
|
162 |
/// \brief Constants for selecting the pivot rule. |
|
163 |
/// |
|
164 |
/// Enum type containing constants for selecting the pivot rule for |
|
165 |
/// the \ref run() function. |
|
166 |
/// |
|
96 | 167 |
/// \ref NetworkSimplex provides five different pivot rule |
97 | 168 |
/// implementations that significantly affect the running time |
98 | 169 |
/// of the algorithm. |
... | ... |
@@ -131,58 +202,6 @@ |
131 | 202 |
ALTERING_LIST |
132 | 203 |
}; |
133 | 204 |
|
134 |
/// \brief Enum type for selecting the problem type. |
|
135 |
/// |
|
136 |
/// Enum type for selecting the problem type, i.e. the direction of |
|
137 |
/// the inequalities in the supply/demand constraints of the |
|
138 |
/// \ref min_cost_flow "minimum cost flow problem". |
|
139 |
/// |
|
140 |
/// The default problem type is \c GEQ, since this form is supported |
|
141 |
/// by other minimum cost flow algorithms and the \ref Circulation |
|
142 |
/// algorithm as well. |
|
143 |
/// The \c LEQ problem type can be selected using the \ref problemType() |
|
144 |
/// function. |
|
145 |
/// |
|
146 |
/// Note that the equality form is a special case of both problem type. |
|
147 |
enum ProblemType { |
|
148 |
|
|
149 |
/// This option means that there are "<em>greater or equal</em>" |
|
150 |
/// constraints in the defintion, i.e. the exact formulation of the |
|
151 |
/// problem is the following. |
|
152 |
/** |
|
153 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
154 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq |
|
155 |
sup(u) \quad \forall u\in V \f] |
|
156 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
157 |
*/ |
|
158 |
/// It means that the total demand must be greater or equal to the |
|
159 |
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or |
|
160 |
/// negative) and all the supplies have to be carried out from |
|
161 |
/// the supply nodes, but there could be demands that are not |
|
162 |
/// satisfied. |
|
163 |
GEQ, |
|
164 |
/// It is just an alias for the \c GEQ option. |
|
165 |
CARRY_SUPPLIES = GEQ, |
|
166 |
|
|
167 |
/// This option means that there are "<em>less or equal</em>" |
|
168 |
/// constraints in the defintion, i.e. the exact formulation of the |
|
169 |
/// problem is the following. |
|
170 |
/** |
|
171 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
172 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq |
|
173 |
sup(u) \quad \forall u\in V \f] |
|
174 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
175 |
*/ |
|
176 |
/// It means that the total demand must be less or equal to the |
|
177 |
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or |
|
178 |
/// positive) and all the demands have to be satisfied, but there |
|
179 |
/// could be supplies that are not carried out from the supply |
|
180 |
/// nodes. |
|
181 |
LEQ, |
|
182 |
/// It is just an alias for the \c LEQ option. |
|
183 |
SATISFY_DEMANDS = LEQ |
|
184 |
}; |
|
185 |
|
|
186 | 205 |
private: |
187 | 206 |
|
188 | 207 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
... | ... |
@@ -220,7 +239,9 @@ |
220 | 239 |
bool _pstsup; |
221 | 240 |
Node _psource, _ptarget; |
222 | 241 |
Flow _pstflow; |
223 |
|
|
242 |
SupplyType _stype; |
|
243 |
|
|
244 |
Flow _sum_supply; |
|
224 | 245 |
|
225 | 246 |
// Result maps |
226 | 247 |
FlowMap *_flow_map; |
... | ... |
@@ -259,6 +280,15 @@ |
259 | 280 |
int stem, par_stem, new_stem; |
260 | 281 |
Flow delta; |
261 | 282 |
|
283 |
public: |
|
284 |
|
|
285 |
/// \brief Constant for infinite upper bounds (capacities). |
|
286 |
/// |
|
287 |
/// Constant for infinite upper bounds (capacities). |
|
288 |
/// It is \c std::numeric_limits<Flow>::infinity() if available, |
|
289 |
/// \c std::numeric_limits<Flow>::max() otherwise. |
|
290 |
const Flow INF; |
|
291 |
|
|
262 | 292 |
private: |
263 | 293 |
|
264 | 294 |
// Implementation of the First Eligible pivot rule |
... | ... |
@@ -661,17 +691,19 @@ |
661 | 691 |
NetworkSimplex(const GR& graph) : |
662 | 692 |
_graph(graph), |
663 | 693 |
_plower(NULL), _pupper(NULL), _pcost(NULL), |
664 |
_psupply(NULL), _pstsup(false), |
|
694 |
_psupply(NULL), _pstsup(false), _stype(GEQ), |
|
665 | 695 |
_flow_map(NULL), _potential_map(NULL), |
666 | 696 |
_local_flow(false), _local_potential(false), |
667 |
_node_id(graph) |
|
697 |
_node_id(graph), |
|
698 |
INF(std::numeric_limits<Flow>::has_infinity ? |
|
699 |
std::numeric_limits<Flow>::infinity() : |
|
700 |
std::numeric_limits<Flow>::max()) |
|
668 | 701 |
{ |
669 |
LEMON_ASSERT(std::numeric_limits<Flow>::is_integer && |
|
670 |
std::numeric_limits<Flow>::is_signed, |
|
671 |
"The flow type of NetworkSimplex must be signed integer"); |
|
672 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_integer && |
|
673 |
std::numeric_limits<Cost>::is_signed, |
|
674 |
"The cost type of NetworkSimplex must be signed integer"); |
|
702 |
// Check the value types |
|
703 |
LEMON_ASSERT(std::numeric_limits<Flow>::is_signed, |
|
704 |
"The flow type of NetworkSimplex must be signed"); |
|
705 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
|
706 |
"The cost type of NetworkSimplex must be signed"); |
|
675 | 707 |
} |
676 | 708 |
|
677 | 709 |
/// Destructor. |
... | ... |
@@ -689,17 +721,16 @@ |
689 | 721 |
/// \brief Set the lower bounds on the arcs. |
690 | 722 |
/// |
691 | 723 |
/// This function sets the lower bounds on the arcs. |
692 |
/// If neither this function nor \ref boundMaps() is used before |
|
693 |
/// calling \ref run(), the lower bounds will be set to zero |
|
694 |
/// |
|
724 |
/// If it is not used before calling \ref run(), the lower bounds |
|
725 |
/// will be set to zero on all arcs. |
|
695 | 726 |
/// |
696 | 727 |
/// \param map An arc map storing the lower bounds. |
697 | 728 |
/// Its \c Value type must be convertible to the \c Flow type |
698 | 729 |
/// of the algorithm. |
699 | 730 |
/// |
700 | 731 |
/// \return <tt>(*this)</tt> |
701 |
template <typename LOWER> |
|
702 |
NetworkSimplex& lowerMap(const LOWER& map) { |
|
732 |
template <typename LowerMap> |
|
733 |
NetworkSimplex& lowerMap(const LowerMap& map) { |
|
703 | 734 |
delete _plower; |
704 | 735 |
_plower = new FlowArcMap(_graph); |
705 | 736 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
... | ... |
@@ -711,18 +742,17 @@ |
711 | 742 |
/// \brief Set the upper bounds (capacities) on the arcs. |
712 | 743 |
/// |
713 | 744 |
/// This function sets the upper bounds (capacities) on the arcs. |
714 |
/// If none of the functions \ref upperMap(), \ref capacityMap() |
|
715 |
/// and \ref boundMaps() is used before calling \ref run(), |
|
716 |
/// the upper bounds (capacities) will be set to |
|
717 |
/// \c std::numeric_limits<Flow>::max() on all arcs. |
|
745 |
/// If it is not used before calling \ref run(), the upper bounds |
|
746 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
|
747 |
/// unbounded from above on each arc). |
|
718 | 748 |
/// |
719 | 749 |
/// \param map An arc map storing the upper bounds. |
720 | 750 |
/// Its \c Value type must be convertible to the \c Flow type |
721 | 751 |
/// of the algorithm. |
722 | 752 |
/// |
723 | 753 |
/// \return <tt>(*this)</tt> |
724 |
template<typename UPPER> |
|
725 |
NetworkSimplex& upperMap(const UPPER& map) { |
|
754 |
template<typename UpperMap> |
|
755 |
NetworkSimplex& upperMap(const UpperMap& map) { |
|
726 | 756 |
delete _pupper; |
727 | 757 |
_pupper = new FlowArcMap(_graph); |
728 | 758 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
... | ... |
@@ -731,43 +761,6 @@ |
731 | 761 |
return *this; |
732 | 762 |
} |
733 | 763 |
|
734 |
/// \brief Set the upper bounds (capacities) on the arcs. |
|
735 |
/// |
|
736 |
/// This function sets the upper bounds (capacities) on the arcs. |
|
737 |
/// It is just an alias for \ref upperMap(). |
|
738 |
/// |
|
739 |
/// \return <tt>(*this)</tt> |
|
740 |
template<typename CAP> |
|
741 |
NetworkSimplex& capacityMap(const CAP& map) { |
|
742 |
return upperMap(map); |
|
743 |
} |
|
744 |
|
|
745 |
/// \brief Set the lower and upper bounds on the arcs. |
|
746 |
/// |
|
747 |
/// This function sets the lower and upper bounds on the arcs. |
|
748 |
/// If neither this function nor \ref lowerMap() is used before |
|
749 |
/// calling \ref run(), the lower bounds will be set to zero |
|
750 |
/// on all arcs. |
|
751 |
/// If none of the functions \ref upperMap(), \ref capacityMap() |
|
752 |
/// and \ref boundMaps() is used before calling \ref run(), |
|
753 |
/// the upper bounds (capacities) will be set to |
|
754 |
/// \c std::numeric_limits<Flow>::max() on all arcs. |
|
755 |
/// |
|
756 |
/// \param lower An arc map storing the lower bounds. |
|
757 |
/// \param upper An arc map storing the upper bounds. |
|
758 |
/// |
|
759 |
/// The \c Value type of the maps must be convertible to the |
|
760 |
/// \c Flow type of the algorithm. |
|
761 |
/// |
|
762 |
/// \note This function is just a shortcut of calling \ref lowerMap() |
|
763 |
/// and \ref upperMap() separately. |
|
764 |
/// |
|
765 |
/// \return <tt>(*this)</tt> |
|
766 |
template <typename LOWER, typename UPPER> |
|
767 |
NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) { |
|
768 |
return lowerMap(lower).upperMap(upper); |
|
769 |
} |
|
770 |
|
|
771 | 764 |
/// \brief Set the costs of the arcs. |
772 | 765 |
/// |
773 | 766 |
/// This function sets the costs of the arcs. |
... | ... |
@@ -779,8 +772,8 @@ |
779 | 772 |
/// of the algorithm. |
780 | 773 |
/// |
781 | 774 |
/// \return <tt>(*this)</tt> |
782 |
template<typename COST> |
|
783 |
NetworkSimplex& costMap(const COST& map) { |
|
775 |
template<typename CostMap> |
|
776 |
NetworkSimplex& costMap(const CostMap& map) { |
|
784 | 777 |
delete _pcost; |
785 | 778 |
_pcost = new CostArcMap(_graph); |
786 | 779 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
... | ... |
@@ -801,8 +794,8 @@ |
801 | 794 |
/// of the algorithm. |
802 | 795 |
/// |
803 | 796 |
/// \return <tt>(*this)</tt> |
804 |
template<typename SUP> |
|
805 |
NetworkSimplex& supplyMap(const SUP& map) { |
|
797 |
template<typename SupplyMap> |
|
798 |
NetworkSimplex& supplyMap(const SupplyMap& map) { |
|
806 | 799 |
delete _psupply; |
807 | 800 |
_pstsup = false; |
808 | 801 |
_psupply = new FlowNodeMap(_graph); |
... | ... |
@@ -820,6 +813,10 @@ |
820 | 813 |
/// calling \ref run(), the supply of each node will be set to zero. |
821 | 814 |
/// (It makes sense only if non-zero lower bounds are given.) |
822 | 815 |
/// |
816 |
/// Using this function has the same effect as using \ref supplyMap() |
|
817 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
|
818 |
/// assigned to \c t and all other nodes have zero supply value. |
|
819 |
/// |
|
823 | 820 |
/// \param s The source node. |
824 | 821 |
/// \param t The target node. |
825 | 822 |
/// \param k The required amount of flow from node \c s to node \c t |
... | ... |
@@ -836,17 +833,17 @@ |
836 | 833 |
return *this; |
837 | 834 |
} |
838 | 835 |
|
839 |
/// \brief Set the |
|
836 |
/// \brief Set the type of the supply constraints. |
|
840 | 837 |
/// |
841 |
/// This function sets the problem type for the algorithm. |
|
842 |
/// If it is not used before calling \ref run(), the \ref GEQ problem |
|
838 |
/// This function sets the type of the supply/demand constraints. |
|
839 |
/// If it is not used before calling \ref run(), the \ref GEQ supply |
|
843 | 840 |
/// type will be used. |
844 | 841 |
/// |
845 |
/// For more information see \ref |
|
842 |
/// For more information see \ref SupplyType. |
|
846 | 843 |
/// |
847 | 844 |
/// \return <tt>(*this)</tt> |
848 |
NetworkSimplex& problemType(ProblemType problem_type) { |
|
849 |
_ptype = problem_type; |
|
845 |
NetworkSimplex& supplyType(SupplyType supply_type) { |
|
846 |
_stype = supply_type; |
|
850 | 847 |
return *this; |
851 | 848 |
} |
852 | 849 |
|
... | ... |
@@ -896,13 +893,12 @@ |
896 | 893 |
/// |
897 | 894 |
/// This function runs the algorithm. |
898 | 895 |
/// The paramters can be specified using functions \ref lowerMap(), |
899 |
/// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), |
|
900 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(), |
|
901 |
/// \ref |
|
896 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), |
|
897 |
/// \ref supplyType(), \ref flowMap() and \ref potentialMap(). |
|
902 | 898 |
/// For example, |
903 | 899 |
/// \code |
904 | 900 |
/// NetworkSimplex<ListDigraph> ns(graph); |
905 |
/// ns. |
|
901 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
|
906 | 902 |
/// .supplyMap(sup).run(); |
907 | 903 |
/// \endcode |
908 | 904 |
/// |
... | ... |
@@ -914,17 +910,25 @@ |
914 | 910 |
/// \param pivot_rule The pivot rule that will be used during the |
915 | 911 |
/// algorithm. For more information see \ref PivotRule. |
916 | 912 |
/// |
917 |
/// \return \c true if a feasible flow can be found. |
|
918 |
bool run(PivotRule pivot_rule = BLOCK_SEARCH) { |
|
919 |
|
|
913 |
/// \return \c INFEASIBLE if no feasible flow exists, |
|
914 |
/// \n \c OPTIMAL if the problem has optimal solution |
|
915 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
|
916 |
/// optimal flow and node potentials (primal and dual solutions), |
|
917 |
/// \n \c UNBOUNDED if the objective function of the problem is |
|
918 |
/// unbounded, i.e. there is a directed cycle having negative total |
|
919 |
/// cost and infinite upper bound. |
|
920 |
/// |
|
921 |
/// \see ProblemType, PivotRule |
|
922 |
ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { |
|
923 |
if (!init()) return INFEASIBLE; |
|
924 |
return start(pivot_rule); |
|
920 | 925 |
} |
921 | 926 |
|
922 | 927 |
/// \brief Reset all the parameters that have been given before. |
923 | 928 |
/// |
924 | 929 |
/// This function resets all the paramaters that have been given |
925 | 930 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
926 |
/// \ref capacityMap(), \ref boundMaps(), \ref costMap(), |
|
927 |
/// \ref supplyMap(), \ref stSupply(), \ref problemType(), |
|
931 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(), |
|
928 | 932 |
/// \ref flowMap() and \ref potentialMap(). |
929 | 933 |
/// |
930 | 934 |
/// It is useful for multiple run() calls. If this function is not |
... | ... |
@@ -936,7 +940,7 @@ |
936 | 940 |
/// NetworkSimplex<ListDigraph> ns(graph); |
937 | 941 |
/// |
938 | 942 |
/// // First run |
939 |
/// ns.lowerMap(lower). |
|
943 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
|
940 | 944 |
/// .supplyMap(sup).run(); |
941 | 945 |
/// |
942 | 946 |
/// // Run again with modified cost map (reset() is not called, |
... | ... |
@@ -947,7 +951,7 @@ |
947 | 951 |
/// // Run again from scratch using reset() |
948 | 952 |
/// // (the lower bounds will be set to zero on all arcs) |
949 | 953 |
/// ns.reset(); |
950 |
/// ns. |
|
954 |
/// ns.upperMap(capacity).costMap(cost) |
|
951 | 955 |
/// .supplyMap(sup).run(); |
952 | 956 |
/// \endcode |
953 | 957 |
/// |
... | ... |
@@ -962,7 +966,7 @@ |
962 | 966 |
_pcost = NULL; |
963 | 967 |
_psupply = NULL; |
964 | 968 |
_pstsup = false; |
965 |
|
|
969 |
_stype = GEQ; |
|
966 | 970 |
if (_local_flow) delete _flow_map; |
967 | 971 |
if (_local_potential) delete _potential_map; |
968 | 972 |
_flow_map = NULL; |
... | ... |
@@ -985,7 +989,7 @@ |
985 | 989 |
/// \brief Return the total cost of the found flow. |
986 | 990 |
/// |
987 | 991 |
/// This function returns the total cost of the found flow. |
988 |
/// |
|
992 |
/// Its complexity is O(e). |
|
989 | 993 |
/// |
990 | 994 |
/// \note The return type of the function can be specified as a |
991 | 995 |
/// template parameter. For example, |
... | ... |
@@ -997,9 +1001,9 @@ |
997 | 1001 |
/// function. |
998 | 1002 |
/// |
999 | 1003 |
/// \pre \ref run() must be called before using this function. |
1000 |
template <typename Num> |
|
1001 |
Num totalCost() const { |
|
1002 |
|
|
1004 |
template <typename Value> |
|
1005 |
Value totalCost() const { |
|
1006 |
Value c = 0; |
|
1003 | 1007 |
if (_pcost) { |
1004 | 1008 |
for (ArcIt e(_graph); e != INVALID; ++e) |
1005 | 1009 |
c += (*_flow_map)[e] * (*_pcost)[e]; |
... | ... |
@@ -1050,7 +1054,7 @@ |
1050 | 1054 |
/// |
1051 | 1055 |
/// This function returns a const reference to a node map storing |
1052 | 1056 |
/// the found potentials, which form the dual solution of the |
1053 |
/// \ref min_cost_flow "minimum cost flow" |
|
1057 |
/// \ref min_cost_flow "minimum cost flow problem". |
|
1054 | 1058 |
/// |
1055 | 1059 |
/// \pre \ref run() must be called before using this function. |
1056 | 1060 |
const PotentialMap& potentialMap() const { |
... | ... |
@@ -1101,7 +1105,7 @@ |
1101 | 1105 |
|
1102 | 1106 |
// Initialize node related data |
1103 | 1107 |
bool valid_supply = true; |
1104 |
|
|
1108 |
_sum_supply = 0; |
|
1105 | 1109 |
if (!_pstsup && !_psupply) { |
1106 | 1110 |
_pstsup = true; |
1107 | 1111 |
_psource = _ptarget = NodeIt(_graph); |
... | ... |
@@ -1112,10 +1116,10 @@ |
1112 | 1116 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
1113 | 1117 |
_node_id[n] = i; |
1114 | 1118 |
_supply[i] = (*_psupply)[n]; |
1115 |
|
|
1119 |
_sum_supply += _supply[i]; |
|
1116 | 1120 |
} |
1117 |
valid_supply = (_ptype == GEQ && sum_supply <= 0) || |
|
1118 |
(_ptype == LEQ && sum_supply >= 0); |
|
1121 |
valid_supply = (_stype == GEQ && _sum_supply <= 0) || |
|
1122 |
(_stype == LEQ && _sum_supply >= 0); |
|
1119 | 1123 |
} else { |
1120 | 1124 |
int i = 0; |
1121 | 1125 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
... | ... |
@@ -1127,92 +1131,18 @@ |
1127 | 1131 |
} |
1128 | 1132 |
if (!valid_supply) return false; |
1129 | 1133 |
|
1130 |
// Infinite capacity value |
|
1131 |
Flow inf_cap = |
|
1132 |
std::numeric_limits<Flow>::has_infinity ? |
|
1133 |
std::numeric_limits<Flow>::infinity() : |
|
1134 |
std::numeric_limits<Flow>::max(); |
|
1135 |
|
|
1136 | 1134 |
// Initialize artifical cost |
1137 |
Cost |
|
1135 |
Cost ART_COST; |
|
1138 | 1136 |
if (std::numeric_limits<Cost>::is_exact) { |
1139 |
|
|
1137 |
ART_COST = std::numeric_limits<Cost>::max() / 4 + 1; |
|
1140 | 1138 |
} else { |
1141 |
|
|
1139 |
ART_COST = std::numeric_limits<Cost>::min(); |
|
1142 | 1140 |
for (int i = 0; i != _arc_num; ++i) { |
1143 |
if (_cost[i] > |
|
1141 |
if (_cost[i] > ART_COST) ART_COST = _cost[i]; |
|
1144 | 1142 |
} |
1145 |
|
|
1143 |
ART_COST = (ART_COST + 1) * _node_num; |
|
1146 | 1144 |
} |
1147 | 1145 |
|
1148 |
// Run Circulation to check if a feasible solution exists |
|
1149 |
typedef ConstMap<Arc, Flow> ConstArcMap; |
|
1150 |
ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap); |
|
1151 |
FlowNodeMap *csup = NULL; |
|
1152 |
bool local_csup = false; |
|
1153 |
if (_psupply) { |
|
1154 |
csup = _psupply; |
|
1155 |
} else { |
|
1156 |
csup = new FlowNodeMap(_graph, 0); |
|
1157 |
(*csup)[_psource] = _pstflow; |
|
1158 |
(*csup)[_ptarget] = -_pstflow; |
|
1159 |
local_csup = true; |
|
1160 |
} |
|
1161 |
bool circ_result = false; |
|
1162 |
if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) { |
|
1163 |
// GEQ problem type |
|
1164 |
if (_plower) { |
|
1165 |
if (_pupper) { |
|
1166 |
Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap> |
|
1167 |
circ(_graph, *_plower, *_pupper, *csup); |
|
1168 |
circ_result = circ.run(); |
|
1169 |
} else { |
|
1170 |
Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap> |
|
1171 |
circ(_graph, *_plower, inf_arc_map, *csup); |
|
1172 |
circ_result = circ.run(); |
|
1173 |
} |
|
1174 |
} else { |
|
1175 |
if (_pupper) { |
|
1176 |
Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap> |
|
1177 |
circ(_graph, zero_arc_map, *_pupper, *csup); |
|
1178 |
circ_result = circ.run(); |
|
1179 |
} else { |
|
1180 |
Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap> |
|
1181 |
circ(_graph, zero_arc_map, inf_arc_map, *csup); |
|
1182 |
circ_result = circ.run(); |
|
1183 |
} |
|
1184 |
} |
|
1185 |
} else { |
|
1186 |
// LEQ problem type |
|
1187 |
typedef ReverseDigraph<const GR> RevGraph; |
|
1188 |
typedef NegMap<FlowNodeMap> NegNodeMap; |
|
1189 |
RevGraph rgraph(_graph); |
|
1190 |
NegNodeMap neg_csup(*csup); |
|
1191 |
if (_plower) { |
|
1192 |
if (_pupper) { |
|
1193 |
Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap> |
|
1194 |
circ(rgraph, *_plower, *_pupper, neg_csup); |
|
1195 |
circ_result = circ.run(); |
|
1196 |
} else { |
|
1197 |
Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap> |
|
1198 |
circ(rgraph, *_plower, inf_arc_map, neg_csup); |
|
1199 |
circ_result = circ.run(); |
|
1200 |
} |
|
1201 |
} else { |
|
1202 |
if (_pupper) { |
|
1203 |
Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap> |
|
1204 |
circ(rgraph, zero_arc_map, *_pupper, neg_csup); |
|
1205 |
circ_result = circ.run(); |
|
1206 |
} else { |
|
1207 |
Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap> |
|
1208 |
circ(rgraph, zero_arc_map, inf_arc_map, neg_csup); |
|
1209 |
circ_result = circ.run(); |
|
1210 |
} |
|
1211 |
} |
|
1212 |
} |
|
1213 |
if (local_csup) delete csup; |
|
1214 |
if (!circ_result) return false; |
|
1215 |
|
|
1216 | 1146 |
// Set data for the artificial root node |
1217 | 1147 |
_root = _node_num; |
1218 | 1148 |
_parent[_root] = -1; |
... | ... |
@@ -1221,11 +1151,11 @@ |
1221 | 1151 |
_rev_thread[0] = _root; |
1222 | 1152 |
_succ_num[_root] = all_node_num; |
1223 | 1153 |
_last_succ[_root] = _root - 1; |
1224 |
_supply[_root] = -sum_supply; |
|
1225 |
if (sum_supply < 0) { |
|
1226 |
|
|
1154 |
_supply[_root] = -_sum_supply; |
|
1155 |
if (_sum_supply < 0) { |
|
1156 |
_pi[_root] = -ART_COST; |
|
1227 | 1157 |
} else { |
1228 |
_pi[_root] = |
|
1158 |
_pi[_root] = ART_COST; |
|
1229 | 1159 |
} |
1230 | 1160 |
|
1231 | 1161 |
// Store the arcs in a mixed order |
... | ... |
@@ -1260,7 +1190,7 @@ |
1260 | 1190 |
_cap[i] = (*_pupper)[_arc_ref[i]]; |
1261 | 1191 |
} else { |
1262 | 1192 |
for (int i = 0; i != _arc_num; ++i) |
1263 |
_cap[i] = |
|
1193 |
_cap[i] = INF; |
|
1264 | 1194 |
} |
1265 | 1195 |
if (_pcost) { |
1266 | 1196 |
for (int i = 0; i != _arc_num; ++i) |
... | ... |
@@ -1275,8 +1205,17 @@ |
1275 | 1205 |
if (_plower) { |
1276 | 1206 |
for (int i = 0; i != _arc_num; ++i) { |
1277 | 1207 |
Flow c = (*_plower)[_arc_ref[i]]; |
1278 |
if (c |
|
1208 |
if (c > 0) { |
|
1209 |
if (_cap[i] < INF) _cap[i] -= c; |
|
1210 |
_supply[_source[i]] -= c; |
|
1211 |
_supply[_target[i]] += c; |
|
1212 |
} |
|
1213 |
else if (c < 0) { |
|
1214 |
if (_cap[i] < INF + c) { |
|
1279 | 1215 |
_cap[i] -= c; |
1216 |
} else { |
|
1217 |
_cap[i] = INF; |
|
1218 |
} |
|
1280 | 1219 |
_supply[_source[i]] -= c; |
1281 | 1220 |
_supply[_target[i]] += c; |
1282 | 1221 |
} |
... | ... |
@@ -1291,17 +1230,17 @@ |
1291 | 1230 |
_last_succ[u] = u; |
1292 | 1231 |
_parent[u] = _root; |
1293 | 1232 |
_pred[u] = e; |
1294 |
_cost[e] = art_cost; |
|
1295 |
_cap[e] = inf_cap; |
|
1233 |
_cost[e] = ART_COST; |
|
1234 |
_cap[e] = INF; |
|
1296 | 1235 |
_state[e] = STATE_TREE; |
1297 |
if (_supply[u] > 0 || (_supply[u] == 0 && |
|
1236 |
if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) { |
|
1298 | 1237 |
_flow[e] = _supply[u]; |
1299 | 1238 |
_forward[u] = true; |
1300 |
_pi[u] = - |
|
1239 |
_pi[u] = -ART_COST + _pi[_root]; |
|
1301 | 1240 |
} else { |
1302 | 1241 |
_flow[e] = -_supply[u]; |
1303 | 1242 |
_forward[u] = false; |
1304 |
_pi[u] = |
|
1243 |
_pi[u] = ART_COST + _pi[_root]; |
|
1305 | 1244 |
} |
1306 | 1245 |
} |
1307 | 1246 |
|
... | ... |
@@ -1342,7 +1281,8 @@ |
1342 | 1281 |
// Search the cycle along the path form the first node to the root |
1343 | 1282 |
for (int u = first; u != join; u = _parent[u]) { |
1344 | 1283 |
e = _pred[u]; |
1345 |
d = _forward[u] ? |
|
1284 |
d = _forward[u] ? |
|
1285 |
_flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]); |
|
1346 | 1286 |
if (d < delta) { |
1347 | 1287 |
delta = d; |
1348 | 1288 |
u_out = u; |
... | ... |
@@ -1352,7 +1292,8 @@ |
1352 | 1292 |
// Search the cycle along the path form the second node to the root |
1353 | 1293 |
for (int u = second; u != join; u = _parent[u]) { |
1354 | 1294 |
e = _pred[u]; |
1355 |
d = _forward[u] ? |
|
1295 |
d = _forward[u] ? |
|
1296 |
(_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e]; |
|
1356 | 1297 |
if (d <= delta) { |
1357 | 1298 |
delta = d; |
1358 | 1299 |
u_out = u; |
... | ... |
@@ -1526,7 +1467,7 @@ |
1526 | 1467 |
} |
1527 | 1468 |
|
1528 | 1469 |
// Execute the algorithm |
1529 |
|
|
1470 |
ProblemType start(PivotRule pivot_rule) { |
|
1530 | 1471 |
// Select the pivot rule implementation |
1531 | 1472 |
switch (pivot_rule) { |
1532 | 1473 |
case FIRST_ELIGIBLE: |
... | ... |
@@ -1540,17 +1481,18 @@ |
1540 | 1481 |
case ALTERING_LIST: |
1541 | 1482 |
return start<AlteringListPivotRule>(); |
1542 | 1483 |
} |
1543 |
return |
|
1484 |
return INFEASIBLE; // avoid warning |
|
1544 | 1485 |
} |
1545 | 1486 |
|
1546 | 1487 |
template <typename PivotRuleImpl> |
1547 |
|
|
1488 |
ProblemType start() { |
|
1548 | 1489 |
PivotRuleImpl pivot(*this); |
1549 | 1490 |
|
1550 | 1491 |
// Execute the Network Simplex algorithm |
1551 | 1492 |
while (pivot.findEnteringArc()) { |
1552 | 1493 |
findJoinNode(); |
1553 | 1494 |
bool change = findLeavingArc(); |
1495 |
if (delta >= INF) return UNBOUNDED; |
|
1554 | 1496 |
changeFlow(change); |
1555 | 1497 |
if (change) { |
1556 | 1498 |
updateTreeStructure(); |
... | ... |
@@ -1558,6 +1500,23 @@ |
1558 | 1500 |
} |
1559 | 1501 |
} |
1560 | 1502 |
|
1503 |
// Check feasibility |
|
1504 |
if (_sum_supply < 0) { |
|
1505 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1506 |
if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE; |
|
1507 |
} |
|
1508 |
} |
|
1509 |
else if (_sum_supply > 0) { |
|
1510 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1511 |
if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE; |
|
1512 |
} |
|
1513 |
} |
|
1514 |
else { |
|
1515 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1516 |
if (_flow[e] != 0) return INFEASIBLE; |
|
1517 |
} |
|
1518 |
} |
|
1519 |
|
|
1561 | 1520 |
// Copy flow values to _flow_map |
1562 | 1521 |
if (_plower) { |
1563 | 1522 |
for (int i = 0; i != _arc_num; ++i) { |
... | ... |
@@ -1574,7 +1533,7 @@ |
1574 | 1533 |
_potential_map->set(n, _pi[_node_id[n]]); |
1575 | 1534 |
} |
1576 | 1535 |
|
1577 |
return |
|
1536 |
return OPTIMAL; |
|
1578 | 1537 |
} |
1579 | 1538 |
|
1580 | 1539 |
}; //class NetworkSimplex |
... | ... |
@@ -18,6 +18,7 @@ |
18 | 18 |
|
19 | 19 |
#include <iostream> |
20 | 20 |
#include <fstream> |
21 |
#include <limits> |
|
21 | 22 |
|
22 | 23 |
#include <lemon/list_graph.h> |
23 | 24 |
#include <lemon/lgf_reader.h> |
... | ... |
@@ -33,50 +34,50 @@ |
33 | 34 |
|
34 | 35 |
char test_lgf[] = |
35 | 36 |
"@nodes\n" |
36 |
"label sup1 sup2 sup3 sup4 sup5\n" |
|
37 |
" 1 20 27 0 20 30\n" |
|
38 |
" 2 -4 0 0 -8 -3\n" |
|
39 |
" 3 0 0 0 0 0\n" |
|
40 |
" 4 0 0 0 0 0\n" |
|
41 |
" 5 9 0 0 6 11\n" |
|
42 |
" 6 -6 0 0 -5 -6\n" |
|
43 |
" 7 0 0 0 0 0\n" |
|
44 |
" 8 0 0 0 0 3\n" |
|
45 |
" 9 3 0 0 0 0\n" |
|
46 |
" 10 -2 0 0 -7 -2\n" |
|
47 |
" 11 0 0 0 -10 0\n" |
|
48 |
" |
|
37 |
"label sup1 sup2 sup3 sup4 sup5 sup6\n" |
|
38 |
" 1 20 27 0 30 20 30\n" |
|
39 |
" 2 -4 0 0 0 -8 -3\n" |
|
40 |
" 3 0 0 0 0 0 0\n" |
|
41 |
" 4 0 0 0 0 0 0\n" |
|
42 |
" 5 9 0 0 0 6 11\n" |
|
43 |
" 6 -6 0 0 0 -5 -6\n" |
|
44 |
" 7 0 0 0 0 0 0\n" |
|
45 |
" 8 0 0 0 0 0 3\n" |
|
46 |
" 9 3 0 0 0 0 0\n" |
|
47 |
" 10 -2 0 0 0 -7 -2\n" |
|
48 |
" 11 0 0 0 0 -10 0\n" |
|
49 |
" 12 -20 -27 0 -30 -30 -20\n" |
|
49 | 50 |
"\n" |
50 | 51 |
"@arcs\n" |
51 |
" cost cap low1 low2\n" |
|
52 |
" 1 2 70 11 0 8\n" |
|
53 |
" 1 3 150 3 0 1\n" |
|
54 |
" 1 4 80 15 0 2\n" |
|
55 |
" 2 8 80 12 0 0\n" |
|
56 |
" 3 5 140 5 0 3\n" |
|
57 |
" 4 6 60 10 0 1\n" |
|
58 |
" 4 7 80 2 0 0\n" |
|
59 |
" 4 8 110 3 0 0\n" |
|
60 |
" 5 7 60 14 0 0\n" |
|
61 |
" 5 11 120 12 0 0\n" |
|
62 |
" 6 3 0 3 0 0\n" |
|
63 |
" 6 9 140 4 0 0\n" |
|
64 |
" 6 10 90 8 0 0\n" |
|
65 |
" 7 1 30 5 0 0\n" |
|
66 |
" 8 12 60 16 0 4\n" |
|
67 |
" 9 12 50 6 0 0\n" |
|
68 |
"10 12 70 13 0 5\n" |
|
69 |
"10 2 100 7 0 0\n" |
|
70 |
"10 7 60 10 0 0\n" |
|
71 |
"11 10 20 14 0 6\n" |
|
72 |
"12 11 30 10 0 0\n" |
|
52 |
" cost cap low1 low2 low3\n" |
|
53 |
" 1 2 70 11 0 8 8\n" |
|
54 |
" 1 3 150 3 0 1 0\n" |
|
55 |
" 1 4 80 15 0 2 2\n" |
|
56 |
" 2 8 80 12 0 0 0\n" |
|
57 |
" 3 5 140 5 0 3 1\n" |
|
58 |
" 4 6 60 10 0 1 0\n" |
|
59 |
" 4 7 80 2 0 0 0\n" |
|
60 |
" 4 8 110 3 0 0 0\n" |
|
61 |
" 5 7 60 14 0 0 0\n" |
|
62 |
" 5 11 120 12 0 0 0\n" |
|
63 |
" 6 3 0 3 0 0 0\n" |
|
64 |
" 6 9 140 4 0 0 0\n" |
|
65 |
" 6 10 90 8 0 0 0\n" |
|
66 |
" 7 1 30 5 0 0 -5\n" |
|
67 |
" 8 12 60 16 0 4 3\n" |
|
68 |
" 9 12 50 6 0 0 0\n" |
|
69 |
"10 12 70 13 0 5 2\n" |
|
70 |
"10 2 100 7 0 0 0\n" |
|
71 |
"10 7 60 10 0 0 -3\n" |
|
72 |
"11 10 20 14 0 6 -20\n" |
|
73 |
"12 11 30 10 0 0 -10\n" |
|
73 | 74 |
"\n" |
74 | 75 |
"@attributes\n" |
75 | 76 |
"source 1\n" |
76 | 77 |
"target 12\n"; |
77 | 78 |
|
78 | 79 |
|
79 |
enum |
|
80 |
enum SupplyType { |
|
80 | 81 |
EQ, |
81 | 82 |
GEQ, |
82 | 83 |
LEQ |
... | ... |
@@ -98,8 +99,6 @@ |
98 | 99 |
b = mcf.reset() |
99 | 100 |
.lowerMap(lower) |
100 | 101 |
.upperMap(upper) |
101 |
.capacityMap(upper) |
|
102 |
.boundMaps(lower, upper) |
|
103 | 102 |
.costMap(cost) |
104 | 103 |
.supplyMap(sup) |
105 | 104 |
.stSupply(n, n, k) |
... | ... |
@@ -112,10 +111,12 @@ |
112 | 111 |
const typename MCF::FlowMap &fm = const_mcf.flowMap(); |
113 | 112 |
const typename MCF::PotentialMap &pm = const_mcf.potentialMap(); |
114 | 113 |
|
115 |
|
|
114 |
c = const_mcf.totalCost(); |
|
116 | 115 |
double x = const_mcf.template totalCost<double>(); |
117 | 116 |
v = const_mcf.flow(a); |
118 |
|
|
117 |
c = const_mcf.potential(n); |
|
118 |
|
|
119 |
v = const_mcf.INF; |
|
119 | 120 |
|
120 | 121 |
ignore_unused_variable_warning(fm); |
121 | 122 |
ignore_unused_variable_warning(pm); |
... | ... |
@@ -137,6 +138,7 @@ |
137 | 138 |
const Arc &a; |
138 | 139 |
const Flow &k; |
139 | 140 |
Flow v; |
141 |
Cost c; |
|
140 | 142 |
bool b; |
141 | 143 |
|
142 | 144 |
typename MCF::FlowMap &flow; |
... | ... |
@@ -151,7 +153,7 @@ |
151 | 153 |
typename SM, typename FM > |
152 | 154 |
bool checkFlow( const GR& gr, const LM& lower, const UM& upper, |
153 | 155 |
const SM& supply, const FM& flow, |
154 |
|
|
156 |
SupplyType type = EQ ) |
|
155 | 157 |
{ |
156 | 158 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
157 | 159 |
|
... | ... |
@@ -208,16 +210,17 @@ |
208 | 210 |
// Run a minimum cost flow algorithm and check the results |
209 | 211 |
template < typename MCF, typename GR, |
210 | 212 |
typename LM, typename UM, |
211 |
typename CM, typename SM > |
|
212 |
void checkMcf( const MCF& mcf, bool mcf_result, |
|
213 |
typename CM, typename SM, |
|
214 |
typename PT > |
|
215 |
void checkMcf( const MCF& mcf, PT mcf_result, |
|
213 | 216 |
const GR& gr, const LM& lower, const UM& upper, |
214 | 217 |
const CM& cost, const SM& supply, |
215 |
bool |
|
218 |
PT result, bool optimal, typename CM::Value total, |
|
216 | 219 |
const std::string &test_id = "", |
217 |
|
|
220 |
SupplyType type = EQ ) |
|
218 | 221 |
{ |
219 | 222 |
check(mcf_result == result, "Wrong result " + test_id); |
220 |
if ( |
|
223 |
if (optimal) { |
|
221 | 224 |
check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type), |
222 | 225 |
"The flow is not feasible " + test_id); |
223 | 226 |
check(mcf.totalCost() == total, "The flow is not optimal " + test_id); |
... | ... |
@@ -244,8 +247,8 @@ |
244 | 247 |
|
245 | 248 |
// Read the test digraph |
246 | 249 |
Digraph gr; |
247 |
Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr); |
|
248 |
Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr); |
|
250 |
Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr); |
|
251 |
Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr); |
|
249 | 252 |
ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max()); |
250 | 253 |
Node v, w; |
251 | 254 |
|
... | ... |
@@ -255,15 +258,57 @@ |
255 | 258 |
.arcMap("cap", u) |
256 | 259 |
.arcMap("low1", l1) |
257 | 260 |
.arcMap("low2", l2) |
261 |
.arcMap("low3", l3) |
|
258 | 262 |
.nodeMap("sup1", s1) |
259 | 263 |
.nodeMap("sup2", s2) |
260 | 264 |
.nodeMap("sup3", s3) |
261 | 265 |
.nodeMap("sup4", s4) |
262 | 266 |
.nodeMap("sup5", s5) |
267 |
.nodeMap("sup6", s6) |
|
263 | 268 |
.node("source", v) |
264 | 269 |
.node("target", w) |
265 | 270 |
.run(); |
266 | 271 |
|
272 |
// Build a test digraph for testing negative costs |
|
273 |
Digraph ngr; |
|
274 |
Node n1 = ngr.addNode(); |
|
275 |
Node n2 = ngr.addNode(); |
|
276 |
Node n3 = ngr.addNode(); |
|
277 |
Node n4 = ngr.addNode(); |
|
278 |
Node n5 = ngr.addNode(); |
|
279 |
Node n6 = ngr.addNode(); |
|
280 |
Node n7 = ngr.addNode(); |
|
281 |
|
|
282 |
Arc a1 = ngr.addArc(n1, n2); |
|
283 |
Arc a2 = ngr.addArc(n1, n3); |
|
284 |
Arc a3 = ngr.addArc(n2, n4); |
|
285 |
Arc a4 = ngr.addArc(n3, n4); |
|
286 |
Arc a5 = ngr.addArc(n3, n2); |
|
287 |
Arc a6 = ngr.addArc(n5, n3); |
|
288 |
Arc a7 = ngr.addArc(n5, n6); |
|
289 |
Arc a8 = ngr.addArc(n6, n7); |
|
290 |
Arc a9 = ngr.addArc(n7, n5); |
|
291 |
|
|
292 |
Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0); |
|
293 |
ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000); |
|
294 |
Digraph::NodeMap<int> ns(ngr, 0); |
|
295 |
|
|
296 |
nl2[a7] = 1000; |
|
297 |
nl2[a8] = -1000; |
|
298 |
|
|
299 |
ns[n1] = 100; |
|
300 |
ns[n4] = -100; |
|
301 |
|
|
302 |
nc[a1] = 100; |
|
303 |
nc[a2] = 30; |
|
304 |
nc[a3] = 20; |
|
305 |
nc[a4] = 80; |
|
306 |
nc[a5] = 50; |
|
307 |
nc[a6] = 10; |
|
308 |
nc[a7] = 80; |
|
309 |
nc[a8] = 30; |
|
310 |
nc[a9] = -120; |
|
311 |
|
|
267 | 312 |
// A. Test NetworkSimplex with the default pivot rule |
268 | 313 |
{ |
269 | 314 |
NetworkSimplex<Digraph> mcf(gr); |
... | ... |
@@ -271,63 +316,77 @@ |
271 | 316 |
// Check the equality form |
272 | 317 |
mcf.upperMap(u).costMap(c); |
273 | 318 |
checkMcf(mcf, mcf.supplyMap(s1).run(), |
274 |
gr, l1, u, c, s1, true, 5240, "#A1"); |
|
319 |
gr, l1, u, c, s1, mcf.OPTIMAL, true, 5240, "#A1"); |
|
275 | 320 |
checkMcf(mcf, mcf.stSupply(v, w, 27).run(), |
276 |
gr, l1, u, c, s2, true, 7620, "#A2"); |
|
321 |
gr, l1, u, c, s2, mcf.OPTIMAL, true, 7620, "#A2"); |
|
277 | 322 |
mcf.lowerMap(l2); |
278 | 323 |
checkMcf(mcf, mcf.supplyMap(s1).run(), |
279 |
gr, l2, u, c, s1, true, 5970, "#A3"); |
|
324 |
gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#A3"); |
|
280 | 325 |
checkMcf(mcf, mcf.stSupply(v, w, 27).run(), |
281 |
gr, l2, u, c, s2, true, 8010, "#A4"); |
|
326 |
gr, l2, u, c, s2, mcf.OPTIMAL, true, 8010, "#A4"); |
|
282 | 327 |
mcf.reset(); |
283 | 328 |
checkMcf(mcf, mcf.supplyMap(s1).run(), |
284 |
gr, l1, cu, cc, s1, true, 74, "#A5"); |
|
329 |
gr, l1, cu, cc, s1, mcf.OPTIMAL, true, 74, "#A5"); |
|
285 | 330 |
checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(), |
286 |
gr, l2, cu, cc, s2, true, 94, "#A6"); |
|
331 |
gr, l2, cu, cc, s2, mcf.OPTIMAL, true, 94, "#A6"); |
|
287 | 332 |
mcf.reset(); |
288 | 333 |
checkMcf(mcf, mcf.run(), |
289 |
gr, l1, cu, cc, s3, true, 0, "#A7"); |
|
290 |
checkMcf(mcf, mcf.boundMaps(l2, u).run(), |
|
291 |
gr, |
|
334 |
gr, l1, cu, cc, s3, mcf.OPTIMAL, true, 0, "#A7"); |
|
335 |
checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(), |
|
336 |
gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8"); |
|
337 |
mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4); |
|
338 |
checkMcf(mcf, mcf.run(), |
|
339 |
gr, l3, u, c, s4, mcf.OPTIMAL, true, 6360, "#A9"); |
|
292 | 340 |
|
293 | 341 |
// Check the GEQ form |
294 |
mcf.reset().upperMap(u).costMap(c).supplyMap( |
|
342 |
mcf.reset().upperMap(u).costMap(c).supplyMap(s5); |
|
295 | 343 |
checkMcf(mcf, mcf.run(), |
296 |
gr, l1, u, c, s4, true, 3530, "#A9", GEQ); |
|
297 |
mcf.problemType(mcf.GEQ); |
|
344 |
gr, l1, u, c, s5, mcf.OPTIMAL, true, 3530, "#A10", GEQ); |
|
345 |
mcf.supplyType(mcf.GEQ); |
|
298 | 346 |
checkMcf(mcf, mcf.lowerMap(l2).run(), |
299 |
gr, l2, u, c, s4, true, 4540, "#A10", GEQ); |
|
300 |
mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5); |
|
347 |
gr, l2, u, c, s5, mcf.OPTIMAL, true, 4540, "#A11", GEQ); |
|
348 |
mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6); |
|
301 | 349 |
checkMcf(mcf, mcf.run(), |
302 |
gr, l2, u, c, |
|
350 |
gr, l2, u, c, s6, mcf.INFEASIBLE, false, 0, "#A12", GEQ); |
|
303 | 351 |
|
304 | 352 |
// Check the LEQ form |
305 |
mcf.reset().problemType(mcf.LEQ); |
|
306 |
mcf.upperMap(u).costMap(c).supplyMap(s5); |
|
353 |
mcf.reset().supplyType(mcf.LEQ); |
|
354 |
mcf.upperMap(u).costMap(c).supplyMap(s6); |
|
307 | 355 |
checkMcf(mcf, mcf.run(), |
308 |
gr, l1, u, c, |
|
356 |
gr, l1, u, c, s6, mcf.OPTIMAL, true, 5080, "#A13", LEQ); |
|
309 | 357 |
checkMcf(mcf, mcf.lowerMap(l2).run(), |
310 |
gr, l2, u, c, s5, true, 5930, "#A13", LEQ); |
|
311 |
mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4); |
|
358 |
gr, l2, u, c, s6, mcf.OPTIMAL, true, 5930, "#A14", LEQ); |
|
359 |
mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5); |
|
312 | 360 |
checkMcf(mcf, mcf.run(), |
313 |
gr, l2, u, c, |
|
361 |
gr, l2, u, c, s5, mcf.INFEASIBLE, false, 0, "#A15", LEQ); |
|
362 |
|
|
363 |
// Check negative costs |
|
364 |
NetworkSimplex<Digraph> nmcf(ngr); |
|
365 |
nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns); |
|
366 |
checkMcf(nmcf, nmcf.run(), |
|
367 |
ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A16"); |
|
368 |
checkMcf(nmcf, nmcf.upperMap(nu2).run(), |
|
369 |
ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true, -40000, "#A17"); |
|
370 |
nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns); |
|
371 |
checkMcf(nmcf, nmcf.run(), |
|
372 |
ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false, 0, "#A18"); |
|
314 | 373 |
} |
315 | 374 |
|
316 | 375 |
// B. Test NetworkSimplex with each pivot rule |
317 | 376 |
{ |
318 | 377 |
NetworkSimplex<Digraph> mcf(gr); |
319 |
mcf.supplyMap(s1).costMap(c). |
|
378 |
mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2); |
|
320 | 379 |
|
321 | 380 |
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE), |
322 |
gr, l2, u, c, s1, true, 5970, "#B1"); |
|
381 |
gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B1"); |
|
323 | 382 |
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE), |
324 |
gr, l2, u, c, s1, true, 5970, "#B2"); |
|
383 |
gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B2"); |
|
325 | 384 |
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH), |
326 |
gr, l2, u, c, s1, true, 5970, "#B3"); |
|
385 |
gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B3"); |
|
327 | 386 |
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST), |
328 |
gr, l2, u, c, s1, true, 5970, "#B4"); |
|
387 |
gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B4"); |
|
329 | 388 |
checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST), |
330 |
gr, l2, u, c, s1, true, 5970, "#B5"); |
|
389 |
gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B5"); |
|
331 | 390 |
} |
332 | 391 |
|
333 | 392 |
return 0; |
... | ... |
@@ -119,8 +119,8 @@ |
119 | 119 |
|
120 | 120 |
ti.restart(); |
121 | 121 |
NetworkSimplex<Digraph, Value> ns(g); |
122 |
ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup); |
|
123 |
if (sum_sup > 0) ns.problemType(ns.LEQ); |
|
122 |
ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup); |
|
123 |
if (sum_sup > 0) ns.supplyType(ns.LEQ); |
|
124 | 124 |
if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n'; |
125 | 125 |
ti.restart(); |
126 | 126 |
bool res = ns.run(); |
0 comments (0 inline)