0
3
1
1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*- |
|
2 |
* |
|
3 |
* This file is a part of LEMON, a generic C++ optimization library. |
|
4 |
* |
|
5 |
* Copyright (C) 2003-2009 |
|
6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
|
7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES). |
|
8 |
* |
|
9 |
* Permission to use, modify and distribute this software is granted |
|
10 |
* provided that this copyright notice appears in all copies. For |
|
11 |
* precise terms see the accompanying LICENSE file. |
|
12 |
* |
|
13 |
* This software is provided "AS IS" with no warranty of any kind, |
|
14 |
* express or implied, and with no claim as to its suitability for any |
|
15 |
* purpose. |
|
16 |
* |
|
17 |
*/ |
|
18 |
|
|
19 |
namespace lemon { |
|
20 |
|
|
21 |
/** |
|
22 |
\page min_cost_flow Minimum Cost Flow Problem |
|
23 |
|
|
24 |
\section mcf_def Definition (GEQ form) |
|
25 |
|
|
26 |
The \e minimum \e cost \e flow \e problem is to find a feasible flow of |
|
27 |
minimum total cost from a set of supply nodes to a set of demand nodes |
|
28 |
in a network with capacity constraints (lower and upper bounds) |
|
29 |
and arc costs. |
|
30 |
|
|
31 |
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$, |
|
32 |
\f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and |
|
33 |
upper bounds for the flow values on the arcs, for which |
|
34 |
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, |
|
35 |
\f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow |
|
36 |
on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the |
|
37 |
signed supply values of the nodes. |
|
38 |
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ |
|
39 |
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with |
|
40 |
\f$-sup(u)\f$ demand. |
|
41 |
A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution |
|
42 |
of the following optimization problem. |
|
43 |
|
|
44 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
45 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq |
|
46 |
sup(u) \quad \forall u\in V \f] |
|
47 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
48 |
|
|
49 |
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be |
|
50 |
zero or negative in order to have a feasible solution (since the sum |
|
51 |
of the expressions on the left-hand side of the inequalities is zero). |
|
52 |
It means that the total demand must be greater or equal to the total |
|
53 |
supply and all the supplies have to be carried out from the supply nodes, |
|
54 |
but there could be demands that are not satisfied. |
|
55 |
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand |
|
56 |
constraints have to be satisfied with equality, i.e. all demands |
|
57 |
have to be satisfied and all supplies have to be used. |
|
58 |
|
|
59 |
|
|
60 |
\section mcf_algs Algorithms |
|
61 |
|
|
62 |
LEMON contains several algorithms for solving this problem, for more |
|
63 |
information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms". |
|
64 |
|
|
65 |
A feasible solution for this problem can be found using \ref Circulation. |
|
66 |
|
|
67 |
|
|
68 |
\section mcf_dual Dual Solution |
|
69 |
|
|
70 |
The dual solution of the minimum cost flow problem is represented by |
|
71 |
node potentials \f$\pi: V\rightarrow\mathbf{R}\f$. |
|
72 |
An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal |
|
73 |
if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials |
|
74 |
the following \e complementary \e slackness optimality conditions hold. |
|
75 |
|
|
76 |
- For all \f$uv\in A\f$ arcs: |
|
77 |
- if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; |
|
78 |
- if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; |
|
79 |
- if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. |
|
80 |
- For all \f$u\in V\f$ nodes: |
|
81 |
- \f$\pi(u)<=0\f$; |
|
82 |
- if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$, |
|
83 |
then \f$\pi(u)=0\f$. |
|
84 |
|
|
85 |
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc |
|
86 |
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e. |
|
87 |
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f] |
|
88 |
|
|
89 |
All algorithms provide dual solution (node potentials), as well, |
|
90 |
if an optimal flow is found. |
|
91 |
|
|
92 |
|
|
93 |
\section mcf_eq Equality Form |
|
94 |
|
|
95 |
The above \ref mcf_def "definition" is actually more general than the |
|
96 |
usual formulation of the minimum cost flow problem, in which strict |
|
97 |
equalities are required in the supply/demand contraints. |
|
98 |
|
|
99 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
100 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) = |
|
101 |
sup(u) \quad \forall u\in V \f] |
|
102 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
103 |
|
|
104 |
However if the sum of the supply values is zero, then these two problems |
|
105 |
are equivalent. |
|
106 |
The \ref min_cost_flow_algs "algorithms" in LEMON support the general |
|
107 |
form, so if you need the equality form, you have to ensure this additional |
|
108 |
contraint manually. |
|
109 |
|
|
110 |
|
|
111 |
\section mcf_leq Opposite Inequalites (LEQ Form) |
|
112 |
|
|
113 |
Another possible definition of the minimum cost flow problem is |
|
114 |
when there are <em>"less or equal"</em> (LEQ) supply/demand constraints, |
|
115 |
instead of the <em>"greater or equal"</em> (GEQ) constraints. |
|
116 |
|
|
117 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
118 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq |
|
119 |
sup(u) \quad \forall u\in V \f] |
|
120 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
121 |
|
|
122 |
It means that the total demand must be less or equal to the |
|
123 |
total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or |
|
124 |
positive) and all the demands have to be satisfied, but there |
|
125 |
could be supplies that are not carried out from the supply |
|
126 |
nodes. |
|
127 |
The equality form is also a special case of this form, of course. |
|
128 |
|
|
129 |
You could easily transform this case to the \ref mcf_def "GEQ form" |
|
130 |
of the problem by reversing the direction of the arcs and taking the |
|
131 |
negative of the supply values (e.g. using \ref ReverseDigraph and |
|
132 |
\ref NegMap adaptors). |
|
133 |
However \ref NetworkSimplex algorithm also supports this form directly |
|
134 |
for the sake of convenience. |
|
135 |
|
|
136 |
Note that the optimality conditions for this supply constraint type are |
|
137 |
slightly differ from the conditions that are discussed for the GEQ form, |
|
138 |
namely the potentials have to be non-negative instead of non-positive. |
|
139 |
An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem |
|
140 |
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ |
|
141 |
node potentials the following conditions hold. |
|
142 |
|
|
143 |
- For all \f$uv\in A\f$ arcs: |
|
144 |
- if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; |
|
145 |
- if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; |
|
146 |
- if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. |
|
147 |
- For all \f$u\in V\f$ nodes: |
|
148 |
- \f$\pi(u)>=0\f$; |
|
149 |
- if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$, |
|
150 |
then \f$\pi(u)=0\f$. |
|
151 |
|
|
152 |
*/ |
|
153 |
} |
1 | 1 |
EXTRA_DIST += \ |
2 | 2 |
doc/Doxyfile.in \ |
3 | 3 |
doc/DoxygenLayout.xml \ |
4 | 4 |
doc/coding_style.dox \ |
5 | 5 |
doc/dirs.dox \ |
6 | 6 |
doc/groups.dox \ |
7 | 7 |
doc/lgf.dox \ |
8 | 8 |
doc/license.dox \ |
9 | 9 |
doc/mainpage.dox \ |
10 | 10 |
doc/migration.dox \ |
11 |
doc/min_cost_flow.dox \ |
|
11 | 12 |
doc/named-param.dox \ |
12 | 13 |
doc/namespaces.dox \ |
13 | 14 |
doc/html \ |
14 | 15 |
doc/CMakeLists.txt |
15 | 16 |
|
16 | 17 |
DOC_EPS_IMAGES18 = \ |
17 | 18 |
grid_graph.eps \ |
18 | 19 |
nodeshape_0.eps \ |
19 | 20 |
nodeshape_1.eps \ |
20 | 21 |
nodeshape_2.eps \ |
21 | 22 |
nodeshape_3.eps \ |
22 | 23 |
nodeshape_4.eps |
23 | 24 |
|
24 | 25 |
DOC_EPS_IMAGES27 = \ |
25 | 26 |
bipartite_matching.eps \ |
26 | 27 |
bipartite_partitions.eps \ |
27 | 28 |
connected_components.eps \ |
28 | 29 |
edge_biconnected_components.eps \ |
29 | 30 |
node_biconnected_components.eps \ |
30 | 31 |
strongly_connected_components.eps |
31 | 32 |
|
32 | 33 |
DOC_EPS_IMAGES = \ |
33 | 34 |
$(DOC_EPS_IMAGES18) \ |
34 | 35 |
$(DOC_EPS_IMAGES27) |
35 | 36 |
|
36 | 37 |
DOC_PNG_IMAGES = \ |
37 | 38 |
$(DOC_EPS_IMAGES:%.eps=doc/gen-images/%.png) |
38 | 39 |
|
39 | 40 |
EXTRA_DIST += $(DOC_EPS_IMAGES:%=doc/images/%) |
40 | 41 |
|
41 | 42 |
doc/html: |
42 | 43 |
$(MAKE) $(AM_MAKEFLAGS) html |
43 | 44 |
|
44 | 45 |
GS_COMMAND=gs -dNOPAUSE -dBATCH -q -dEPSCrop -dTextAlphaBits=4 -dGraphicsAlphaBits=4 |
45 | 46 |
|
46 | 47 |
$(DOC_EPS_IMAGES18:%.eps=doc/gen-images/%.png): doc/gen-images/%.png: doc/images/%.eps |
47 | 48 |
-mkdir doc/gen-images |
48 | 49 |
if test ${gs_found} = yes; then \ |
49 | 50 |
$(GS_COMMAND) -sDEVICE=pngalpha -r18 -sOutputFile=$@ $<; \ |
50 | 51 |
else \ |
51 | 52 |
echo; \ |
52 | 53 |
echo "Ghostscript not found."; \ |
53 | 54 |
echo; \ |
54 | 55 |
exit 1; \ |
55 | 56 |
fi |
56 | 57 |
|
57 | 58 |
$(DOC_EPS_IMAGES27:%.eps=doc/gen-images/%.png): doc/gen-images/%.png: doc/images/%.eps |
58 | 59 |
-mkdir doc/gen-images |
59 | 60 |
if test ${gs_found} = yes; then \ |
60 | 61 |
$(GS_COMMAND) -sDEVICE=pngalpha -r27 -sOutputFile=$@ $<; \ |
61 | 62 |
else \ |
62 | 63 |
echo; \ |
63 | 64 |
echo "Ghostscript not found."; \ |
64 | 65 |
echo; \ |
65 | 66 |
exit 1; \ |
66 | 67 |
fi |
67 | 68 |
|
68 | 69 |
html-local: $(DOC_PNG_IMAGES) |
69 | 70 |
if test ${doxygen_found} = yes; then \ |
70 | 71 |
cd doc; \ |
71 | 72 |
doxygen Doxyfile; \ |
72 | 73 |
cd ..; \ |
73 | 74 |
else \ |
74 | 75 |
echo; \ |
75 | 76 |
echo "Doxygen not found."; \ |
76 | 77 |
echo; \ |
77 | 78 |
exit 1; \ |
78 | 79 |
fi |
79 | 80 |
|
80 | 81 |
clean-local: |
81 | 82 |
-rm -rf doc/html |
82 | 83 |
-rm -f doc/doxygen.log |
83 | 84 |
-rm -f $(DOC_PNG_IMAGES) |
84 | 85 |
-rm -rf doc/gen-images |
85 | 86 |
|
86 | 87 |
update-external-tags: |
87 | 88 |
wget -O doc/libstdc++.tag.tmp http://gcc.gnu.org/onlinedocs/libstdc++/latest-doxygen/libstdc++.tag && \ |
88 | 89 |
mv doc/libstdc++.tag.tmp doc/libstdc++.tag || \ |
89 | 90 |
rm doc/libstdc++.tag.tmp |
90 | 91 |
|
91 | 92 |
install-html-local: doc/html |
92 | 93 |
@$(NORMAL_INSTALL) |
93 | 94 |
$(mkinstalldirs) $(DESTDIR)$(htmldir)/docs |
94 | 95 |
for p in doc/html/*.{html,css,png,map,gif,tag} ; do \ |
95 | 96 |
f="`echo $$p | sed -e 's|^.*/||'`"; \ |
96 | 97 |
echo " $(INSTALL_DATA) $$p $(DESTDIR)$(htmldir)/docs/$$f"; \ |
97 | 98 |
$(INSTALL_DATA) $$p $(DESTDIR)$(htmldir)/docs/$$f; \ |
98 | 99 |
done |
99 | 100 |
|
100 | 101 |
uninstall-local: |
101 | 102 |
@$(NORMAL_UNINSTALL) |
102 | 103 |
for p in doc/html/*.{html,css,png,map,gif,tag} ; do \ |
103 | 104 |
f="`echo $$p | sed -e 's|^.*/||'`"; \ |
104 | 105 |
echo " rm -f $(DESTDIR)$(htmldir)/docs/$$f"; \ |
105 | 106 |
rm -f $(DESTDIR)$(htmldir)/docs/$$f; \ |
106 | 107 |
done |
107 | 108 |
|
108 | 109 |
.PHONY: update-external-tags |
1 | 1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library. |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2009 |
6 | 6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
7 | 7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES). |
8 | 8 |
* |
9 | 9 |
* Permission to use, modify and distribute this software is granted |
10 | 10 |
* provided that this copyright notice appears in all copies. For |
11 | 11 |
* precise terms see the accompanying LICENSE file. |
12 | 12 |
* |
13 | 13 |
* This software is provided "AS IS" with no warranty of any kind, |
14 | 14 |
* express or implied, and with no claim as to its suitability for any |
15 | 15 |
* purpose. |
16 | 16 |
* |
17 | 17 |
*/ |
18 | 18 |
|
19 | 19 |
namespace lemon { |
20 | 20 |
|
21 | 21 |
/** |
22 | 22 |
@defgroup datas Data Structures |
23 | 23 |
This group contains the several data structures implemented in LEMON. |
24 | 24 |
*/ |
25 | 25 |
|
26 | 26 |
/** |
27 | 27 |
@defgroup graphs Graph Structures |
28 | 28 |
@ingroup datas |
29 | 29 |
\brief Graph structures implemented in LEMON. |
30 | 30 |
|
31 | 31 |
The implementation of combinatorial algorithms heavily relies on |
32 | 32 |
efficient graph implementations. LEMON offers data structures which are |
33 | 33 |
planned to be easily used in an experimental phase of implementation studies, |
34 | 34 |
and thereafter the program code can be made efficient by small modifications. |
35 | 35 |
|
36 | 36 |
The most efficient implementation of diverse applications require the |
37 | 37 |
usage of different physical graph implementations. These differences |
38 | 38 |
appear in the size of graph we require to handle, memory or time usage |
39 | 39 |
limitations or in the set of operations through which the graph can be |
40 | 40 |
accessed. LEMON provides several physical graph structures to meet |
41 | 41 |
the diverging requirements of the possible users. In order to save on |
42 | 42 |
running time or on memory usage, some structures may fail to provide |
43 | 43 |
some graph features like arc/edge or node deletion. |
44 | 44 |
|
45 | 45 |
Alteration of standard containers need a very limited number of |
46 | 46 |
operations, these together satisfy the everyday requirements. |
47 | 47 |
In the case of graph structures, different operations are needed which do |
48 | 48 |
not alter the physical graph, but gives another view. If some nodes or |
49 | 49 |
arcs have to be hidden or the reverse oriented graph have to be used, then |
50 | 50 |
this is the case. It also may happen that in a flow implementation |
51 | 51 |
the residual graph can be accessed by another algorithm, or a node-set |
52 | 52 |
is to be shrunk for another algorithm. |
53 | 53 |
LEMON also provides a variety of graphs for these requirements called |
54 | 54 |
\ref graph_adaptors "graph adaptors". Adaptors cannot be used alone but only |
55 | 55 |
in conjunction with other graph representations. |
56 | 56 |
|
57 | 57 |
You are free to use the graph structure that fit your requirements |
58 | 58 |
the best, most graph algorithms and auxiliary data structures can be used |
59 | 59 |
with any graph structure. |
60 | 60 |
|
61 | 61 |
<b>See also:</b> \ref graph_concepts "Graph Structure Concepts". |
62 | 62 |
*/ |
63 | 63 |
|
64 | 64 |
/** |
65 | 65 |
@defgroup graph_adaptors Adaptor Classes for Graphs |
66 | 66 |
@ingroup graphs |
67 | 67 |
\brief Adaptor classes for digraphs and graphs |
68 | 68 |
|
69 | 69 |
This group contains several useful adaptor classes for digraphs and graphs. |
70 | 70 |
|
71 | 71 |
The main parts of LEMON are the different graph structures, generic |
72 | 72 |
graph algorithms, graph concepts, which couple them, and graph |
73 | 73 |
adaptors. While the previous notions are more or less clear, the |
74 | 74 |
latter one needs further explanation. Graph adaptors are graph classes |
75 | 75 |
which serve for considering graph structures in different ways. |
76 | 76 |
|
77 | 77 |
A short example makes this much clearer. Suppose that we have an |
78 | 78 |
instance \c g of a directed graph type, say ListDigraph and an algorithm |
79 | 79 |
\code |
80 | 80 |
template <typename Digraph> |
81 | 81 |
int algorithm(const Digraph&); |
82 | 82 |
\endcode |
83 | 83 |
is needed to run on the reverse oriented graph. It may be expensive |
84 | 84 |
(in time or in memory usage) to copy \c g with the reversed |
85 | 85 |
arcs. In this case, an adaptor class is used, which (according |
86 | 86 |
to LEMON \ref concepts::Digraph "digraph concepts") works as a digraph. |
87 | 87 |
The adaptor uses the original digraph structure and digraph operations when |
88 | 88 |
methods of the reversed oriented graph are called. This means that the adaptor |
89 | 89 |
have minor memory usage, and do not perform sophisticated algorithmic |
90 | 90 |
actions. The purpose of it is to give a tool for the cases when a |
91 | 91 |
graph have to be used in a specific alteration. If this alteration is |
92 | 92 |
obtained by a usual construction like filtering the node or the arc set or |
93 | 93 |
considering a new orientation, then an adaptor is worthwhile to use. |
94 | 94 |
To come back to the reverse oriented graph, in this situation |
95 | 95 |
\code |
96 | 96 |
template<typename Digraph> class ReverseDigraph; |
97 | 97 |
\endcode |
98 | 98 |
template class can be used. The code looks as follows |
99 | 99 |
\code |
100 | 100 |
ListDigraph g; |
101 | 101 |
ReverseDigraph<ListDigraph> rg(g); |
102 | 102 |
int result = algorithm(rg); |
103 | 103 |
\endcode |
104 | 104 |
During running the algorithm, the original digraph \c g is untouched. |
105 | 105 |
This techniques give rise to an elegant code, and based on stable |
106 | 106 |
graph adaptors, complex algorithms can be implemented easily. |
107 | 107 |
|
108 | 108 |
In flow, circulation and matching problems, the residual |
109 | 109 |
graph is of particular importance. Combining an adaptor implementing |
110 | 110 |
this with shortest path algorithms or minimum mean cycle algorithms, |
111 | 111 |
a range of weighted and cardinality optimization algorithms can be |
112 | 112 |
obtained. For other examples, the interested user is referred to the |
113 | 113 |
detailed documentation of particular adaptors. |
114 | 114 |
|
115 | 115 |
The behavior of graph adaptors can be very different. Some of them keep |
116 | 116 |
capabilities of the original graph while in other cases this would be |
117 | 117 |
meaningless. This means that the concepts that they meet depend |
118 | 118 |
on the graph adaptor, and the wrapped graph. |
119 | 119 |
For example, if an arc of a reversed digraph is deleted, this is carried |
120 | 120 |
out by deleting the corresponding arc of the original digraph, thus the |
121 | 121 |
adaptor modifies the original digraph. |
122 | 122 |
However in case of a residual digraph, this operation has no sense. |
123 | 123 |
|
124 | 124 |
Let us stand one more example here to simplify your work. |
125 | 125 |
ReverseDigraph has constructor |
126 | 126 |
\code |
127 | 127 |
ReverseDigraph(Digraph& digraph); |
128 | 128 |
\endcode |
129 | 129 |
This means that in a situation, when a <tt>const %ListDigraph&</tt> |
130 | 130 |
reference to a graph is given, then it have to be instantiated with |
131 | 131 |
<tt>Digraph=const %ListDigraph</tt>. |
132 | 132 |
\code |
133 | 133 |
int algorithm1(const ListDigraph& g) { |
134 | 134 |
ReverseDigraph<const ListDigraph> rg(g); |
135 | 135 |
return algorithm2(rg); |
136 | 136 |
} |
137 | 137 |
\endcode |
138 | 138 |
*/ |
139 | 139 |
|
140 | 140 |
/** |
141 | 141 |
@defgroup maps Maps |
142 | 142 |
@ingroup datas |
143 | 143 |
\brief Map structures implemented in LEMON. |
144 | 144 |
|
145 | 145 |
This group contains the map structures implemented in LEMON. |
146 | 146 |
|
147 | 147 |
LEMON provides several special purpose maps and map adaptors that e.g. combine |
148 | 148 |
new maps from existing ones. |
149 | 149 |
|
150 | 150 |
<b>See also:</b> \ref map_concepts "Map Concepts". |
151 | 151 |
*/ |
152 | 152 |
|
153 | 153 |
/** |
154 | 154 |
@defgroup graph_maps Graph Maps |
155 | 155 |
@ingroup maps |
156 | 156 |
\brief Special graph-related maps. |
157 | 157 |
|
158 | 158 |
This group contains maps that are specifically designed to assign |
159 | 159 |
values to the nodes and arcs/edges of graphs. |
160 | 160 |
|
161 | 161 |
If you are looking for the standard graph maps (\c NodeMap, \c ArcMap, |
162 | 162 |
\c EdgeMap), see the \ref graph_concepts "Graph Structure Concepts". |
163 | 163 |
*/ |
164 | 164 |
|
165 | 165 |
/** |
166 | 166 |
\defgroup map_adaptors Map Adaptors |
167 | 167 |
\ingroup maps |
168 | 168 |
\brief Tools to create new maps from existing ones |
169 | 169 |
|
170 | 170 |
This group contains map adaptors that are used to create "implicit" |
171 | 171 |
maps from other maps. |
172 | 172 |
|
173 | 173 |
Most of them are \ref concepts::ReadMap "read-only maps". |
174 | 174 |
They can make arithmetic and logical operations between one or two maps |
175 | 175 |
(negation, shifting, addition, multiplication, logical 'and', 'or', |
176 | 176 |
'not' etc.) or e.g. convert a map to another one of different Value type. |
177 | 177 |
|
178 | 178 |
The typical usage of this classes is passing implicit maps to |
179 | 179 |
algorithms. If a function type algorithm is called then the function |
180 | 180 |
type map adaptors can be used comfortable. For example let's see the |
181 | 181 |
usage of map adaptors with the \c graphToEps() function. |
182 | 182 |
\code |
183 | 183 |
Color nodeColor(int deg) { |
184 | 184 |
if (deg >= 2) { |
185 | 185 |
return Color(0.5, 0.0, 0.5); |
186 | 186 |
} else if (deg == 1) { |
187 | 187 |
return Color(1.0, 0.5, 1.0); |
188 | 188 |
} else { |
189 | 189 |
return Color(0.0, 0.0, 0.0); |
190 | 190 |
} |
191 | 191 |
} |
192 | 192 |
|
193 | 193 |
Digraph::NodeMap<int> degree_map(graph); |
194 | 194 |
|
195 | 195 |
graphToEps(graph, "graph.eps") |
196 | 196 |
.coords(coords).scaleToA4().undirected() |
197 | 197 |
.nodeColors(composeMap(functorToMap(nodeColor), degree_map)) |
198 | 198 |
.run(); |
199 | 199 |
\endcode |
200 | 200 |
The \c functorToMap() function makes an \c int to \c Color map from the |
201 | 201 |
\c nodeColor() function. The \c composeMap() compose the \c degree_map |
202 | 202 |
and the previously created map. The composed map is a proper function to |
203 | 203 |
get the color of each node. |
204 | 204 |
|
205 | 205 |
The usage with class type algorithms is little bit harder. In this |
206 | 206 |
case the function type map adaptors can not be used, because the |
207 | 207 |
function map adaptors give back temporary objects. |
208 | 208 |
\code |
209 | 209 |
Digraph graph; |
210 | 210 |
|
211 | 211 |
typedef Digraph::ArcMap<double> DoubleArcMap; |
212 | 212 |
DoubleArcMap length(graph); |
213 | 213 |
DoubleArcMap speed(graph); |
214 | 214 |
|
215 | 215 |
typedef DivMap<DoubleArcMap, DoubleArcMap> TimeMap; |
216 | 216 |
TimeMap time(length, speed); |
217 | 217 |
|
218 | 218 |
Dijkstra<Digraph, TimeMap> dijkstra(graph, time); |
219 | 219 |
dijkstra.run(source, target); |
220 | 220 |
\endcode |
221 | 221 |
We have a length map and a maximum speed map on the arcs of a digraph. |
222 | 222 |
The minimum time to pass the arc can be calculated as the division of |
223 | 223 |
the two maps which can be done implicitly with the \c DivMap template |
224 | 224 |
class. We use the implicit minimum time map as the length map of the |
225 | 225 |
\c Dijkstra algorithm. |
226 | 226 |
*/ |
227 | 227 |
|
228 | 228 |
/** |
229 | 229 |
@defgroup matrices Matrices |
230 | 230 |
@ingroup datas |
231 | 231 |
\brief Two dimensional data storages implemented in LEMON. |
232 | 232 |
|
233 | 233 |
This group contains two dimensional data storages implemented in LEMON. |
234 | 234 |
*/ |
235 | 235 |
|
236 | 236 |
/** |
237 | 237 |
@defgroup paths Path Structures |
238 | 238 |
@ingroup datas |
239 | 239 |
\brief %Path structures implemented in LEMON. |
240 | 240 |
|
241 | 241 |
This group contains the path structures implemented in LEMON. |
242 | 242 |
|
243 | 243 |
LEMON provides flexible data structures to work with paths. |
244 | 244 |
All of them have similar interfaces and they can be copied easily with |
245 | 245 |
assignment operators and copy constructors. This makes it easy and |
246 | 246 |
efficient to have e.g. the Dijkstra algorithm to store its result in |
247 | 247 |
any kind of path structure. |
248 | 248 |
|
249 | 249 |
\sa lemon::concepts::Path |
250 | 250 |
*/ |
251 | 251 |
|
252 | 252 |
/** |
253 | 253 |
@defgroup auxdat Auxiliary Data Structures |
254 | 254 |
@ingroup datas |
255 | 255 |
\brief Auxiliary data structures implemented in LEMON. |
256 | 256 |
|
257 | 257 |
This group contains some data structures implemented in LEMON in |
258 | 258 |
order to make it easier to implement combinatorial algorithms. |
259 | 259 |
*/ |
260 | 260 |
|
261 | 261 |
/** |
262 | 262 |
@defgroup algs Algorithms |
263 | 263 |
\brief This group contains the several algorithms |
264 | 264 |
implemented in LEMON. |
265 | 265 |
|
266 | 266 |
This group contains the several algorithms |
267 | 267 |
implemented in LEMON. |
268 | 268 |
*/ |
269 | 269 |
|
270 | 270 |
/** |
271 | 271 |
@defgroup search Graph Search |
272 | 272 |
@ingroup algs |
273 | 273 |
\brief Common graph search algorithms. |
274 | 274 |
|
275 | 275 |
This group contains the common graph search algorithms, namely |
276 | 276 |
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS). |
277 | 277 |
*/ |
278 | 278 |
|
279 | 279 |
/** |
280 | 280 |
@defgroup shortest_path Shortest Path Algorithms |
281 | 281 |
@ingroup algs |
282 | 282 |
\brief Algorithms for finding shortest paths. |
283 | 283 |
|
284 | 284 |
This group contains the algorithms for finding shortest paths in digraphs. |
285 | 285 |
|
286 | 286 |
- \ref Dijkstra algorithm for finding shortest paths from a source node |
287 | 287 |
when all arc lengths are non-negative. |
288 | 288 |
- \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths |
289 | 289 |
from a source node when arc lenghts can be either positive or negative, |
290 | 290 |
but the digraph should not contain directed cycles with negative total |
291 | 291 |
length. |
292 | 292 |
- \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms |
293 | 293 |
for solving the \e all-pairs \e shortest \e paths \e problem when arc |
294 | 294 |
lenghts can be either positive or negative, but the digraph should |
295 | 295 |
not contain directed cycles with negative total length. |
296 | 296 |
- \ref Suurballe A successive shortest path algorithm for finding |
297 | 297 |
arc-disjoint paths between two nodes having minimum total length. |
298 | 298 |
*/ |
299 | 299 |
|
300 | 300 |
/** |
301 | 301 |
@defgroup max_flow Maximum Flow Algorithms |
302 | 302 |
@ingroup algs |
303 | 303 |
\brief Algorithms for finding maximum flows. |
304 | 304 |
|
305 | 305 |
This group contains the algorithms for finding maximum flows and |
306 | 306 |
feasible circulations. |
307 | 307 |
|
308 | 308 |
The \e maximum \e flow \e problem is to find a flow of maximum value between |
309 | 309 |
a single source and a single target. Formally, there is a \f$G=(V,A)\f$ |
310 | 310 |
digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and |
311 | 311 |
\f$s, t \in V\f$ source and target nodes. |
312 | 312 |
A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the |
313 | 313 |
following optimization problem. |
314 | 314 |
|
315 | 315 |
\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f] |
316 | 316 |
\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu) |
317 | 317 |
\quad \forall u\in V\setminus\{s,t\} \f] |
318 | 318 |
\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f] |
319 | 319 |
|
320 | 320 |
LEMON contains several algorithms for solving maximum flow problems: |
321 | 321 |
- \ref EdmondsKarp Edmonds-Karp algorithm. |
322 | 322 |
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm. |
323 | 323 |
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees. |
324 | 324 |
- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees. |
325 | 325 |
|
326 | 326 |
In most cases the \ref Preflow "Preflow" algorithm provides the |
327 | 327 |
fastest method for computing a maximum flow. All implementations |
328 | 328 |
also provide functions to query the minimum cut, which is the dual |
329 | 329 |
problem of maximum flow. |
330 | 330 |
|
331 | 331 |
\ref Circulation is a preflow push-relabel algorithm implemented directly |
332 | 332 |
for finding feasible circulations, which is a somewhat different problem, |
333 | 333 |
but it is strongly related to maximum flow. |
334 | 334 |
For more information, see \ref Circulation. |
335 | 335 |
*/ |
336 | 336 |
|
337 | 337 |
/** |
338 |
@defgroup |
|
338 |
@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms |
|
339 | 339 |
@ingroup algs |
340 | 340 |
|
341 | 341 |
\brief Algorithms for finding minimum cost flows and circulations. |
342 | 342 |
|
343 | 343 |
This group contains the algorithms for finding minimum cost flows and |
344 |
circulations. |
|
344 |
circulations. For more information about this problem and its dual |
|
345 |
solution see \ref min_cost_flow "Minimum Cost Flow Problem". |
|
345 | 346 |
|
346 |
The \e minimum \e cost \e flow \e problem is to find a feasible flow of |
|
347 |
minimum total cost from a set of supply nodes to a set of demand nodes |
|
348 |
in a network with capacity constraints (lower and upper bounds) |
|
349 |
and arc costs. |
|
350 |
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$, |
|
351 |
\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and |
|
352 |
upper bounds for the flow values on the arcs, for which |
|
353 |
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, |
|
354 |
\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow |
|
355 |
on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the |
|
356 |
signed supply values of the nodes. |
|
357 |
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ |
|
358 |
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with |
|
359 |
\f$-sup(u)\f$ demand. |
|
360 |
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution |
|
361 |
of the following optimization problem. |
|
362 |
|
|
363 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
364 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq |
|
365 |
sup(u) \quad \forall u\in V \f] |
|
366 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
367 |
|
|
368 |
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be |
|
369 |
zero or negative in order to have a feasible solution (since the sum |
|
370 |
of the expressions on the left-hand side of the inequalities is zero). |
|
371 |
It means that the total demand must be greater or equal to the total |
|
372 |
supply and all the supplies have to be carried out from the supply nodes, |
|
373 |
but there could be demands that are not satisfied. |
|
374 |
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand |
|
375 |
constraints have to be satisfied with equality, i.e. all demands |
|
376 |
have to be satisfied and all supplies have to be used. |
|
377 |
|
|
378 |
If you need the opposite inequalities in the supply/demand constraints |
|
379 |
(i.e. the total demand is less than the total supply and all the demands |
|
380 |
have to be satisfied while there could be supplies that are not used), |
|
381 |
then you could easily transform the problem to the above form by reversing |
|
382 |
the direction of the arcs and taking the negative of the supply values |
|
383 |
(e.g. using \ref ReverseDigraph and \ref NegMap adaptors). |
|
384 |
However \ref NetworkSimplex algorithm also supports this form directly |
|
385 |
for the sake of convenience. |
|
386 |
|
|
387 |
A feasible solution for this problem can be found using \ref Circulation. |
|
388 |
|
|
389 |
Note that the above formulation is actually more general than the usual |
|
390 |
definition of the minimum cost flow problem, in which strict equalities |
|
391 |
are required in the supply/demand contraints, i.e. |
|
392 |
|
|
393 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) = |
|
394 |
sup(u) \quad \forall u\in V. \f] |
|
395 |
|
|
396 |
However if the sum of the supply values is zero, then these two problems |
|
397 |
are equivalent. So if you need the equality form, you have to ensure this |
|
398 |
additional contraint for the algorithms. |
|
399 |
|
|
400 |
The dual solution of the minimum cost flow problem is represented by node |
|
401 |
potentials \f$\pi: V\rightarrow\mathbf{Z}\f$. |
|
402 |
An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem |
|
403 |
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$ |
|
404 |
node potentials the following \e complementary \e slackness optimality |
|
405 |
conditions hold. |
|
406 |
|
|
407 |
- For all \f$uv\in A\f$ arcs: |
|
408 |
- if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; |
|
409 |
- if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; |
|
410 |
- if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. |
|
411 |
- For all \f$u\in V\f$ nodes: |
|
412 |
- if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$, |
|
413 |
then \f$\pi(u)=0\f$. |
|
414 |
|
|
415 |
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc |
|
416 |
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e. |
|
417 |
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f] |
|
418 |
|
|
419 |
All algorithms provide dual solution (node potentials) as well, |
|
420 |
if an optimal flow is found. |
|
421 |
|
|
422 |
LEMON contains several algorithms for |
|
347 |
LEMON contains several algorithms for this problem. |
|
423 | 348 |
- \ref NetworkSimplex Primal Network Simplex algorithm with various |
424 | 349 |
pivot strategies. |
425 | 350 |
- \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on |
426 | 351 |
cost scaling. |
427 | 352 |
- \ref CapacityScaling Successive Shortest %Path algorithm with optional |
428 | 353 |
capacity scaling. |
429 | 354 |
- \ref CancelAndTighten The Cancel and Tighten algorithm. |
430 | 355 |
- \ref CycleCanceling Cycle-Canceling algorithms. |
431 | 356 |
|
432 |
Most of these implementations support the general inequality form of the |
|
433 |
minimum cost flow problem, but CancelAndTighten and CycleCanceling |
|
434 |
only support the equality form due to the primal method they use. |
|
435 |
|
|
436 | 357 |
In general NetworkSimplex is the most efficient implementation, |
437 | 358 |
but in special cases other algorithms could be faster. |
438 | 359 |
For example, if the total supply and/or capacities are rather small, |
439 | 360 |
CapacityScaling is usually the fastest algorithm (without effective scaling). |
440 | 361 |
*/ |
441 | 362 |
|
442 | 363 |
/** |
443 | 364 |
@defgroup min_cut Minimum Cut Algorithms |
444 | 365 |
@ingroup algs |
445 | 366 |
|
446 | 367 |
\brief Algorithms for finding minimum cut in graphs. |
447 | 368 |
|
448 | 369 |
This group contains the algorithms for finding minimum cut in graphs. |
449 | 370 |
|
450 | 371 |
The \e minimum \e cut \e problem is to find a non-empty and non-complete |
451 | 372 |
\f$X\f$ subset of the nodes with minimum overall capacity on |
452 | 373 |
outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a |
453 | 374 |
\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum |
454 | 375 |
cut is the \f$X\f$ solution of the next optimization problem: |
455 | 376 |
|
456 | 377 |
\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}} |
457 | 378 |
\sum_{uv\in A, u\in X, v\not\in X}cap(uv) \f] |
458 | 379 |
|
459 | 380 |
LEMON contains several algorithms related to minimum cut problems: |
460 | 381 |
|
461 | 382 |
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut |
462 | 383 |
in directed graphs. |
463 | 384 |
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for |
464 | 385 |
calculating minimum cut in undirected graphs. |
465 | 386 |
- \ref GomoryHu "Gomory-Hu tree computation" for calculating |
466 | 387 |
all-pairs minimum cut in undirected graphs. |
467 | 388 |
|
468 | 389 |
If you want to find minimum cut just between two distinict nodes, |
469 | 390 |
see the \ref max_flow "maximum flow problem". |
470 | 391 |
*/ |
471 | 392 |
|
472 | 393 |
/** |
473 | 394 |
@defgroup graph_properties Connectivity and Other Graph Properties |
474 | 395 |
@ingroup algs |
475 | 396 |
\brief Algorithms for discovering the graph properties |
476 | 397 |
|
477 | 398 |
This group contains the algorithms for discovering the graph properties |
478 | 399 |
like connectivity, bipartiteness, euler property, simplicity etc. |
479 | 400 |
|
480 | 401 |
\image html edge_biconnected_components.png |
481 | 402 |
\image latex edge_biconnected_components.eps "bi-edge-connected components" width=\textwidth |
482 | 403 |
*/ |
483 | 404 |
|
484 | 405 |
/** |
485 | 406 |
@defgroup planar Planarity Embedding and Drawing |
486 | 407 |
@ingroup algs |
487 | 408 |
\brief Algorithms for planarity checking, embedding and drawing |
488 | 409 |
|
489 | 410 |
This group contains the algorithms for planarity checking, |
490 | 411 |
embedding and drawing. |
491 | 412 |
|
492 | 413 |
\image html planar.png |
493 | 414 |
\image latex planar.eps "Plane graph" width=\textwidth |
494 | 415 |
*/ |
495 | 416 |
|
496 | 417 |
/** |
497 | 418 |
@defgroup matching Matching Algorithms |
498 | 419 |
@ingroup algs |
499 | 420 |
\brief Algorithms for finding matchings in graphs and bipartite graphs. |
500 | 421 |
|
501 | 422 |
This group contains the algorithms for calculating |
502 | 423 |
matchings in graphs and bipartite graphs. The general matching problem is |
503 | 424 |
finding a subset of the edges for which each node has at most one incident |
504 | 425 |
edge. |
505 | 426 |
|
506 | 427 |
There are several different algorithms for calculate matchings in |
507 | 428 |
graphs. The matching problems in bipartite graphs are generally |
508 | 429 |
easier than in general graphs. The goal of the matching optimization |
509 | 430 |
can be finding maximum cardinality, maximum weight or minimum cost |
510 | 431 |
matching. The search can be constrained to find perfect or |
511 | 432 |
maximum cardinality matching. |
512 | 433 |
|
513 | 434 |
The matching algorithms implemented in LEMON: |
514 | 435 |
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm |
515 | 436 |
for calculating maximum cardinality matching in bipartite graphs. |
516 | 437 |
- \ref PrBipartiteMatching Push-relabel algorithm |
517 | 438 |
for calculating maximum cardinality matching in bipartite graphs. |
518 | 439 |
- \ref MaxWeightedBipartiteMatching |
519 | 440 |
Successive shortest path algorithm for calculating maximum weighted |
520 | 441 |
matching and maximum weighted bipartite matching in bipartite graphs. |
521 | 442 |
- \ref MinCostMaxBipartiteMatching |
522 | 443 |
Successive shortest path algorithm for calculating minimum cost maximum |
523 | 444 |
matching in bipartite graphs. |
524 | 445 |
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating |
525 | 446 |
maximum cardinality matching in general graphs. |
526 | 447 |
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating |
527 | 448 |
maximum weighted matching in general graphs. |
528 | 449 |
- \ref MaxWeightedPerfectMatching |
529 | 450 |
Edmond's blossom shrinking algorithm for calculating maximum weighted |
530 | 451 |
perfect matching in general graphs. |
531 | 452 |
|
532 | 453 |
\image html bipartite_matching.png |
533 | 454 |
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth |
534 | 455 |
*/ |
535 | 456 |
|
536 | 457 |
/** |
537 | 458 |
@defgroup spantree Minimum Spanning Tree Algorithms |
538 | 459 |
@ingroup algs |
539 | 460 |
\brief Algorithms for finding minimum cost spanning trees and arborescences. |
540 | 461 |
|
541 | 462 |
This group contains the algorithms for finding minimum cost spanning |
542 | 463 |
trees and arborescences. |
543 | 464 |
*/ |
544 | 465 |
|
545 | 466 |
/** |
546 | 467 |
@defgroup auxalg Auxiliary Algorithms |
547 | 468 |
@ingroup algs |
548 | 469 |
\brief Auxiliary algorithms implemented in LEMON. |
549 | 470 |
|
550 | 471 |
This group contains some algorithms implemented in LEMON |
551 | 472 |
in order to make it easier to implement complex algorithms. |
552 | 473 |
*/ |
553 | 474 |
|
554 | 475 |
/** |
555 | 476 |
@defgroup approx Approximation Algorithms |
556 | 477 |
@ingroup algs |
557 | 478 |
\brief Approximation algorithms. |
558 | 479 |
|
559 | 480 |
This group contains the approximation and heuristic algorithms |
560 | 481 |
implemented in LEMON. |
561 | 482 |
*/ |
562 | 483 |
|
563 | 484 |
/** |
564 | 485 |
@defgroup gen_opt_group General Optimization Tools |
565 | 486 |
\brief This group contains some general optimization frameworks |
566 | 487 |
implemented in LEMON. |
567 | 488 |
|
568 | 489 |
This group contains some general optimization frameworks |
569 | 490 |
implemented in LEMON. |
570 | 491 |
*/ |
571 | 492 |
|
572 | 493 |
/** |
573 | 494 |
@defgroup lp_group Lp and Mip Solvers |
574 | 495 |
@ingroup gen_opt_group |
575 | 496 |
\brief Lp and Mip solver interfaces for LEMON. |
576 | 497 |
|
577 | 498 |
This group contains Lp and Mip solver interfaces for LEMON. The |
578 | 499 |
various LP solvers could be used in the same manner with this |
579 | 500 |
interface. |
580 | 501 |
*/ |
581 | 502 |
|
582 | 503 |
/** |
583 | 504 |
@defgroup lp_utils Tools for Lp and Mip Solvers |
584 | 505 |
@ingroup lp_group |
585 | 506 |
\brief Helper tools to the Lp and Mip solvers. |
586 | 507 |
|
587 | 508 |
This group adds some helper tools to general optimization framework |
588 | 509 |
implemented in LEMON. |
589 | 510 |
*/ |
590 | 511 |
|
591 | 512 |
/** |
592 | 513 |
@defgroup metah Metaheuristics |
593 | 514 |
@ingroup gen_opt_group |
594 | 515 |
\brief Metaheuristics for LEMON library. |
595 | 516 |
|
596 | 517 |
This group contains some metaheuristic optimization tools. |
597 | 518 |
*/ |
598 | 519 |
|
599 | 520 |
/** |
600 | 521 |
@defgroup utils Tools and Utilities |
601 | 522 |
\brief Tools and utilities for programming in LEMON |
602 | 523 |
|
603 | 524 |
Tools and utilities for programming in LEMON. |
604 | 525 |
*/ |
605 | 526 |
|
606 | 527 |
/** |
607 | 528 |
@defgroup gutils Basic Graph Utilities |
608 | 529 |
@ingroup utils |
609 | 530 |
\brief Simple basic graph utilities. |
610 | 531 |
|
611 | 532 |
This group contains some simple basic graph utilities. |
612 | 533 |
*/ |
613 | 534 |
|
614 | 535 |
/** |
615 | 536 |
@defgroup misc Miscellaneous Tools |
616 | 537 |
@ingroup utils |
617 | 538 |
\brief Tools for development, debugging and testing. |
618 | 539 |
|
619 | 540 |
This group contains several useful tools for development, |
620 | 541 |
debugging and testing. |
621 | 542 |
*/ |
622 | 543 |
|
623 | 544 |
/** |
624 | 545 |
@defgroup timecount Time Measuring and Counting |
625 | 546 |
@ingroup misc |
626 | 547 |
\brief Simple tools for measuring the performance of algorithms. |
627 | 548 |
|
628 | 549 |
This group contains simple tools for measuring the performance |
629 | 550 |
of algorithms. |
630 | 551 |
*/ |
631 | 552 |
|
632 | 553 |
/** |
633 | 554 |
@defgroup exceptions Exceptions |
634 | 555 |
@ingroup utils |
635 | 556 |
\brief Exceptions defined in LEMON. |
636 | 557 |
|
637 | 558 |
This group contains the exceptions defined in LEMON. |
638 | 559 |
*/ |
639 | 560 |
|
640 | 561 |
/** |
641 | 562 |
@defgroup io_group Input-Output |
642 | 563 |
\brief Graph Input-Output methods |
643 | 564 |
|
644 | 565 |
This group contains the tools for importing and exporting graphs |
645 | 566 |
and graph related data. Now it supports the \ref lgf-format |
646 | 567 |
"LEMON Graph Format", the \c DIMACS format and the encapsulated |
647 | 568 |
postscript (EPS) format. |
648 | 569 |
*/ |
649 | 570 |
|
650 | 571 |
/** |
651 | 572 |
@defgroup lemon_io LEMON Graph Format |
652 | 573 |
@ingroup io_group |
653 | 574 |
\brief Reading and writing LEMON Graph Format. |
654 | 575 |
|
655 | 576 |
This group contains methods for reading and writing |
656 | 577 |
\ref lgf-format "LEMON Graph Format". |
657 | 578 |
*/ |
658 | 579 |
|
659 | 580 |
/** |
660 | 581 |
@defgroup eps_io Postscript Exporting |
661 | 582 |
@ingroup io_group |
662 | 583 |
\brief General \c EPS drawer and graph exporter |
663 | 584 |
|
664 | 585 |
This group contains general \c EPS drawing methods and special |
665 | 586 |
graph exporting tools. |
666 | 587 |
*/ |
667 | 588 |
|
668 | 589 |
/** |
669 | 590 |
@defgroup dimacs_group DIMACS format |
670 | 591 |
@ingroup io_group |
671 | 592 |
\brief Read and write files in DIMACS format |
672 | 593 |
|
673 | 594 |
Tools to read a digraph from or write it to a file in DIMACS format data. |
674 | 595 |
*/ |
675 | 596 |
|
676 | 597 |
/** |
677 | 598 |
@defgroup nauty_group NAUTY Format |
678 | 599 |
@ingroup io_group |
679 | 600 |
\brief Read \e Nauty format |
680 | 601 |
|
681 | 602 |
Tool to read graphs from \e Nauty format data. |
682 | 603 |
*/ |
683 | 604 |
|
684 | 605 |
/** |
685 | 606 |
@defgroup concept Concepts |
686 | 607 |
\brief Skeleton classes and concept checking classes |
687 | 608 |
|
688 | 609 |
This group contains the data/algorithm skeletons and concept checking |
689 | 610 |
classes implemented in LEMON. |
690 | 611 |
|
691 | 612 |
The purpose of the classes in this group is fourfold. |
692 | 613 |
|
693 | 614 |
- These classes contain the documentations of the %concepts. In order |
694 | 615 |
to avoid document multiplications, an implementation of a concept |
695 | 616 |
simply refers to the corresponding concept class. |
696 | 617 |
|
697 | 618 |
- These classes declare every functions, <tt>typedef</tt>s etc. an |
698 | 619 |
implementation of the %concepts should provide, however completely |
699 | 620 |
without implementations and real data structures behind the |
700 | 621 |
interface. On the other hand they should provide nothing else. All |
701 | 622 |
the algorithms working on a data structure meeting a certain concept |
702 | 623 |
should compile with these classes. (Though it will not run properly, |
703 | 624 |
of course.) In this way it is easily to check if an algorithm |
704 | 625 |
doesn't use any extra feature of a certain implementation. |
705 | 626 |
|
706 | 627 |
- The concept descriptor classes also provide a <em>checker class</em> |
707 | 628 |
that makes it possible to check whether a certain implementation of a |
708 | 629 |
concept indeed provides all the required features. |
709 | 630 |
|
710 | 631 |
- Finally, They can serve as a skeleton of a new implementation of a concept. |
711 | 632 |
*/ |
712 | 633 |
|
713 | 634 |
/** |
714 | 635 |
@defgroup graph_concepts Graph Structure Concepts |
715 | 636 |
@ingroup concept |
716 | 637 |
\brief Skeleton and concept checking classes for graph structures |
717 | 638 |
|
718 | 639 |
This group contains the skeletons and concept checking classes of LEMON's |
719 | 640 |
graph structures and helper classes used to implement these. |
720 | 641 |
*/ |
721 | 642 |
|
722 | 643 |
/** |
723 | 644 |
@defgroup map_concepts Map Concepts |
724 | 645 |
@ingroup concept |
725 | 646 |
\brief Skeleton and concept checking classes for maps |
726 | 647 |
|
727 | 648 |
This group contains the skeletons and concept checking classes of maps. |
728 | 649 |
*/ |
729 | 650 |
|
730 | 651 |
/** |
731 | 652 |
\anchor demoprograms |
732 | 653 |
|
733 | 654 |
@defgroup demos Demo Programs |
734 | 655 |
|
735 | 656 |
Some demo programs are listed here. Their full source codes can be found in |
736 | 657 |
the \c demo subdirectory of the source tree. |
737 | 658 |
|
738 | 659 |
In order to compile them, use the <tt>make demo</tt> or the |
739 | 660 |
<tt>make check</tt> commands. |
740 | 661 |
*/ |
741 | 662 |
|
742 | 663 |
/** |
743 | 664 |
@defgroup tools Standalone Utility Applications |
744 | 665 |
|
745 | 666 |
Some utility applications are listed here. |
746 | 667 |
|
747 | 668 |
The standard compilation procedure (<tt>./configure;make</tt>) will compile |
748 | 669 |
them, as well. |
749 | 670 |
*/ |
750 | 671 |
|
751 | 672 |
} |
1 | 1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library. |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2009 |
6 | 6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
7 | 7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES). |
8 | 8 |
* |
9 | 9 |
* Permission to use, modify and distribute this software is granted |
10 | 10 |
* provided that this copyright notice appears in all copies. For |
11 | 11 |
* precise terms see the accompanying LICENSE file. |
12 | 12 |
* |
13 | 13 |
* This software is provided "AS IS" with no warranty of any kind, |
14 | 14 |
* express or implied, and with no claim as to its suitability for any |
15 | 15 |
* purpose. |
16 | 16 |
* |
17 | 17 |
*/ |
18 | 18 |
|
19 | 19 |
#ifndef LEMON_NETWORK_SIMPLEX_H |
20 | 20 |
#define LEMON_NETWORK_SIMPLEX_H |
21 | 21 |
|
22 |
/// \ingroup |
|
22 |
/// \ingroup min_cost_flow_algs |
|
23 | 23 |
/// |
24 | 24 |
/// \file |
25 | 25 |
/// \brief Network Simplex algorithm for finding a minimum cost flow. |
26 | 26 |
|
27 | 27 |
#include <vector> |
28 | 28 |
#include <limits> |
29 | 29 |
#include <algorithm> |
30 | 30 |
|
31 | 31 |
#include <lemon/core.h> |
32 | 32 |
#include <lemon/math.h> |
33 | 33 |
|
34 | 34 |
namespace lemon { |
35 | 35 |
|
36 |
/// \addtogroup |
|
36 |
/// \addtogroup min_cost_flow_algs |
|
37 | 37 |
/// @{ |
38 | 38 |
|
39 | 39 |
/// \brief Implementation of the primal Network Simplex algorithm |
40 | 40 |
/// for finding a \ref min_cost_flow "minimum cost flow". |
41 | 41 |
/// |
42 | 42 |
/// \ref NetworkSimplex implements the primal Network Simplex algorithm |
43 | 43 |
/// for finding a \ref min_cost_flow "minimum cost flow". |
44 | 44 |
/// This algorithm is a specialized version of the linear programming |
45 | 45 |
/// simplex method directly for the minimum cost flow problem. |
46 | 46 |
/// It is one of the most efficient solution methods. |
47 | 47 |
/// |
48 | 48 |
/// In general this class is the fastest implementation available |
49 | 49 |
/// in LEMON for the minimum cost flow problem. |
50 | 50 |
/// Moreover it supports both directions of the supply/demand inequality |
51 | 51 |
/// constraints. For more information see \ref SupplyType. |
52 | 52 |
/// |
53 | 53 |
/// Most of the parameters of the problem (except for the digraph) |
54 | 54 |
/// can be given using separate functions, and the algorithm can be |
55 | 55 |
/// executed using the \ref run() function. If some parameters are not |
56 | 56 |
/// specified, then default values will be used. |
57 | 57 |
/// |
58 | 58 |
/// \tparam GR The digraph type the algorithm runs on. |
59 | 59 |
/// \tparam V The value type used for flow amounts, capacity bounds |
60 | 60 |
/// and supply values in the algorithm. By default it is \c int. |
61 | 61 |
/// \tparam C The value type used for costs and potentials in the |
62 | 62 |
/// algorithm. By default it is the same as \c V. |
63 | 63 |
/// |
64 | 64 |
/// \warning Both value types must be signed and all input data must |
65 | 65 |
/// be integer. |
66 | 66 |
/// |
67 | 67 |
/// \note %NetworkSimplex provides five different pivot rule |
68 | 68 |
/// implementations, from which the most efficient one is used |
69 | 69 |
/// by default. For more information see \ref PivotRule. |
70 | 70 |
template <typename GR, typename V = int, typename C = V> |
71 | 71 |
class NetworkSimplex |
72 | 72 |
{ |
73 | 73 |
public: |
74 | 74 |
|
75 | 75 |
/// The type of the flow amounts, capacity bounds and supply values |
76 | 76 |
typedef V Value; |
77 | 77 |
/// The type of the arc costs |
78 | 78 |
typedef C Cost; |
79 | 79 |
|
80 | 80 |
public: |
81 | 81 |
|
82 | 82 |
/// \brief Problem type constants for the \c run() function. |
83 | 83 |
/// |
84 | 84 |
/// Enum type containing the problem type constants that can be |
85 | 85 |
/// returned by the \ref run() function of the algorithm. |
86 | 86 |
enum ProblemType { |
87 | 87 |
/// The problem has no feasible solution (flow). |
88 | 88 |
INFEASIBLE, |
89 | 89 |
/// The problem has optimal solution (i.e. it is feasible and |
90 | 90 |
/// bounded), and the algorithm has found optimal flow and node |
91 | 91 |
/// potentials (primal and dual solutions). |
92 | 92 |
OPTIMAL, |
93 | 93 |
/// The objective function of the problem is unbounded, i.e. |
94 | 94 |
/// there is a directed cycle having negative total cost and |
95 | 95 |
/// infinite upper bound. |
96 | 96 |
UNBOUNDED |
97 | 97 |
}; |
98 | 98 |
|
99 | 99 |
/// \brief Constants for selecting the type of the supply constraints. |
100 | 100 |
/// |
101 | 101 |
/// Enum type containing constants for selecting the supply type, |
102 | 102 |
/// i.e. the direction of the inequalities in the supply/demand |
103 | 103 |
/// constraints of the \ref min_cost_flow "minimum cost flow problem". |
104 | 104 |
/// |
105 |
/// The default supply type is \c GEQ, since this form is supported |
|
106 |
/// by other minimum cost flow algorithms and the \ref Circulation |
|
107 |
/// algorithm, as well. |
|
108 |
/// The \c LEQ problem type can be selected using the \ref supplyType() |
|
109 |
/// function. |
|
110 |
/// |
|
111 |
/// |
|
105 |
/// The default supply type is \c GEQ, the \c LEQ type can be |
|
106 |
/// selected using \ref supplyType(). |
|
107 |
/// The equality form is a special case of both supply types. |
|
112 | 108 |
enum SupplyType { |
113 |
|
|
114 | 109 |
/// This option means that there are <em>"greater or equal"</em> |
115 |
/// supply/demand constraints in the definition, i.e. the exact |
|
116 |
/// formulation of the problem is the following. |
|
117 |
/** |
|
118 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
119 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq |
|
120 |
sup(u) \quad \forall u\in V \f] |
|
121 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
122 |
*/ |
|
123 |
/// It means that the total demand must be greater or equal to the |
|
124 |
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or |
|
125 |
/// negative) and all the supplies have to be carried out from |
|
126 |
/// the supply nodes, but there could be demands that are not |
|
127 |
/// |
|
110 |
/// supply/demand constraints in the definition of the problem. |
|
128 | 111 |
GEQ, |
129 |
/// It is just an alias for the \c GEQ option. |
|
130 |
CARRY_SUPPLIES = GEQ, |
|
131 |
|
|
132 | 112 |
/// This option means that there are <em>"less or equal"</em> |
133 |
/// supply/demand constraints in the definition, i.e. the exact |
|
134 |
/// formulation of the problem is the following. |
|
135 |
/** |
|
136 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
|
137 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq |
|
138 |
sup(u) \quad \forall u\in V \f] |
|
139 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
|
140 |
*/ |
|
141 |
/// It means that the total demand must be less or equal to the |
|
142 |
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or |
|
143 |
/// positive) and all the demands have to be satisfied, but there |
|
144 |
/// could be supplies that are not carried out from the supply |
|
145 |
/// nodes. |
|
146 |
LEQ, |
|
147 |
/// It is just an alias for the \c LEQ option. |
|
148 |
SATISFY_DEMANDS = LEQ |
|
113 |
/// supply/demand constraints in the definition of the problem. |
|
114 |
LEQ |
|
149 | 115 |
}; |
150 | 116 |
|
151 | 117 |
/// \brief Constants for selecting the pivot rule. |
152 | 118 |
/// |
153 | 119 |
/// Enum type containing constants for selecting the pivot rule for |
154 | 120 |
/// the \ref run() function. |
155 | 121 |
/// |
156 | 122 |
/// \ref NetworkSimplex provides five different pivot rule |
157 | 123 |
/// implementations that significantly affect the running time |
158 | 124 |
/// of the algorithm. |
159 | 125 |
/// By default \ref BLOCK_SEARCH "Block Search" is used, which |
160 | 126 |
/// proved to be the most efficient and the most robust on various |
161 | 127 |
/// test inputs according to our benchmark tests. |
162 | 128 |
/// However another pivot rule can be selected using the \ref run() |
163 | 129 |
/// function with the proper parameter. |
164 | 130 |
enum PivotRule { |
165 | 131 |
|
166 | 132 |
/// The First Eligible pivot rule. |
167 | 133 |
/// The next eligible arc is selected in a wraparound fashion |
168 | 134 |
/// in every iteration. |
169 | 135 |
FIRST_ELIGIBLE, |
170 | 136 |
|
171 | 137 |
/// The Best Eligible pivot rule. |
172 | 138 |
/// The best eligible arc is selected in every iteration. |
173 | 139 |
BEST_ELIGIBLE, |
174 | 140 |
|
175 | 141 |
/// The Block Search pivot rule. |
176 | 142 |
/// A specified number of arcs are examined in every iteration |
177 | 143 |
/// in a wraparound fashion and the best eligible arc is selected |
178 | 144 |
/// from this block. |
179 | 145 |
BLOCK_SEARCH, |
180 | 146 |
|
181 | 147 |
/// The Candidate List pivot rule. |
182 | 148 |
/// In a major iteration a candidate list is built from eligible arcs |
183 | 149 |
/// in a wraparound fashion and in the following minor iterations |
184 | 150 |
/// the best eligible arc is selected from this list. |
185 | 151 |
CANDIDATE_LIST, |
186 | 152 |
|
187 | 153 |
/// The Altering Candidate List pivot rule. |
188 | 154 |
/// It is a modified version of the Candidate List method. |
189 | 155 |
/// It keeps only the several best eligible arcs from the former |
190 | 156 |
/// candidate list and extends this list in every iteration. |
191 | 157 |
ALTERING_LIST |
192 | 158 |
}; |
193 | 159 |
|
194 | 160 |
private: |
195 | 161 |
|
196 | 162 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
197 | 163 |
|
198 | 164 |
typedef std::vector<Arc> ArcVector; |
199 | 165 |
typedef std::vector<Node> NodeVector; |
200 | 166 |
typedef std::vector<int> IntVector; |
201 | 167 |
typedef std::vector<bool> BoolVector; |
202 | 168 |
typedef std::vector<Value> ValueVector; |
203 | 169 |
typedef std::vector<Cost> CostVector; |
204 | 170 |
|
205 | 171 |
// State constants for arcs |
206 | 172 |
enum ArcStateEnum { |
207 | 173 |
STATE_UPPER = -1, |
208 | 174 |
STATE_TREE = 0, |
209 | 175 |
STATE_LOWER = 1 |
210 | 176 |
}; |
211 | 177 |
|
212 | 178 |
private: |
213 | 179 |
|
214 | 180 |
// Data related to the underlying digraph |
215 | 181 |
const GR &_graph; |
216 | 182 |
int _node_num; |
217 | 183 |
int _arc_num; |
184 |
int _all_arc_num; |
|
185 |
int _search_arc_num; |
|
218 | 186 |
|
219 | 187 |
// Parameters of the problem |
220 | 188 |
bool _have_lower; |
221 | 189 |
SupplyType _stype; |
222 | 190 |
Value _sum_supply; |
223 | 191 |
|
224 | 192 |
// Data structures for storing the digraph |
225 | 193 |
IntNodeMap _node_id; |
226 | 194 |
IntArcMap _arc_id; |
227 | 195 |
IntVector _source; |
228 | 196 |
IntVector _target; |
229 | 197 |
|
230 | 198 |
// Node and arc data |
231 | 199 |
ValueVector _lower; |
232 | 200 |
ValueVector _upper; |
233 | 201 |
ValueVector _cap; |
234 | 202 |
CostVector _cost; |
235 | 203 |
ValueVector _supply; |
236 | 204 |
ValueVector _flow; |
237 | 205 |
CostVector _pi; |
238 | 206 |
|
239 | 207 |
// Data for storing the spanning tree structure |
240 | 208 |
IntVector _parent; |
241 | 209 |
IntVector _pred; |
242 | 210 |
IntVector _thread; |
243 | 211 |
IntVector _rev_thread; |
244 | 212 |
IntVector _succ_num; |
245 | 213 |
IntVector _last_succ; |
246 | 214 |
IntVector _dirty_revs; |
247 | 215 |
BoolVector _forward; |
248 | 216 |
IntVector _state; |
249 | 217 |
int _root; |
250 | 218 |
|
251 | 219 |
// Temporary data used in the current pivot iteration |
252 | 220 |
int in_arc, join, u_in, v_in, u_out, v_out; |
253 | 221 |
int first, second, right, last; |
254 | 222 |
int stem, par_stem, new_stem; |
255 | 223 |
Value delta; |
256 | 224 |
|
257 | 225 |
public: |
258 | 226 |
|
259 | 227 |
/// \brief Constant for infinite upper bounds (capacities). |
260 | 228 |
/// |
261 | 229 |
/// Constant for infinite upper bounds (capacities). |
262 | 230 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
263 | 231 |
/// \c std::numeric_limits<Value>::max() otherwise. |
264 | 232 |
const Value INF; |
265 | 233 |
|
266 | 234 |
private: |
267 | 235 |
|
268 | 236 |
// Implementation of the First Eligible pivot rule |
269 | 237 |
class FirstEligiblePivotRule |
270 | 238 |
{ |
271 | 239 |
private: |
272 | 240 |
|
273 | 241 |
// References to the NetworkSimplex class |
274 | 242 |
const IntVector &_source; |
275 | 243 |
const IntVector &_target; |
276 | 244 |
const CostVector &_cost; |
277 | 245 |
const IntVector &_state; |
278 | 246 |
const CostVector &_pi; |
279 | 247 |
int &_in_arc; |
280 |
int |
|
248 |
int _search_arc_num; |
|
281 | 249 |
|
282 | 250 |
// Pivot rule data |
283 | 251 |
int _next_arc; |
284 | 252 |
|
285 | 253 |
public: |
286 | 254 |
|
287 | 255 |
// Constructor |
288 | 256 |
FirstEligiblePivotRule(NetworkSimplex &ns) : |
289 | 257 |
_source(ns._source), _target(ns._target), |
290 | 258 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
291 |
_in_arc(ns.in_arc), |
|
259 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
|
260 |
_next_arc(0) |
|
292 | 261 |
{} |
293 | 262 |
|
294 | 263 |
// Find next entering arc |
295 | 264 |
bool findEnteringArc() { |
296 | 265 |
Cost c; |
297 |
for (int e = _next_arc; e < |
|
266 |
for (int e = _next_arc; e < _search_arc_num; ++e) { |
|
298 | 267 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
299 | 268 |
if (c < 0) { |
300 | 269 |
_in_arc = e; |
301 | 270 |
_next_arc = e + 1; |
302 | 271 |
return true; |
303 | 272 |
} |
304 | 273 |
} |
305 | 274 |
for (int e = 0; e < _next_arc; ++e) { |
306 | 275 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
307 | 276 |
if (c < 0) { |
308 | 277 |
_in_arc = e; |
309 | 278 |
_next_arc = e + 1; |
310 | 279 |
return true; |
311 | 280 |
} |
312 | 281 |
} |
313 | 282 |
return false; |
314 | 283 |
} |
315 | 284 |
|
316 | 285 |
}; //class FirstEligiblePivotRule |
317 | 286 |
|
318 | 287 |
|
319 | 288 |
// Implementation of the Best Eligible pivot rule |
320 | 289 |
class BestEligiblePivotRule |
321 | 290 |
{ |
322 | 291 |
private: |
323 | 292 |
|
324 | 293 |
// References to the NetworkSimplex class |
325 | 294 |
const IntVector &_source; |
326 | 295 |
const IntVector &_target; |
327 | 296 |
const CostVector &_cost; |
328 | 297 |
const IntVector &_state; |
329 | 298 |
const CostVector &_pi; |
330 | 299 |
int &_in_arc; |
331 |
int |
|
300 |
int _search_arc_num; |
|
332 | 301 |
|
333 | 302 |
public: |
334 | 303 |
|
335 | 304 |
// Constructor |
336 | 305 |
BestEligiblePivotRule(NetworkSimplex &ns) : |
337 | 306 |
_source(ns._source), _target(ns._target), |
338 | 307 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
339 |
_in_arc(ns.in_arc), |
|
308 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) |
|
340 | 309 |
{} |
341 | 310 |
|
342 | 311 |
// Find next entering arc |
343 | 312 |
bool findEnteringArc() { |
344 | 313 |
Cost c, min = 0; |
345 |
for (int e = 0; e < |
|
314 |
for (int e = 0; e < _search_arc_num; ++e) { |
|
346 | 315 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
347 | 316 |
if (c < min) { |
348 | 317 |
min = c; |
349 | 318 |
_in_arc = e; |
350 | 319 |
} |
351 | 320 |
} |
352 | 321 |
return min < 0; |
353 | 322 |
} |
354 | 323 |
|
355 | 324 |
}; //class BestEligiblePivotRule |
356 | 325 |
|
357 | 326 |
|
358 | 327 |
// Implementation of the Block Search pivot rule |
359 | 328 |
class BlockSearchPivotRule |
360 | 329 |
{ |
361 | 330 |
private: |
362 | 331 |
|
363 | 332 |
// References to the NetworkSimplex class |
364 | 333 |
const IntVector &_source; |
365 | 334 |
const IntVector &_target; |
366 | 335 |
const CostVector &_cost; |
367 | 336 |
const IntVector &_state; |
368 | 337 |
const CostVector &_pi; |
369 | 338 |
int &_in_arc; |
370 |
int |
|
339 |
int _search_arc_num; |
|
371 | 340 |
|
372 | 341 |
// Pivot rule data |
373 | 342 |
int _block_size; |
374 | 343 |
int _next_arc; |
375 | 344 |
|
376 | 345 |
public: |
377 | 346 |
|
378 | 347 |
// Constructor |
379 | 348 |
BlockSearchPivotRule(NetworkSimplex &ns) : |
380 | 349 |
_source(ns._source), _target(ns._target), |
381 | 350 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
382 |
_in_arc(ns.in_arc), |
|
351 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
|
352 |
_next_arc(0) |
|
383 | 353 |
{ |
384 | 354 |
// The main parameters of the pivot rule |
385 |
const double BLOCK_SIZE_FACTOR = |
|
355 |
const double BLOCK_SIZE_FACTOR = 0.5; |
|
386 | 356 |
const int MIN_BLOCK_SIZE = 10; |
387 | 357 |
|
388 | 358 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * |
389 |
std::sqrt(double( |
|
359 |
std::sqrt(double(_search_arc_num))), |
|
390 | 360 |
MIN_BLOCK_SIZE ); |
391 | 361 |
} |
392 | 362 |
|
393 | 363 |
// Find next entering arc |
394 | 364 |
bool findEnteringArc() { |
395 | 365 |
Cost c, min = 0; |
396 | 366 |
int cnt = _block_size; |
397 | 367 |
int e, min_arc = _next_arc; |
398 |
for (e = _next_arc; e < |
|
368 |
for (e = _next_arc; e < _search_arc_num; ++e) { |
|
399 | 369 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
400 | 370 |
if (c < min) { |
401 | 371 |
min = c; |
402 | 372 |
min_arc = e; |
403 | 373 |
} |
404 | 374 |
if (--cnt == 0) { |
405 | 375 |
if (min < 0) break; |
406 | 376 |
cnt = _block_size; |
407 | 377 |
} |
408 | 378 |
} |
409 | 379 |
if (min == 0 || cnt > 0) { |
410 | 380 |
for (e = 0; e < _next_arc; ++e) { |
411 | 381 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
412 | 382 |
if (c < min) { |
413 | 383 |
min = c; |
414 | 384 |
min_arc = e; |
415 | 385 |
} |
416 | 386 |
if (--cnt == 0) { |
417 | 387 |
if (min < 0) break; |
418 | 388 |
cnt = _block_size; |
419 | 389 |
} |
420 | 390 |
} |
421 | 391 |
} |
422 | 392 |
if (min >= 0) return false; |
423 | 393 |
_in_arc = min_arc; |
424 | 394 |
_next_arc = e; |
425 | 395 |
return true; |
426 | 396 |
} |
427 | 397 |
|
428 | 398 |
}; //class BlockSearchPivotRule |
429 | 399 |
|
430 | 400 |
|
431 | 401 |
// Implementation of the Candidate List pivot rule |
432 | 402 |
class CandidateListPivotRule |
433 | 403 |
{ |
434 | 404 |
private: |
435 | 405 |
|
436 | 406 |
// References to the NetworkSimplex class |
437 | 407 |
const IntVector &_source; |
438 | 408 |
const IntVector &_target; |
439 | 409 |
const CostVector &_cost; |
440 | 410 |
const IntVector &_state; |
441 | 411 |
const CostVector &_pi; |
442 | 412 |
int &_in_arc; |
443 |
int |
|
413 |
int _search_arc_num; |
|
444 | 414 |
|
445 | 415 |
// Pivot rule data |
446 | 416 |
IntVector _candidates; |
447 | 417 |
int _list_length, _minor_limit; |
448 | 418 |
int _curr_length, _minor_count; |
449 | 419 |
int _next_arc; |
450 | 420 |
|
451 | 421 |
public: |
452 | 422 |
|
453 | 423 |
/// Constructor |
454 | 424 |
CandidateListPivotRule(NetworkSimplex &ns) : |
455 | 425 |
_source(ns._source), _target(ns._target), |
456 | 426 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
457 |
_in_arc(ns.in_arc), |
|
427 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
|
428 |
_next_arc(0) |
|
458 | 429 |
{ |
459 | 430 |
// The main parameters of the pivot rule |
460 | 431 |
const double LIST_LENGTH_FACTOR = 1.0; |
461 | 432 |
const int MIN_LIST_LENGTH = 10; |
462 | 433 |
const double MINOR_LIMIT_FACTOR = 0.1; |
463 | 434 |
const int MIN_MINOR_LIMIT = 3; |
464 | 435 |
|
465 | 436 |
_list_length = std::max( int(LIST_LENGTH_FACTOR * |
466 |
std::sqrt(double( |
|
437 |
std::sqrt(double(_search_arc_num))), |
|
467 | 438 |
MIN_LIST_LENGTH ); |
468 | 439 |
_minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), |
469 | 440 |
MIN_MINOR_LIMIT ); |
470 | 441 |
_curr_length = _minor_count = 0; |
471 | 442 |
_candidates.resize(_list_length); |
472 | 443 |
} |
473 | 444 |
|
474 | 445 |
/// Find next entering arc |
475 | 446 |
bool findEnteringArc() { |
476 | 447 |
Cost min, c; |
477 | 448 |
int e, min_arc = _next_arc; |
478 | 449 |
if (_curr_length > 0 && _minor_count < _minor_limit) { |
479 | 450 |
// Minor iteration: select the best eligible arc from the |
480 | 451 |
// current candidate list |
481 | 452 |
++_minor_count; |
482 | 453 |
min = 0; |
483 | 454 |
for (int i = 0; i < _curr_length; ++i) { |
484 | 455 |
e = _candidates[i]; |
485 | 456 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
486 | 457 |
if (c < min) { |
487 | 458 |
min = c; |
488 | 459 |
min_arc = e; |
489 | 460 |
} |
490 | 461 |
if (c >= 0) { |
491 | 462 |
_candidates[i--] = _candidates[--_curr_length]; |
492 | 463 |
} |
493 | 464 |
} |
494 | 465 |
if (min < 0) { |
495 | 466 |
_in_arc = min_arc; |
496 | 467 |
return true; |
497 | 468 |
} |
498 | 469 |
} |
499 | 470 |
|
500 | 471 |
// Major iteration: build a new candidate list |
501 | 472 |
min = 0; |
502 | 473 |
_curr_length = 0; |
503 |
for (e = _next_arc; e < |
|
474 |
for (e = _next_arc; e < _search_arc_num; ++e) { |
|
504 | 475 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
505 | 476 |
if (c < 0) { |
506 | 477 |
_candidates[_curr_length++] = e; |
507 | 478 |
if (c < min) { |
508 | 479 |
min = c; |
509 | 480 |
min_arc = e; |
510 | 481 |
} |
511 | 482 |
if (_curr_length == _list_length) break; |
512 | 483 |
} |
513 | 484 |
} |
514 | 485 |
if (_curr_length < _list_length) { |
515 | 486 |
for (e = 0; e < _next_arc; ++e) { |
516 | 487 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
517 | 488 |
if (c < 0) { |
518 | 489 |
_candidates[_curr_length++] = e; |
519 | 490 |
if (c < min) { |
520 | 491 |
min = c; |
521 | 492 |
min_arc = e; |
522 | 493 |
} |
523 | 494 |
if (_curr_length == _list_length) break; |
524 | 495 |
} |
525 | 496 |
} |
526 | 497 |
} |
527 | 498 |
if (_curr_length == 0) return false; |
528 | 499 |
_minor_count = 1; |
529 | 500 |
_in_arc = min_arc; |
530 | 501 |
_next_arc = e; |
531 | 502 |
return true; |
532 | 503 |
} |
533 | 504 |
|
534 | 505 |
}; //class CandidateListPivotRule |
535 | 506 |
|
536 | 507 |
|
537 | 508 |
// Implementation of the Altering Candidate List pivot rule |
538 | 509 |
class AlteringListPivotRule |
539 | 510 |
{ |
540 | 511 |
private: |
541 | 512 |
|
542 | 513 |
// References to the NetworkSimplex class |
543 | 514 |
const IntVector &_source; |
544 | 515 |
const IntVector &_target; |
545 | 516 |
const CostVector &_cost; |
546 | 517 |
const IntVector &_state; |
547 | 518 |
const CostVector &_pi; |
548 | 519 |
int &_in_arc; |
549 |
int |
|
520 |
int _search_arc_num; |
|
550 | 521 |
|
551 | 522 |
// Pivot rule data |
552 | 523 |
int _block_size, _head_length, _curr_length; |
553 | 524 |
int _next_arc; |
554 | 525 |
IntVector _candidates; |
555 | 526 |
CostVector _cand_cost; |
556 | 527 |
|
557 | 528 |
// Functor class to compare arcs during sort of the candidate list |
558 | 529 |
class SortFunc |
559 | 530 |
{ |
560 | 531 |
private: |
561 | 532 |
const CostVector &_map; |
562 | 533 |
public: |
563 | 534 |
SortFunc(const CostVector &map) : _map(map) {} |
564 | 535 |
bool operator()(int left, int right) { |
565 | 536 |
return _map[left] > _map[right]; |
566 | 537 |
} |
567 | 538 |
}; |
568 | 539 |
|
569 | 540 |
SortFunc _sort_func; |
570 | 541 |
|
571 | 542 |
public: |
572 | 543 |
|
573 | 544 |
// Constructor |
574 | 545 |
AlteringListPivotRule(NetworkSimplex &ns) : |
575 | 546 |
_source(ns._source), _target(ns._target), |
576 | 547 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
577 |
_in_arc(ns.in_arc), _arc_num(ns._arc_num), |
|
578 |
_next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost) |
|
548 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
|
549 |
_next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) |
|
579 | 550 |
{ |
580 | 551 |
// The main parameters of the pivot rule |
581 | 552 |
const double BLOCK_SIZE_FACTOR = 1.5; |
582 | 553 |
const int MIN_BLOCK_SIZE = 10; |
583 | 554 |
const double HEAD_LENGTH_FACTOR = 0.1; |
584 | 555 |
const int MIN_HEAD_LENGTH = 3; |
585 | 556 |
|
586 | 557 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * |
587 |
std::sqrt(double( |
|
558 |
std::sqrt(double(_search_arc_num))), |
|
588 | 559 |
MIN_BLOCK_SIZE ); |
589 | 560 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), |
590 | 561 |
MIN_HEAD_LENGTH ); |
591 | 562 |
_candidates.resize(_head_length + _block_size); |
592 | 563 |
_curr_length = 0; |
593 | 564 |
} |
594 | 565 |
|
595 | 566 |
// Find next entering arc |
596 | 567 |
bool findEnteringArc() { |
597 | 568 |
// Check the current candidate list |
598 | 569 |
int e; |
599 | 570 |
for (int i = 0; i < _curr_length; ++i) { |
600 | 571 |
e = _candidates[i]; |
601 | 572 |
_cand_cost[e] = _state[e] * |
602 | 573 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
603 | 574 |
if (_cand_cost[e] >= 0) { |
604 | 575 |
_candidates[i--] = _candidates[--_curr_length]; |
605 | 576 |
} |
606 | 577 |
} |
607 | 578 |
|
608 | 579 |
// Extend the list |
609 | 580 |
int cnt = _block_size; |
610 | 581 |
int last_arc = 0; |
611 | 582 |
int limit = _head_length; |
612 | 583 |
|
613 |
for (int e = _next_arc; e < |
|
584 |
for (int e = _next_arc; e < _search_arc_num; ++e) { |
|
614 | 585 |
_cand_cost[e] = _state[e] * |
615 | 586 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
616 | 587 |
if (_cand_cost[e] < 0) { |
617 | 588 |
_candidates[_curr_length++] = e; |
618 | 589 |
last_arc = e; |
619 | 590 |
} |
620 | 591 |
if (--cnt == 0) { |
621 | 592 |
if (_curr_length > limit) break; |
622 | 593 |
limit = 0; |
623 | 594 |
cnt = _block_size; |
624 | 595 |
} |
625 | 596 |
} |
626 | 597 |
if (_curr_length <= limit) { |
627 | 598 |
for (int e = 0; e < _next_arc; ++e) { |
628 | 599 |
_cand_cost[e] = _state[e] * |
629 | 600 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
630 | 601 |
if (_cand_cost[e] < 0) { |
631 | 602 |
_candidates[_curr_length++] = e; |
632 | 603 |
last_arc = e; |
633 | 604 |
} |
634 | 605 |
if (--cnt == 0) { |
635 | 606 |
if (_curr_length > limit) break; |
636 | 607 |
limit = 0; |
637 | 608 |
cnt = _block_size; |
638 | 609 |
} |
639 | 610 |
} |
640 | 611 |
} |
641 | 612 |
if (_curr_length == 0) return false; |
642 | 613 |
_next_arc = last_arc + 1; |
643 | 614 |
|
644 | 615 |
// Make heap of the candidate list (approximating a partial sort) |
645 | 616 |
make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
646 | 617 |
_sort_func ); |
647 | 618 |
|
648 | 619 |
// Pop the first element of the heap |
649 | 620 |
_in_arc = _candidates[0]; |
650 | 621 |
pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
651 | 622 |
_sort_func ); |
652 | 623 |
_curr_length = std::min(_head_length, _curr_length - 1); |
653 | 624 |
return true; |
654 | 625 |
} |
655 | 626 |
|
656 | 627 |
}; //class AlteringListPivotRule |
657 | 628 |
|
658 | 629 |
public: |
659 | 630 |
|
660 | 631 |
/// \brief Constructor. |
661 | 632 |
/// |
662 | 633 |
/// The constructor of the class. |
663 | 634 |
/// |
664 | 635 |
/// \param graph The digraph the algorithm runs on. |
665 | 636 |
NetworkSimplex(const GR& graph) : |
666 | 637 |
_graph(graph), _node_id(graph), _arc_id(graph), |
667 | 638 |
INF(std::numeric_limits<Value>::has_infinity ? |
668 | 639 |
std::numeric_limits<Value>::infinity() : |
669 | 640 |
std::numeric_limits<Value>::max()) |
670 | 641 |
{ |
671 | 642 |
// Check the value types |
672 | 643 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
673 | 644 |
"The flow type of NetworkSimplex must be signed"); |
674 | 645 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
675 | 646 |
"The cost type of NetworkSimplex must be signed"); |
676 | 647 |
|
677 | 648 |
// Resize vectors |
678 | 649 |
_node_num = countNodes(_graph); |
679 | 650 |
_arc_num = countArcs(_graph); |
680 | 651 |
int all_node_num = _node_num + 1; |
681 |
int |
|
652 |
int max_arc_num = _arc_num + 2 * _node_num; |
|
682 | 653 |
|
683 |
_source.resize(all_arc_num); |
|
684 |
_target.resize(all_arc_num); |
|
654 |
_source.resize(max_arc_num); |
|
655 |
_target.resize(max_arc_num); |
|
685 | 656 |
|
686 |
_lower.resize(all_arc_num); |
|
687 |
_upper.resize(all_arc_num); |
|
688 |
_cap.resize(all_arc_num); |
|
689 |
_cost.resize(all_arc_num); |
|
657 |
_lower.resize(_arc_num); |
|
658 |
_upper.resize(_arc_num); |
|
659 |
_cap.resize(max_arc_num); |
|
660 |
_cost.resize(max_arc_num); |
|
690 | 661 |
_supply.resize(all_node_num); |
691 |
_flow.resize( |
|
662 |
_flow.resize(max_arc_num); |
|
692 | 663 |
_pi.resize(all_node_num); |
693 | 664 |
|
694 | 665 |
_parent.resize(all_node_num); |
695 | 666 |
_pred.resize(all_node_num); |
696 | 667 |
_forward.resize(all_node_num); |
697 | 668 |
_thread.resize(all_node_num); |
698 | 669 |
_rev_thread.resize(all_node_num); |
699 | 670 |
_succ_num.resize(all_node_num); |
700 | 671 |
_last_succ.resize(all_node_num); |
701 |
_state.resize( |
|
672 |
_state.resize(max_arc_num); |
|
702 | 673 |
|
703 | 674 |
// Copy the graph (store the arcs in a mixed order) |
704 | 675 |
int i = 0; |
705 | 676 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
706 | 677 |
_node_id[n] = i; |
707 | 678 |
} |
708 | 679 |
int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
709 | 680 |
i = 0; |
710 | 681 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
711 | 682 |
_arc_id[a] = i; |
712 | 683 |
_source[i] = _node_id[_graph.source(a)]; |
713 | 684 |
_target[i] = _node_id[_graph.target(a)]; |
714 | 685 |
if ((i += k) >= _arc_num) i = (i % k) + 1; |
715 | 686 |
} |
716 | 687 |
|
717 | 688 |
// Initialize maps |
718 | 689 |
for (int i = 0; i != _node_num; ++i) { |
719 | 690 |
_supply[i] = 0; |
720 | 691 |
} |
721 | 692 |
for (int i = 0; i != _arc_num; ++i) { |
722 | 693 |
_lower[i] = 0; |
723 | 694 |
_upper[i] = INF; |
724 | 695 |
_cost[i] = 1; |
725 | 696 |
} |
726 | 697 |
_have_lower = false; |
727 | 698 |
_stype = GEQ; |
728 | 699 |
} |
729 | 700 |
|
730 | 701 |
/// \name Parameters |
731 | 702 |
/// The parameters of the algorithm can be specified using these |
732 | 703 |
/// functions. |
733 | 704 |
|
734 | 705 |
/// @{ |
735 | 706 |
|
736 | 707 |
/// \brief Set the lower bounds on the arcs. |
737 | 708 |
/// |
738 | 709 |
/// This function sets the lower bounds on the arcs. |
739 | 710 |
/// If it is not used before calling \ref run(), the lower bounds |
740 | 711 |
/// will be set to zero on all arcs. |
741 | 712 |
/// |
742 | 713 |
/// \param map An arc map storing the lower bounds. |
743 | 714 |
/// Its \c Value type must be convertible to the \c Value type |
744 | 715 |
/// of the algorithm. |
745 | 716 |
/// |
746 | 717 |
/// \return <tt>(*this)</tt> |
747 | 718 |
template <typename LowerMap> |
748 | 719 |
NetworkSimplex& lowerMap(const LowerMap& map) { |
749 | 720 |
_have_lower = true; |
750 | 721 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
751 | 722 |
_lower[_arc_id[a]] = map[a]; |
752 | 723 |
} |
753 | 724 |
return *this; |
754 | 725 |
} |
755 | 726 |
|
756 | 727 |
/// \brief Set the upper bounds (capacities) on the arcs. |
757 | 728 |
/// |
758 | 729 |
/// This function sets the upper bounds (capacities) on the arcs. |
759 | 730 |
/// If it is not used before calling \ref run(), the upper bounds |
760 | 731 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
761 | 732 |
/// unbounded from above on each arc). |
762 | 733 |
/// |
763 | 734 |
/// \param map An arc map storing the upper bounds. |
764 | 735 |
/// Its \c Value type must be convertible to the \c Value type |
765 | 736 |
/// of the algorithm. |
766 | 737 |
/// |
767 | 738 |
/// \return <tt>(*this)</tt> |
768 | 739 |
template<typename UpperMap> |
769 | 740 |
NetworkSimplex& upperMap(const UpperMap& map) { |
770 | 741 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
771 | 742 |
_upper[_arc_id[a]] = map[a]; |
772 | 743 |
} |
773 | 744 |
return *this; |
774 | 745 |
} |
775 | 746 |
|
776 | 747 |
/// \brief Set the costs of the arcs. |
777 | 748 |
/// |
778 | 749 |
/// This function sets the costs of the arcs. |
779 | 750 |
/// If it is not used before calling \ref run(), the costs |
780 | 751 |
/// will be set to \c 1 on all arcs. |
781 | 752 |
/// |
782 | 753 |
/// \param map An arc map storing the costs. |
783 | 754 |
/// Its \c Value type must be convertible to the \c Cost type |
784 | 755 |
/// of the algorithm. |
785 | 756 |
/// |
786 | 757 |
/// \return <tt>(*this)</tt> |
787 | 758 |
template<typename CostMap> |
788 | 759 |
NetworkSimplex& costMap(const CostMap& map) { |
789 | 760 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
790 | 761 |
_cost[_arc_id[a]] = map[a]; |
791 | 762 |
} |
792 | 763 |
return *this; |
793 | 764 |
} |
794 | 765 |
|
795 | 766 |
/// \brief Set the supply values of the nodes. |
796 | 767 |
/// |
797 | 768 |
/// This function sets the supply values of the nodes. |
798 | 769 |
/// If neither this function nor \ref stSupply() is used before |
799 | 770 |
/// calling \ref run(), the supply of each node will be set to zero. |
800 | 771 |
/// (It makes sense only if non-zero lower bounds are given.) |
801 | 772 |
/// |
802 | 773 |
/// \param map A node map storing the supply values. |
803 | 774 |
/// Its \c Value type must be convertible to the \c Value type |
804 | 775 |
/// of the algorithm. |
805 | 776 |
/// |
806 | 777 |
/// \return <tt>(*this)</tt> |
807 | 778 |
template<typename SupplyMap> |
808 | 779 |
NetworkSimplex& supplyMap(const SupplyMap& map) { |
809 | 780 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
810 | 781 |
_supply[_node_id[n]] = map[n]; |
811 | 782 |
} |
812 | 783 |
return *this; |
813 | 784 |
} |
814 | 785 |
|
815 | 786 |
/// \brief Set single source and target nodes and a supply value. |
816 | 787 |
/// |
817 | 788 |
/// This function sets a single source node and a single target node |
818 | 789 |
/// and the required flow value. |
819 | 790 |
/// If neither this function nor \ref supplyMap() is used before |
820 | 791 |
/// calling \ref run(), the supply of each node will be set to zero. |
821 | 792 |
/// (It makes sense only if non-zero lower bounds are given.) |
822 | 793 |
/// |
823 | 794 |
/// Using this function has the same effect as using \ref supplyMap() |
824 | 795 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
825 | 796 |
/// assigned to \c t and all other nodes have zero supply value. |
826 | 797 |
/// |
827 | 798 |
/// \param s The source node. |
828 | 799 |
/// \param t The target node. |
829 | 800 |
/// \param k The required amount of flow from node \c s to node \c t |
830 | 801 |
/// (i.e. the supply of \c s and the demand of \c t). |
831 | 802 |
/// |
832 | 803 |
/// \return <tt>(*this)</tt> |
833 | 804 |
NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { |
834 | 805 |
for (int i = 0; i != _node_num; ++i) { |
835 | 806 |
_supply[i] = 0; |
836 | 807 |
} |
837 | 808 |
_supply[_node_id[s]] = k; |
838 | 809 |
_supply[_node_id[t]] = -k; |
839 | 810 |
return *this; |
840 | 811 |
} |
841 | 812 |
|
842 | 813 |
/// \brief Set the type of the supply constraints. |
843 | 814 |
/// |
844 | 815 |
/// This function sets the type of the supply/demand constraints. |
845 | 816 |
/// If it is not used before calling \ref run(), the \ref GEQ supply |
846 | 817 |
/// type will be used. |
847 | 818 |
/// |
848 | 819 |
/// For more information see \ref SupplyType. |
849 | 820 |
/// |
850 | 821 |
/// \return <tt>(*this)</tt> |
851 | 822 |
NetworkSimplex& supplyType(SupplyType supply_type) { |
852 | 823 |
_stype = supply_type; |
853 | 824 |
return *this; |
854 | 825 |
} |
855 | 826 |
|
856 | 827 |
/// @} |
857 | 828 |
|
858 | 829 |
/// \name Execution Control |
859 | 830 |
/// The algorithm can be executed using \ref run(). |
860 | 831 |
|
861 | 832 |
/// @{ |
862 | 833 |
|
863 | 834 |
/// \brief Run the algorithm. |
864 | 835 |
/// |
865 | 836 |
/// This function runs the algorithm. |
866 | 837 |
/// The paramters can be specified using functions \ref lowerMap(), |
867 | 838 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), |
868 | 839 |
/// \ref supplyType(). |
869 | 840 |
/// For example, |
870 | 841 |
/// \code |
871 | 842 |
/// NetworkSimplex<ListDigraph> ns(graph); |
872 | 843 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
873 | 844 |
/// .supplyMap(sup).run(); |
874 | 845 |
/// \endcode |
875 | 846 |
/// |
876 | 847 |
/// This function can be called more than once. All the parameters |
877 | 848 |
/// that have been given are kept for the next call, unless |
878 | 849 |
/// \ref reset() is called, thus only the modified parameters |
879 | 850 |
/// have to be set again. See \ref reset() for examples. |
880 | 851 |
/// However the underlying digraph must not be modified after this |
881 | 852 |
/// class have been constructed, since it copies and extends the graph. |
882 | 853 |
/// |
883 | 854 |
/// \param pivot_rule The pivot rule that will be used during the |
884 | 855 |
/// algorithm. For more information see \ref PivotRule. |
885 | 856 |
/// |
886 | 857 |
/// \return \c INFEASIBLE if no feasible flow exists, |
887 | 858 |
/// \n \c OPTIMAL if the problem has optimal solution |
888 | 859 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
889 | 860 |
/// optimal flow and node potentials (primal and dual solutions), |
890 | 861 |
/// \n \c UNBOUNDED if the objective function of the problem is |
891 | 862 |
/// unbounded, i.e. there is a directed cycle having negative total |
892 | 863 |
/// cost and infinite upper bound. |
893 | 864 |
/// |
894 | 865 |
/// \see ProblemType, PivotRule |
895 | 866 |
ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { |
896 | 867 |
if (!init()) return INFEASIBLE; |
897 | 868 |
return start(pivot_rule); |
898 | 869 |
} |
899 | 870 |
|
900 | 871 |
/// \brief Reset all the parameters that have been given before. |
901 | 872 |
/// |
902 | 873 |
/// This function resets all the paramaters that have been given |
903 | 874 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
904 | 875 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). |
905 | 876 |
/// |
906 | 877 |
/// It is useful for multiple run() calls. If this function is not |
907 | 878 |
/// used, all the parameters given before are kept for the next |
908 | 879 |
/// \ref run() call. |
909 | 880 |
/// However the underlying digraph must not be modified after this |
910 | 881 |
/// class have been constructed, since it copies and extends the graph. |
911 | 882 |
/// |
912 | 883 |
/// For example, |
913 | 884 |
/// \code |
914 | 885 |
/// NetworkSimplex<ListDigraph> ns(graph); |
915 | 886 |
/// |
916 | 887 |
/// // First run |
917 | 888 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
918 | 889 |
/// .supplyMap(sup).run(); |
919 | 890 |
/// |
920 | 891 |
/// // Run again with modified cost map (reset() is not called, |
921 | 892 |
/// // so only the cost map have to be set again) |
922 | 893 |
/// cost[e] += 100; |
923 | 894 |
/// ns.costMap(cost).run(); |
924 | 895 |
/// |
925 | 896 |
/// // Run again from scratch using reset() |
926 | 897 |
/// // (the lower bounds will be set to zero on all arcs) |
927 | 898 |
/// ns.reset(); |
928 | 899 |
/// ns.upperMap(capacity).costMap(cost) |
929 | 900 |
/// .supplyMap(sup).run(); |
930 | 901 |
/// \endcode |
931 | 902 |
/// |
932 | 903 |
/// \return <tt>(*this)</tt> |
933 | 904 |
NetworkSimplex& reset() { |
934 | 905 |
for (int i = 0; i != _node_num; ++i) { |
935 | 906 |
_supply[i] = 0; |
936 | 907 |
} |
937 | 908 |
for (int i = 0; i != _arc_num; ++i) { |
938 | 909 |
_lower[i] = 0; |
939 | 910 |
_upper[i] = INF; |
940 | 911 |
_cost[i] = 1; |
941 | 912 |
} |
942 | 913 |
_have_lower = false; |
943 | 914 |
_stype = GEQ; |
944 | 915 |
return *this; |
945 | 916 |
} |
946 | 917 |
|
947 | 918 |
/// @} |
948 | 919 |
|
949 | 920 |
/// \name Query Functions |
950 | 921 |
/// The results of the algorithm can be obtained using these |
951 | 922 |
/// functions.\n |
952 | 923 |
/// The \ref run() function must be called before using them. |
953 | 924 |
|
954 | 925 |
/// @{ |
955 | 926 |
|
956 | 927 |
/// \brief Return the total cost of the found flow. |
957 | 928 |
/// |
958 | 929 |
/// This function returns the total cost of the found flow. |
959 | 930 |
/// Its complexity is O(e). |
960 | 931 |
/// |
961 | 932 |
/// \note The return type of the function can be specified as a |
962 | 933 |
/// template parameter. For example, |
963 | 934 |
/// \code |
964 | 935 |
/// ns.totalCost<double>(); |
965 | 936 |
/// \endcode |
966 | 937 |
/// It is useful if the total cost cannot be stored in the \c Cost |
967 | 938 |
/// type of the algorithm, which is the default return type of the |
968 | 939 |
/// function. |
969 | 940 |
/// |
970 | 941 |
/// \pre \ref run() must be called before using this function. |
971 | 942 |
template <typename Number> |
972 | 943 |
Number totalCost() const { |
973 | 944 |
Number c = 0; |
974 | 945 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
975 | 946 |
int i = _arc_id[a]; |
976 | 947 |
c += Number(_flow[i]) * Number(_cost[i]); |
977 | 948 |
} |
978 | 949 |
return c; |
979 | 950 |
} |
980 | 951 |
|
981 | 952 |
#ifndef DOXYGEN |
982 | 953 |
Cost totalCost() const { |
983 | 954 |
return totalCost<Cost>(); |
984 | 955 |
} |
985 | 956 |
#endif |
986 | 957 |
|
987 | 958 |
/// \brief Return the flow on the given arc. |
988 | 959 |
/// |
989 | 960 |
/// This function returns the flow on the given arc. |
990 | 961 |
/// |
991 | 962 |
/// \pre \ref run() must be called before using this function. |
992 | 963 |
Value flow(const Arc& a) const { |
993 | 964 |
return _flow[_arc_id[a]]; |
994 | 965 |
} |
995 | 966 |
|
996 | 967 |
/// \brief Return the flow map (the primal solution). |
997 | 968 |
/// |
998 | 969 |
/// This function copies the flow value on each arc into the given |
999 | 970 |
/// map. The \c Value type of the algorithm must be convertible to |
1000 | 971 |
/// the \c Value type of the map. |
1001 | 972 |
/// |
1002 | 973 |
/// \pre \ref run() must be called before using this function. |
1003 | 974 |
template <typename FlowMap> |
1004 | 975 |
void flowMap(FlowMap &map) const { |
1005 | 976 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
1006 | 977 |
map.set(a, _flow[_arc_id[a]]); |
1007 | 978 |
} |
1008 | 979 |
} |
1009 | 980 |
|
1010 | 981 |
/// \brief Return the potential (dual value) of the given node. |
1011 | 982 |
/// |
1012 | 983 |
/// This function returns the potential (dual value) of the |
1013 | 984 |
/// given node. |
1014 | 985 |
/// |
1015 | 986 |
/// \pre \ref run() must be called before using this function. |
1016 | 987 |
Cost potential(const Node& n) const { |
1017 | 988 |
return _pi[_node_id[n]]; |
1018 | 989 |
} |
1019 | 990 |
|
1020 | 991 |
/// \brief Return the potential map (the dual solution). |
1021 | 992 |
/// |
1022 | 993 |
/// This function copies the potential (dual value) of each node |
1023 | 994 |
/// into the given map. |
1024 | 995 |
/// The \c Cost type of the algorithm must be convertible to the |
1025 | 996 |
/// \c Value type of the map. |
1026 | 997 |
/// |
1027 | 998 |
/// \pre \ref run() must be called before using this function. |
1028 | 999 |
template <typename PotentialMap> |
1029 | 1000 |
void potentialMap(PotentialMap &map) const { |
1030 | 1001 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
1031 | 1002 |
map.set(n, _pi[_node_id[n]]); |
1032 | 1003 |
} |
1033 | 1004 |
} |
1034 | 1005 |
|
1035 | 1006 |
/// @} |
1036 | 1007 |
|
1037 | 1008 |
private: |
1038 | 1009 |
|
1039 | 1010 |
// Initialize internal data structures |
1040 | 1011 |
bool init() { |
1041 | 1012 |
if (_node_num == 0) return false; |
1042 | 1013 |
|
1043 | 1014 |
// Check the sum of supply values |
1044 | 1015 |
_sum_supply = 0; |
1045 | 1016 |
for (int i = 0; i != _node_num; ++i) { |
1046 | 1017 |
_sum_supply += _supply[i]; |
1047 | 1018 |
} |
1048 | 1019 |
if ( !((_stype == GEQ && _sum_supply <= 0) || |
1049 | 1020 |
(_stype == LEQ && _sum_supply >= 0)) ) return false; |
1050 | 1021 |
|
1051 | 1022 |
// Remove non-zero lower bounds |
1052 | 1023 |
if (_have_lower) { |
1053 | 1024 |
for (int i = 0; i != _arc_num; ++i) { |
1054 | 1025 |
Value c = _lower[i]; |
1055 | 1026 |
if (c >= 0) { |
1056 | 1027 |
_cap[i] = _upper[i] < INF ? _upper[i] - c : INF; |
1057 | 1028 |
} else { |
1058 | 1029 |
_cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF; |
1059 | 1030 |
} |
1060 | 1031 |
_supply[_source[i]] -= c; |
1061 | 1032 |
_supply[_target[i]] += c; |
1062 | 1033 |
} |
1063 | 1034 |
} else { |
1064 | 1035 |
for (int i = 0; i != _arc_num; ++i) { |
1065 | 1036 |
_cap[i] = _upper[i]; |
1066 | 1037 |
} |
1067 | 1038 |
} |
1068 | 1039 |
|
1069 | 1040 |
// Initialize artifical cost |
1070 | 1041 |
Cost ART_COST; |
1071 | 1042 |
if (std::numeric_limits<Cost>::is_exact) { |
1072 |
ART_COST = std::numeric_limits<Cost>::max() / |
|
1043 |
ART_COST = std::numeric_limits<Cost>::max() / 2 + 1; |
|
1073 | 1044 |
} else { |
1074 | 1045 |
ART_COST = std::numeric_limits<Cost>::min(); |
1075 | 1046 |
for (int i = 0; i != _arc_num; ++i) { |
1076 | 1047 |
if (_cost[i] > ART_COST) ART_COST = _cost[i]; |
1077 | 1048 |
} |
1078 | 1049 |
ART_COST = (ART_COST + 1) * _node_num; |
1079 | 1050 |
} |
1080 | 1051 |
|
1081 | 1052 |
// Initialize arc maps |
1082 | 1053 |
for (int i = 0; i != _arc_num; ++i) { |
1083 | 1054 |
_flow[i] = 0; |
1084 | 1055 |
_state[i] = STATE_LOWER; |
1085 | 1056 |
} |
1086 | 1057 |
|
1087 | 1058 |
// Set data for the artificial root node |
1088 | 1059 |
_root = _node_num; |
1089 | 1060 |
_parent[_root] = -1; |
1090 | 1061 |
_pred[_root] = -1; |
1091 | 1062 |
_thread[_root] = 0; |
1092 | 1063 |
_rev_thread[0] = _root; |
1093 | 1064 |
_succ_num[_root] = _node_num + 1; |
1094 | 1065 |
_last_succ[_root] = _root - 1; |
1095 | 1066 |
_supply[_root] = -_sum_supply; |
1096 |
_pi[_root] = |
|
1067 |
_pi[_root] = 0; |
|
1097 | 1068 |
|
1098 | 1069 |
// Add artificial arcs and initialize the spanning tree data structure |
1099 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1100 |
_parent[u] = _root; |
|
1101 |
_pred[u] = e; |
|
1102 |
_thread[u] = u + 1; |
|
1103 |
_rev_thread[u + 1] = u; |
|
1104 |
_succ_num[u] = 1; |
|
1105 |
_last_succ[u] = u; |
|
1106 |
_cost[e] = ART_COST; |
|
1107 |
_cap[e] = INF; |
|
1108 |
_state[e] = STATE_TREE; |
|
1109 |
if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) { |
|
1110 |
_flow[e] = _supply[u]; |
|
1111 |
_forward[u] = true; |
|
1112 |
_pi[u] = -ART_COST + _pi[_root]; |
|
1113 |
} else { |
|
1114 |
_flow[e] = -_supply[u]; |
|
1115 |
_forward[u] = false; |
|
1116 |
_pi[u] = ART_COST + _pi[_root]; |
|
1070 |
if (_sum_supply == 0) { |
|
1071 |
// EQ supply constraints |
|
1072 |
_search_arc_num = _arc_num; |
|
1073 |
_all_arc_num = _arc_num + _node_num; |
|
1074 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1075 |
_parent[u] = _root; |
|
1076 |
_pred[u] = e; |
|
1077 |
_thread[u] = u + 1; |
|
1078 |
_rev_thread[u + 1] = u; |
|
1079 |
_succ_num[u] = 1; |
|
1080 |
_last_succ[u] = u; |
|
1081 |
_cap[e] = INF; |
|
1082 |
_state[e] = STATE_TREE; |
|
1083 |
if (_supply[u] >= 0) { |
|
1084 |
_forward[u] = true; |
|
1085 |
_pi[u] = 0; |
|
1086 |
_source[e] = u; |
|
1087 |
_target[e] = _root; |
|
1088 |
_flow[e] = _supply[u]; |
|
1089 |
_cost[e] = 0; |
|
1090 |
} else { |
|
1091 |
_forward[u] = false; |
|
1092 |
_pi[u] = ART_COST; |
|
1093 |
_source[e] = _root; |
|
1094 |
_target[e] = u; |
|
1095 |
_flow[e] = -_supply[u]; |
|
1096 |
_cost[e] = ART_COST; |
|
1097 |
} |
|
1117 | 1098 |
} |
1118 | 1099 |
} |
1100 |
else if (_sum_supply > 0) { |
|
1101 |
// LEQ supply constraints |
|
1102 |
_search_arc_num = _arc_num + _node_num; |
|
1103 |
int f = _arc_num + _node_num; |
|
1104 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1105 |
_parent[u] = _root; |
|
1106 |
_thread[u] = u + 1; |
|
1107 |
_rev_thread[u + 1] = u; |
|
1108 |
_succ_num[u] = 1; |
|
1109 |
_last_succ[u] = u; |
|
1110 |
if (_supply[u] >= 0) { |
|
1111 |
_forward[u] = true; |
|
1112 |
_pi[u] = 0; |
|
1113 |
_pred[u] = e; |
|
1114 |
_source[e] = u; |
|
1115 |
_target[e] = _root; |
|
1116 |
_cap[e] = INF; |
|
1117 |
_flow[e] = _supply[u]; |
|
1118 |
_cost[e] = 0; |
|
1119 |
_state[e] = STATE_TREE; |
|
1120 |
} else { |
|
1121 |
_forward[u] = false; |
|
1122 |
_pi[u] = ART_COST; |
|
1123 |
_pred[u] = f; |
|
1124 |
_source[f] = _root; |
|
1125 |
_target[f] = u; |
|
1126 |
_cap[f] = INF; |
|
1127 |
_flow[f] = -_supply[u]; |
|
1128 |
_cost[f] = ART_COST; |
|
1129 |
_state[f] = STATE_TREE; |
|
1130 |
_source[e] = u; |
|
1131 |
_target[e] = _root; |
|
1132 |
_cap[e] = INF; |
|
1133 |
_flow[e] = 0; |
|
1134 |
_cost[e] = 0; |
|
1135 |
_state[e] = STATE_LOWER; |
|
1136 |
++f; |
|
1137 |
} |
|
1138 |
} |
|
1139 |
_all_arc_num = f; |
|
1140 |
} |
|
1141 |
else { |
|
1142 |
// GEQ supply constraints |
|
1143 |
_search_arc_num = _arc_num + _node_num; |
|
1144 |
int f = _arc_num + _node_num; |
|
1145 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1146 |
_parent[u] = _root; |
|
1147 |
_thread[u] = u + 1; |
|
1148 |
_rev_thread[u + 1] = u; |
|
1149 |
_succ_num[u] = 1; |
|
1150 |
_last_succ[u] = u; |
|
1151 |
if (_supply[u] <= 0) { |
|
1152 |
_forward[u] = false; |
|
1153 |
_pi[u] = 0; |
|
1154 |
_pred[u] = e; |
|
1155 |
_source[e] = _root; |
|
1156 |
_target[e] = u; |
|
1157 |
_cap[e] = INF; |
|
1158 |
_flow[e] = -_supply[u]; |
|
1159 |
_cost[e] = 0; |
|
1160 |
_state[e] = STATE_TREE; |
|
1161 |
} else { |
|
1162 |
_forward[u] = true; |
|
1163 |
_pi[u] = -ART_COST; |
|
1164 |
_pred[u] = f; |
|
1165 |
_source[f] = u; |
|
1166 |
_target[f] = _root; |
|
1167 |
_cap[f] = INF; |
|
1168 |
_flow[f] = _supply[u]; |
|
1169 |
_state[f] = STATE_TREE; |
|
1170 |
_cost[f] = ART_COST; |
|
1171 |
_source[e] = _root; |
|
1172 |
_target[e] = u; |
|
1173 |
_cap[e] = INF; |
|
1174 |
_flow[e] = 0; |
|
1175 |
_cost[e] = 0; |
|
1176 |
_state[e] = STATE_LOWER; |
|
1177 |
++f; |
|
1178 |
} |
|
1179 |
} |
|
1180 |
_all_arc_num = f; |
|
1181 |
} |
|
1119 | 1182 |
|
1120 | 1183 |
return true; |
1121 | 1184 |
} |
1122 | 1185 |
|
1123 | 1186 |
// Find the join node |
1124 | 1187 |
void findJoinNode() { |
1125 | 1188 |
int u = _source[in_arc]; |
1126 | 1189 |
int v = _target[in_arc]; |
1127 | 1190 |
while (u != v) { |
1128 | 1191 |
if (_succ_num[u] < _succ_num[v]) { |
1129 | 1192 |
u = _parent[u]; |
1130 | 1193 |
} else { |
1131 | 1194 |
v = _parent[v]; |
1132 | 1195 |
} |
1133 | 1196 |
} |
1134 | 1197 |
join = u; |
1135 | 1198 |
} |
1136 | 1199 |
|
1137 | 1200 |
// Find the leaving arc of the cycle and returns true if the |
1138 | 1201 |
// leaving arc is not the same as the entering arc |
1139 | 1202 |
bool findLeavingArc() { |
1140 | 1203 |
// Initialize first and second nodes according to the direction |
1141 | 1204 |
// of the cycle |
1142 | 1205 |
if (_state[in_arc] == STATE_LOWER) { |
1143 | 1206 |
first = _source[in_arc]; |
1144 | 1207 |
second = _target[in_arc]; |
1145 | 1208 |
} else { |
1146 | 1209 |
first = _target[in_arc]; |
1147 | 1210 |
second = _source[in_arc]; |
1148 | 1211 |
} |
1149 | 1212 |
delta = _cap[in_arc]; |
1150 | 1213 |
int result = 0; |
1151 | 1214 |
Value d; |
1152 | 1215 |
int e; |
1153 | 1216 |
|
1154 | 1217 |
// Search the cycle along the path form the first node to the root |
1155 | 1218 |
for (int u = first; u != join; u = _parent[u]) { |
1156 | 1219 |
e = _pred[u]; |
1157 | 1220 |
d = _forward[u] ? |
1158 | 1221 |
_flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]); |
1159 | 1222 |
if (d < delta) { |
1160 | 1223 |
delta = d; |
1161 | 1224 |
u_out = u; |
1162 | 1225 |
result = 1; |
1163 | 1226 |
} |
1164 | 1227 |
} |
1165 | 1228 |
// Search the cycle along the path form the second node to the root |
1166 | 1229 |
for (int u = second; u != join; u = _parent[u]) { |
1167 | 1230 |
e = _pred[u]; |
1168 | 1231 |
d = _forward[u] ? |
1169 | 1232 |
(_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e]; |
1170 | 1233 |
if (d <= delta) { |
1171 | 1234 |
delta = d; |
1172 | 1235 |
u_out = u; |
1173 | 1236 |
result = 2; |
1174 | 1237 |
} |
1175 | 1238 |
} |
1176 | 1239 |
|
1177 | 1240 |
if (result == 1) { |
1178 | 1241 |
u_in = first; |
1179 | 1242 |
v_in = second; |
1180 | 1243 |
} else { |
1181 | 1244 |
u_in = second; |
1182 | 1245 |
v_in = first; |
1183 | 1246 |
} |
1184 | 1247 |
return result != 0; |
1185 | 1248 |
} |
1186 | 1249 |
|
1187 | 1250 |
// Change _flow and _state vectors |
1188 | 1251 |
void changeFlow(bool change) { |
1189 | 1252 |
// Augment along the cycle |
1190 | 1253 |
if (delta > 0) { |
1191 | 1254 |
Value val = _state[in_arc] * delta; |
1192 | 1255 |
_flow[in_arc] += val; |
1193 | 1256 |
for (int u = _source[in_arc]; u != join; u = _parent[u]) { |
1194 | 1257 |
_flow[_pred[u]] += _forward[u] ? -val : val; |
1195 | 1258 |
} |
1196 | 1259 |
for (int u = _target[in_arc]; u != join; u = _parent[u]) { |
1197 | 1260 |
_flow[_pred[u]] += _forward[u] ? val : -val; |
1198 | 1261 |
} |
1199 | 1262 |
} |
1200 | 1263 |
// Update the state of the entering and leaving arcs |
1201 | 1264 |
if (change) { |
1202 | 1265 |
_state[in_arc] = STATE_TREE; |
1203 | 1266 |
_state[_pred[u_out]] = |
1204 | 1267 |
(_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; |
1205 | 1268 |
} else { |
1206 | 1269 |
_state[in_arc] = -_state[in_arc]; |
1207 | 1270 |
} |
1208 | 1271 |
} |
1209 | 1272 |
|
1210 | 1273 |
// Update the tree structure |
1211 | 1274 |
void updateTreeStructure() { |
1212 | 1275 |
int u, w; |
1213 | 1276 |
int old_rev_thread = _rev_thread[u_out]; |
1214 | 1277 |
int old_succ_num = _succ_num[u_out]; |
1215 | 1278 |
int old_last_succ = _last_succ[u_out]; |
1216 | 1279 |
v_out = _parent[u_out]; |
1217 | 1280 |
|
1218 | 1281 |
u = _last_succ[u_in]; // the last successor of u_in |
1219 | 1282 |
right = _thread[u]; // the node after it |
1220 | 1283 |
|
1221 | 1284 |
// Handle the case when old_rev_thread equals to v_in |
1222 | 1285 |
// (it also means that join and v_out coincide) |
1223 | 1286 |
if (old_rev_thread == v_in) { |
1224 | 1287 |
last = _thread[_last_succ[u_out]]; |
1225 | 1288 |
} else { |
1226 | 1289 |
last = _thread[v_in]; |
1227 | 1290 |
} |
1228 | 1291 |
|
1229 | 1292 |
// Update _thread and _parent along the stem nodes (i.e. the nodes |
1230 | 1293 |
// between u_in and u_out, whose parent have to be changed) |
1231 | 1294 |
_thread[v_in] = stem = u_in; |
1232 | 1295 |
_dirty_revs.clear(); |
1233 | 1296 |
_dirty_revs.push_back(v_in); |
1234 | 1297 |
par_stem = v_in; |
1235 | 1298 |
while (stem != u_out) { |
1236 | 1299 |
// Insert the next stem node into the thread list |
1237 | 1300 |
new_stem = _parent[stem]; |
1238 | 1301 |
_thread[u] = new_stem; |
1239 | 1302 |
_dirty_revs.push_back(u); |
1240 | 1303 |
|
1241 | 1304 |
// Remove the subtree of stem from the thread list |
1242 | 1305 |
w = _rev_thread[stem]; |
1243 | 1306 |
_thread[w] = right; |
1244 | 1307 |
_rev_thread[right] = w; |
1245 | 1308 |
|
1246 | 1309 |
// Change the parent node and shift stem nodes |
1247 | 1310 |
_parent[stem] = par_stem; |
1248 | 1311 |
par_stem = stem; |
1249 | 1312 |
stem = new_stem; |
1250 | 1313 |
|
1251 | 1314 |
// Update u and right |
1252 | 1315 |
u = _last_succ[stem] == _last_succ[par_stem] ? |
1253 | 1316 |
_rev_thread[par_stem] : _last_succ[stem]; |
1254 | 1317 |
right = _thread[u]; |
1255 | 1318 |
} |
1256 | 1319 |
_parent[u_out] = par_stem; |
1257 | 1320 |
_thread[u] = last; |
1258 | 1321 |
_rev_thread[last] = u; |
1259 | 1322 |
_last_succ[u_out] = u; |
1260 | 1323 |
|
1261 | 1324 |
// Remove the subtree of u_out from the thread list except for |
1262 | 1325 |
// the case when old_rev_thread equals to v_in |
1263 | 1326 |
// (it also means that join and v_out coincide) |
1264 | 1327 |
if (old_rev_thread != v_in) { |
1265 | 1328 |
_thread[old_rev_thread] = right; |
1266 | 1329 |
_rev_thread[right] = old_rev_thread; |
1267 | 1330 |
} |
1268 | 1331 |
|
1269 | 1332 |
// Update _rev_thread using the new _thread values |
1270 | 1333 |
for (int i = 0; i < int(_dirty_revs.size()); ++i) { |
1271 | 1334 |
u = _dirty_revs[i]; |
1272 | 1335 |
_rev_thread[_thread[u]] = u; |
1273 | 1336 |
} |
1274 | 1337 |
|
1275 | 1338 |
// Update _pred, _forward, _last_succ and _succ_num for the |
1276 | 1339 |
// stem nodes from u_out to u_in |
1277 | 1340 |
int tmp_sc = 0, tmp_ls = _last_succ[u_out]; |
1278 | 1341 |
u = u_out; |
1279 | 1342 |
while (u != u_in) { |
1280 | 1343 |
w = _parent[u]; |
1281 | 1344 |
_pred[u] = _pred[w]; |
1282 | 1345 |
_forward[u] = !_forward[w]; |
1283 | 1346 |
tmp_sc += _succ_num[u] - _succ_num[w]; |
1284 | 1347 |
_succ_num[u] = tmp_sc; |
1285 | 1348 |
_last_succ[w] = tmp_ls; |
1286 | 1349 |
u = w; |
1287 | 1350 |
} |
1288 | 1351 |
_pred[u_in] = in_arc; |
1289 | 1352 |
_forward[u_in] = (u_in == _source[in_arc]); |
1290 | 1353 |
_succ_num[u_in] = old_succ_num; |
1291 | 1354 |
|
1292 | 1355 |
// Set limits for updating _last_succ form v_in and v_out |
1293 | 1356 |
// towards the root |
1294 | 1357 |
int up_limit_in = -1; |
1295 | 1358 |
int up_limit_out = -1; |
1296 | 1359 |
if (_last_succ[join] == v_in) { |
1297 | 1360 |
up_limit_out = join; |
1298 | 1361 |
} else { |
1299 | 1362 |
up_limit_in = join; |
1300 | 1363 |
} |
1301 | 1364 |
|
1302 | 1365 |
// Update _last_succ from v_in towards the root |
1303 | 1366 |
for (u = v_in; u != up_limit_in && _last_succ[u] == v_in; |
1304 | 1367 |
u = _parent[u]) { |
1305 | 1368 |
_last_succ[u] = _last_succ[u_out]; |
1306 | 1369 |
} |
1307 | 1370 |
// Update _last_succ from v_out towards the root |
1308 | 1371 |
if (join != old_rev_thread && v_in != old_rev_thread) { |
1309 | 1372 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; |
1310 | 1373 |
u = _parent[u]) { |
1311 | 1374 |
_last_succ[u] = old_rev_thread; |
1312 | 1375 |
} |
1313 | 1376 |
} else { |
1314 | 1377 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; |
1315 | 1378 |
u = _parent[u]) { |
1316 | 1379 |
_last_succ[u] = _last_succ[u_out]; |
1317 | 1380 |
} |
1318 | 1381 |
} |
1319 | 1382 |
|
1320 | 1383 |
// Update _succ_num from v_in to join |
1321 | 1384 |
for (u = v_in; u != join; u = _parent[u]) { |
1322 | 1385 |
_succ_num[u] += old_succ_num; |
1323 | 1386 |
} |
1324 | 1387 |
// Update _succ_num from v_out to join |
1325 | 1388 |
for (u = v_out; u != join; u = _parent[u]) { |
1326 | 1389 |
_succ_num[u] -= old_succ_num; |
1327 | 1390 |
} |
1328 | 1391 |
} |
1329 | 1392 |
|
1330 | 1393 |
// Update potentials |
1331 | 1394 |
void updatePotential() { |
1332 | 1395 |
Cost sigma = _forward[u_in] ? |
1333 | 1396 |
_pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
1334 | 1397 |
_pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
1335 | 1398 |
// Update potentials in the subtree, which has been moved |
1336 | 1399 |
int end = _thread[_last_succ[u_in]]; |
1337 | 1400 |
for (int u = u_in; u != end; u = _thread[u]) { |
1338 | 1401 |
_pi[u] += sigma; |
1339 | 1402 |
} |
1340 | 1403 |
} |
1341 | 1404 |
|
1342 | 1405 |
// Execute the algorithm |
1343 | 1406 |
ProblemType start(PivotRule pivot_rule) { |
1344 | 1407 |
// Select the pivot rule implementation |
1345 | 1408 |
switch (pivot_rule) { |
1346 | 1409 |
case FIRST_ELIGIBLE: |
1347 | 1410 |
return start<FirstEligiblePivotRule>(); |
1348 | 1411 |
case BEST_ELIGIBLE: |
1349 | 1412 |
return start<BestEligiblePivotRule>(); |
1350 | 1413 |
case BLOCK_SEARCH: |
1351 | 1414 |
return start<BlockSearchPivotRule>(); |
1352 | 1415 |
case CANDIDATE_LIST: |
1353 | 1416 |
return start<CandidateListPivotRule>(); |
1354 | 1417 |
case ALTERING_LIST: |
1355 | 1418 |
return start<AlteringListPivotRule>(); |
1356 | 1419 |
} |
1357 | 1420 |
return INFEASIBLE; // avoid warning |
1358 | 1421 |
} |
1359 | 1422 |
|
1360 | 1423 |
template <typename PivotRuleImpl> |
1361 | 1424 |
ProblemType start() { |
1362 | 1425 |
PivotRuleImpl pivot(*this); |
1363 | 1426 |
|
1364 | 1427 |
// Execute the Network Simplex algorithm |
1365 | 1428 |
while (pivot.findEnteringArc()) { |
1366 | 1429 |
findJoinNode(); |
1367 | 1430 |
bool change = findLeavingArc(); |
1368 | 1431 |
if (delta >= INF) return UNBOUNDED; |
1369 | 1432 |
changeFlow(change); |
1370 | 1433 |
if (change) { |
1371 | 1434 |
updateTreeStructure(); |
1372 | 1435 |
updatePotential(); |
1373 | 1436 |
} |
1374 | 1437 |
} |
1375 | 1438 |
|
1376 | 1439 |
// Check feasibility |
1377 |
if (_sum_supply < 0) { |
|
1378 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1379 |
if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE; |
|
1380 |
} |
|
1381 |
} |
|
1382 |
else if (_sum_supply > 0) { |
|
1383 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1384 |
if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE; |
|
1385 |
} |
|
1386 |
} |
|
1387 |
else { |
|
1388 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1389 |
if (_flow[e] != 0) return INFEASIBLE; |
|
1390 |
} |
|
1440 |
for (int e = _search_arc_num; e != _all_arc_num; ++e) { |
|
1441 |
if (_flow[e] != 0) return INFEASIBLE; |
|
1391 | 1442 |
} |
1392 | 1443 |
|
1393 | 1444 |
// Transform the solution and the supply map to the original form |
1394 | 1445 |
if (_have_lower) { |
1395 | 1446 |
for (int i = 0; i != _arc_num; ++i) { |
1396 | 1447 |
Value c = _lower[i]; |
1397 | 1448 |
if (c != 0) { |
1398 | 1449 |
_flow[i] += c; |
1399 | 1450 |
_supply[_source[i]] += c; |
1400 | 1451 |
_supply[_target[i]] -= c; |
1401 | 1452 |
} |
1402 | 1453 |
} |
1403 | 1454 |
} |
1455 |
|
|
1456 |
// Shift potentials to meet the requirements of the GEQ/LEQ type |
|
1457 |
// optimality conditions |
|
1458 |
if (_sum_supply == 0) { |
|
1459 |
if (_stype == GEQ) { |
|
1460 |
Cost max_pot = std::numeric_limits<Cost>::min(); |
|
1461 |
for (int i = 0; i != _node_num; ++i) { |
|
1462 |
if (_pi[i] > max_pot) max_pot = _pi[i]; |
|
1463 |
} |
|
1464 |
if (max_pot > 0) { |
|
1465 |
for (int i = 0; i != _node_num; ++i) |
|
1466 |
_pi[i] -= max_pot; |
|
1467 |
} |
|
1468 |
} else { |
|
1469 |
Cost min_pot = std::numeric_limits<Cost>::max(); |
|
1470 |
for (int i = 0; i != _node_num; ++i) { |
|
1471 |
if (_pi[i] < min_pot) min_pot = _pi[i]; |
|
1472 |
} |
|
1473 |
if (min_pot < 0) { |
|
1474 |
for (int i = 0; i != _node_num; ++i) |
|
1475 |
_pi[i] -= min_pot; |
|
1476 |
} |
|
1477 |
} |
|
1478 |
} |
|
1404 | 1479 |
|
1405 | 1480 |
return OPTIMAL; |
1406 | 1481 |
} |
1407 | 1482 |
|
1408 | 1483 |
}; //class NetworkSimplex |
1409 | 1484 |
|
1410 | 1485 |
///@} |
1411 | 1486 |
|
1412 | 1487 |
} //namespace lemon |
1413 | 1488 |
|
1414 | 1489 |
#endif //LEMON_NETWORK_SIMPLEX_H |
0 comments (0 inline)