0
6
0
34
21
15
9
31
60
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 paths Path Structures |
230 | 230 |
@ingroup datas |
231 | 231 |
\brief %Path structures implemented in LEMON. |
232 | 232 |
|
233 | 233 |
This group contains the path structures implemented in LEMON. |
234 | 234 |
|
235 | 235 |
LEMON provides flexible data structures to work with paths. |
236 | 236 |
All of them have similar interfaces and they can be copied easily with |
237 | 237 |
assignment operators and copy constructors. This makes it easy and |
238 | 238 |
efficient to have e.g. the Dijkstra algorithm to store its result in |
239 | 239 |
any kind of path structure. |
240 | 240 |
|
241 | 241 |
\sa \ref concepts::Path "Path concept" |
242 | 242 |
*/ |
243 | 243 |
|
244 | 244 |
/** |
245 | 245 |
@defgroup heaps Heap Structures |
246 | 246 |
@ingroup datas |
247 | 247 |
\brief %Heap structures implemented in LEMON. |
248 | 248 |
|
249 | 249 |
This group contains the heap structures implemented in LEMON. |
250 | 250 |
|
251 | 251 |
LEMON provides several heap classes. They are efficient implementations |
252 | 252 |
of the abstract data type \e priority \e queue. They store items with |
253 | 253 |
specified values called \e priorities in such a way that finding and |
254 | 254 |
removing the item with minimum priority are efficient. |
255 | 255 |
The basic operations are adding and erasing items, changing the priority |
256 | 256 |
of an item, etc. |
257 | 257 |
|
258 | 258 |
Heaps are crucial in several algorithms, such as Dijkstra and Prim. |
259 | 259 |
The heap implementations have the same interface, thus any of them can be |
260 | 260 |
used easily in such algorithms. |
261 | 261 |
|
262 | 262 |
\sa \ref concepts::Heap "Heap concept" |
263 | 263 |
*/ |
264 | 264 |
|
265 | 265 |
/** |
266 | 266 |
@defgroup matrices Matrices |
267 | 267 |
@ingroup datas |
268 | 268 |
\brief Two dimensional data storages implemented in LEMON. |
269 | 269 |
|
270 | 270 |
This group contains two dimensional data storages implemented in LEMON. |
271 | 271 |
*/ |
272 | 272 |
|
273 | 273 |
/** |
274 | 274 |
@defgroup auxdat Auxiliary Data Structures |
275 | 275 |
@ingroup datas |
276 | 276 |
\brief Auxiliary data structures implemented in LEMON. |
277 | 277 |
|
278 | 278 |
This group contains some data structures implemented in LEMON in |
279 | 279 |
order to make it easier to implement combinatorial algorithms. |
280 | 280 |
*/ |
281 | 281 |
|
282 | 282 |
/** |
283 | 283 |
@defgroup geomdat Geometric Data Structures |
284 | 284 |
@ingroup auxdat |
285 | 285 |
\brief Geometric data structures implemented in LEMON. |
286 | 286 |
|
287 | 287 |
This group contains geometric data structures implemented in LEMON. |
288 | 288 |
|
289 | 289 |
- \ref lemon::dim2::Point "dim2::Point" implements a two dimensional |
290 | 290 |
vector with the usual operations. |
291 | 291 |
- \ref lemon::dim2::Box "dim2::Box" can be used to determine the |
292 | 292 |
rectangular bounding box of a set of \ref lemon::dim2::Point |
293 | 293 |
"dim2::Point"'s. |
294 | 294 |
*/ |
295 | 295 |
|
296 | 296 |
/** |
297 | 297 |
@defgroup matrices Matrices |
298 | 298 |
@ingroup auxdat |
299 | 299 |
\brief Two dimensional data storages implemented in LEMON. |
300 | 300 |
|
301 | 301 |
This group contains two dimensional data storages implemented in LEMON. |
302 | 302 |
*/ |
303 | 303 |
|
304 | 304 |
/** |
305 | 305 |
@defgroup algs Algorithms |
306 | 306 |
\brief This group contains the several algorithms |
307 | 307 |
implemented in LEMON. |
308 | 308 |
|
309 | 309 |
This group contains the several algorithms |
310 | 310 |
implemented in LEMON. |
311 | 311 |
*/ |
312 | 312 |
|
313 | 313 |
/** |
314 | 314 |
@defgroup search Graph Search |
315 | 315 |
@ingroup algs |
316 | 316 |
\brief Common graph search algorithms. |
317 | 317 |
|
318 | 318 |
This group contains the common graph search algorithms, namely |
319 |
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS) |
|
319 |
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS) |
|
320 |
\ref clrs01algorithms. |
|
320 | 321 |
*/ |
321 | 322 |
|
322 | 323 |
/** |
323 | 324 |
@defgroup shortest_path Shortest Path Algorithms |
324 | 325 |
@ingroup algs |
325 | 326 |
\brief Algorithms for finding shortest paths. |
326 | 327 |
|
327 |
This group contains the algorithms for finding shortest paths in digraphs |
|
328 |
This group contains the algorithms for finding shortest paths in digraphs |
|
329 |
\ref clrs01algorithms. |
|
328 | 330 |
|
329 | 331 |
- \ref Dijkstra algorithm for finding shortest paths from a source node |
330 | 332 |
when all arc lengths are non-negative. |
331 | 333 |
- \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths |
332 | 334 |
from a source node when arc lenghts can be either positive or negative, |
333 | 335 |
but the digraph should not contain directed cycles with negative total |
334 | 336 |
length. |
335 | 337 |
- \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms |
336 | 338 |
for solving the \e all-pairs \e shortest \e paths \e problem when arc |
337 | 339 |
lenghts can be either positive or negative, but the digraph should |
338 | 340 |
not contain directed cycles with negative total length. |
339 | 341 |
- \ref Suurballe A successive shortest path algorithm for finding |
340 | 342 |
arc-disjoint paths between two nodes having minimum total length. |
341 | 343 |
*/ |
342 | 344 |
|
343 | 345 |
/** |
344 | 346 |
@defgroup spantree Minimum Spanning Tree Algorithms |
345 | 347 |
@ingroup algs |
346 | 348 |
\brief Algorithms for finding minimum cost spanning trees and arborescences. |
347 | 349 |
|
348 | 350 |
This group contains the algorithms for finding minimum cost spanning |
349 |
trees and arborescences. |
|
351 |
trees and arborescences \ref clrs01algorithms. |
|
350 | 352 |
*/ |
351 | 353 |
|
352 | 354 |
/** |
353 | 355 |
@defgroup max_flow Maximum Flow Algorithms |
354 | 356 |
@ingroup algs |
355 | 357 |
\brief Algorithms for finding maximum flows. |
356 | 358 |
|
357 | 359 |
This group contains the algorithms for finding maximum flows and |
358 |
feasible circulations. |
|
360 |
feasible circulations \ref clrs01algorithms, \ref amo93networkflows. |
|
359 | 361 |
|
360 | 362 |
The \e maximum \e flow \e problem is to find a flow of maximum value between |
361 | 363 |
a single source and a single target. Formally, there is a \f$G=(V,A)\f$ |
362 | 364 |
digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and |
363 | 365 |
\f$s, t \in V\f$ source and target nodes. |
364 | 366 |
A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the |
365 | 367 |
following optimization problem. |
366 | 368 |
|
367 | 369 |
\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f] |
368 | 370 |
\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu) |
369 | 371 |
\quad \forall u\in V\setminus\{s,t\} \f] |
370 | 372 |
\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f] |
371 | 373 |
|
372 | 374 |
LEMON contains several algorithms for solving maximum flow problems: |
373 |
- \ref EdmondsKarp Edmonds-Karp algorithm. |
|
374 |
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm. |
|
375 |
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees. |
|
376 |
- \ref GoldbergTarjan Preflow push-relabel algorithm with dynamic trees. |
|
375 |
- \ref EdmondsKarp Edmonds-Karp algorithm |
|
376 |
\ref edmondskarp72theoretical. |
|
377 |
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm |
|
378 |
\ref goldberg88newapproach. |
|
379 |
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees |
|
380 |
\ref dinic70algorithm, \ref sleator83dynamic. |
|
381 |
- \ref GoldbergTarjan !Preflow push-relabel algorithm with dynamic trees |
|
382 |
\ref goldberg88newapproach, \ref sleator83dynamic. |
|
377 | 383 |
|
378 |
In most cases the \ref Preflow |
|
384 |
In most cases the \ref Preflow algorithm provides the |
|
379 | 385 |
fastest method for computing a maximum flow. All implementations |
380 | 386 |
also provide functions to query the minimum cut, which is the dual |
381 | 387 |
problem of maximum flow. |
382 | 388 |
|
383 | 389 |
\ref Circulation is a preflow push-relabel algorithm implemented directly |
384 | 390 |
for finding feasible circulations, which is a somewhat different problem, |
385 | 391 |
but it is strongly related to maximum flow. |
386 | 392 |
For more information, see \ref Circulation. |
387 | 393 |
*/ |
388 | 394 |
|
389 | 395 |
/** |
390 | 396 |
@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms |
391 | 397 |
@ingroup algs |
392 | 398 |
|
393 | 399 |
\brief Algorithms for finding minimum cost flows and circulations. |
394 | 400 |
|
395 | 401 |
This group contains the algorithms for finding minimum cost flows and |
396 |
circulations. For more information about this problem and its dual |
|
397 |
solution see \ref min_cost_flow "Minimum Cost Flow Problem". |
|
402 |
circulations \ref amo93networkflows. For more information about this |
|
403 |
problem and its dual solution, see \ref min_cost_flow |
|
404 |
"Minimum Cost Flow Problem". |
|
398 | 405 |
|
399 | 406 |
LEMON contains several algorithms for this problem. |
400 | 407 |
- \ref NetworkSimplex Primal Network Simplex algorithm with various |
401 |
pivot strategies. |
|
408 |
pivot strategies \ref dantzig63linearprog, \ref kellyoneill91netsimplex. |
|
402 | 409 |
- \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on |
403 |
cost scaling |
|
410 |
cost scaling \ref goldberg90approximation, \ref goldberg97efficient, |
|
411 |
\ref bunnagel98efficient. |
|
404 | 412 |
- \ref CapacityScaling Successive Shortest %Path algorithm with optional |
405 |
capacity scaling. |
|
406 |
- \ref CancelAndTighten The Cancel and Tighten algorithm. |
|
407 |
|
|
413 |
capacity scaling \ref edmondskarp72theoretical. |
|
414 |
- \ref CancelAndTighten The Cancel and Tighten algorithm |
|
415 |
\ref goldberg89cyclecanceling. |
|
416 |
- \ref CycleCanceling Cycle-Canceling algorithms |
|
417 |
\ref klein67primal, \ref goldberg89cyclecanceling. |
|
408 | 418 |
|
409 | 419 |
In general NetworkSimplex is the most efficient implementation, |
410 | 420 |
but in special cases other algorithms could be faster. |
411 | 421 |
For example, if the total supply and/or capacities are rather small, |
412 | 422 |
CapacityScaling is usually the fastest algorithm (without effective scaling). |
413 | 423 |
*/ |
414 | 424 |
|
415 | 425 |
/** |
416 | 426 |
@defgroup min_cut Minimum Cut Algorithms |
417 | 427 |
@ingroup algs |
418 | 428 |
|
419 | 429 |
\brief Algorithms for finding minimum cut in graphs. |
420 | 430 |
|
421 | 431 |
This group contains the algorithms for finding minimum cut in graphs. |
422 | 432 |
|
423 | 433 |
The \e minimum \e cut \e problem is to find a non-empty and non-complete |
424 | 434 |
\f$X\f$ subset of the nodes with minimum overall capacity on |
425 | 435 |
outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a |
426 | 436 |
\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum |
427 | 437 |
cut is the \f$X\f$ solution of the next optimization problem: |
428 | 438 |
|
429 | 439 |
\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}} |
430 | 440 |
\sum_{uv\in A: u\in X, v\not\in X}cap(uv) \f] |
431 | 441 |
|
432 | 442 |
LEMON contains several algorithms related to minimum cut problems: |
433 | 443 |
|
434 | 444 |
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut |
435 | 445 |
in directed graphs. |
436 | 446 |
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for |
437 | 447 |
calculating minimum cut in undirected graphs. |
438 | 448 |
- \ref GomoryHu "Gomory-Hu tree computation" for calculating |
439 | 449 |
all-pairs minimum cut in undirected graphs. |
440 | 450 |
|
441 | 451 |
If you want to find minimum cut just between two distinict nodes, |
442 | 452 |
see the \ref max_flow "maximum flow problem". |
443 | 453 |
*/ |
444 | 454 |
|
445 | 455 |
/** |
446 | 456 |
@defgroup matching Matching Algorithms |
447 | 457 |
@ingroup algs |
448 | 458 |
\brief Algorithms for finding matchings in graphs and bipartite graphs. |
449 | 459 |
|
450 | 460 |
This group contains the algorithms for calculating |
451 | 461 |
matchings in graphs and bipartite graphs. The general matching problem is |
452 | 462 |
finding a subset of the edges for which each node has at most one incident |
453 | 463 |
edge. |
454 | 464 |
|
455 | 465 |
There are several different algorithms for calculate matchings in |
456 | 466 |
graphs. The matching problems in bipartite graphs are generally |
457 | 467 |
easier than in general graphs. The goal of the matching optimization |
458 | 468 |
can be finding maximum cardinality, maximum weight or minimum cost |
459 | 469 |
matching. The search can be constrained to find perfect or |
460 | 470 |
maximum cardinality matching. |
461 | 471 |
|
462 | 472 |
The matching algorithms implemented in LEMON: |
463 | 473 |
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm |
464 | 474 |
for calculating maximum cardinality matching in bipartite graphs. |
465 | 475 |
- \ref PrBipartiteMatching Push-relabel algorithm |
466 | 476 |
for calculating maximum cardinality matching in bipartite graphs. |
467 | 477 |
- \ref MaxWeightedBipartiteMatching |
468 | 478 |
Successive shortest path algorithm for calculating maximum weighted |
469 | 479 |
matching and maximum weighted bipartite matching in bipartite graphs. |
470 | 480 |
- \ref MinCostMaxBipartiteMatching |
471 | 481 |
Successive shortest path algorithm for calculating minimum cost maximum |
472 | 482 |
matching in bipartite graphs. |
473 | 483 |
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating |
474 | 484 |
maximum cardinality matching in general graphs. |
475 | 485 |
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating |
476 | 486 |
maximum weighted matching in general graphs. |
477 | 487 |
- \ref MaxWeightedPerfectMatching |
478 | 488 |
Edmond's blossom shrinking algorithm for calculating maximum weighted |
479 | 489 |
perfect matching in general graphs. |
480 | 490 |
|
481 | 491 |
\image html bipartite_matching.png |
482 | 492 |
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth |
483 | 493 |
*/ |
484 | 494 |
|
485 | 495 |
/** |
486 | 496 |
@defgroup graph_properties Connectivity and Other Graph Properties |
487 | 497 |
@ingroup algs |
488 | 498 |
\brief Algorithms for discovering the graph properties |
489 | 499 |
|
490 | 500 |
This group contains the algorithms for discovering the graph properties |
491 | 501 |
like connectivity, bipartiteness, euler property, simplicity etc. |
492 | 502 |
|
493 | 503 |
\image html connected_components.png |
494 | 504 |
\image latex connected_components.eps "Connected components" width=\textwidth |
495 | 505 |
*/ |
496 | 506 |
|
497 | 507 |
/** |
498 | 508 |
@defgroup planar Planarity Embedding and Drawing |
499 | 509 |
@ingroup algs |
500 | 510 |
\brief Algorithms for planarity checking, embedding and drawing |
501 | 511 |
|
502 | 512 |
This group contains the algorithms for planarity checking, |
503 | 513 |
embedding and drawing. |
504 | 514 |
|
505 | 515 |
\image html planar.png |
506 | 516 |
\image latex planar.eps "Plane graph" width=\textwidth |
507 | 517 |
*/ |
508 | 518 |
|
509 | 519 |
/** |
510 | 520 |
@defgroup approx Approximation Algorithms |
511 | 521 |
@ingroup algs |
512 | 522 |
\brief Approximation algorithms. |
513 | 523 |
|
514 | 524 |
This group contains the approximation and heuristic algorithms |
515 | 525 |
implemented in LEMON. |
516 | 526 |
*/ |
517 | 527 |
|
518 | 528 |
/** |
519 | 529 |
@defgroup auxalg Auxiliary Algorithms |
520 | 530 |
@ingroup algs |
521 | 531 |
\brief Auxiliary algorithms implemented in LEMON. |
522 | 532 |
|
523 | 533 |
This group contains some algorithms implemented in LEMON |
524 | 534 |
in order to make it easier to implement complex algorithms. |
525 | 535 |
*/ |
526 | 536 |
|
527 | 537 |
/** |
528 | 538 |
@defgroup gen_opt_group General Optimization Tools |
529 | 539 |
\brief This group contains some general optimization frameworks |
530 | 540 |
implemented in LEMON. |
531 | 541 |
|
532 | 542 |
This group contains some general optimization frameworks |
533 | 543 |
implemented in LEMON. |
534 | 544 |
*/ |
535 | 545 |
|
536 | 546 |
/** |
537 |
@defgroup lp_group |
|
547 |
@defgroup lp_group LP and MIP Solvers |
|
538 | 548 |
@ingroup gen_opt_group |
539 |
\brief |
|
549 |
\brief LP and MIP solver interfaces for LEMON. |
|
540 | 550 |
|
541 |
This group contains Lp and Mip solver interfaces for LEMON. The |
|
542 |
various LP solvers could be used in the same manner with this |
|
543 |
|
|
551 |
This group contains LP and MIP solver interfaces for LEMON. |
|
552 |
Various LP solvers could be used in the same manner with this |
|
553 |
high-level interface. |
|
554 |
|
|
555 |
The currently supported solvers are \ref glpk, \ref clp, \ref cbc, |
|
556 |
\ref cplex, \ref soplex. |
|
544 | 557 |
*/ |
545 | 558 |
|
546 | 559 |
/** |
547 | 560 |
@defgroup lp_utils Tools for Lp and Mip Solvers |
548 | 561 |
@ingroup lp_group |
549 | 562 |
\brief Helper tools to the Lp and Mip solvers. |
550 | 563 |
|
551 | 564 |
This group adds some helper tools to general optimization framework |
552 | 565 |
implemented in LEMON. |
553 | 566 |
*/ |
554 | 567 |
|
555 | 568 |
/** |
556 | 569 |
@defgroup metah Metaheuristics |
557 | 570 |
@ingroup gen_opt_group |
558 | 571 |
\brief Metaheuristics for LEMON library. |
559 | 572 |
|
560 | 573 |
This group contains some metaheuristic optimization tools. |
561 | 574 |
*/ |
562 | 575 |
|
563 | 576 |
/** |
564 | 577 |
@defgroup utils Tools and Utilities |
565 | 578 |
\brief Tools and utilities for programming in LEMON |
566 | 579 |
|
567 | 580 |
Tools and utilities for programming in LEMON. |
568 | 581 |
*/ |
569 | 582 |
|
570 | 583 |
/** |
571 | 584 |
@defgroup gutils Basic Graph Utilities |
572 | 585 |
@ingroup utils |
573 | 586 |
\brief Simple basic graph utilities. |
574 | 587 |
|
575 | 588 |
This group contains some simple basic graph utilities. |
576 | 589 |
*/ |
577 | 590 |
|
578 | 591 |
/** |
579 | 592 |
@defgroup misc Miscellaneous Tools |
580 | 593 |
@ingroup utils |
581 | 594 |
\brief Tools for development, debugging and testing. |
582 | 595 |
|
583 | 596 |
This group contains several useful tools for development, |
584 | 597 |
debugging and testing. |
585 | 598 |
*/ |
586 | 599 |
|
587 | 600 |
/** |
588 | 601 |
@defgroup timecount Time Measuring and Counting |
589 | 602 |
@ingroup misc |
590 | 603 |
\brief Simple tools for measuring the performance of algorithms. |
591 | 604 |
|
592 | 605 |
This group contains simple tools for measuring the performance |
593 | 606 |
of algorithms. |
594 | 607 |
*/ |
595 | 608 |
|
596 | 609 |
/** |
597 | 610 |
@defgroup exceptions Exceptions |
598 | 611 |
@ingroup utils |
599 | 612 |
\brief Exceptions defined in LEMON. |
600 | 613 |
|
601 | 614 |
This group contains the exceptions defined in LEMON. |
602 | 615 |
*/ |
603 | 616 |
|
604 | 617 |
/** |
605 | 618 |
@defgroup io_group Input-Output |
606 | 619 |
\brief Graph Input-Output methods |
607 | 620 |
|
608 | 621 |
This group contains the tools for importing and exporting graphs |
609 | 622 |
and graph related data. Now it supports the \ref lgf-format |
610 | 623 |
"LEMON Graph Format", the \c DIMACS format and the encapsulated |
611 | 624 |
postscript (EPS) format. |
612 | 625 |
*/ |
613 | 626 |
|
614 | 627 |
/** |
615 | 628 |
@defgroup lemon_io LEMON Graph Format |
616 | 629 |
@ingroup io_group |
617 | 630 |
\brief Reading and writing LEMON Graph Format. |
618 | 631 |
|
619 | 632 |
This group contains methods for reading and writing |
620 | 633 |
\ref lgf-format "LEMON Graph Format". |
621 | 634 |
*/ |
622 | 635 |
|
623 | 636 |
/** |
624 | 637 |
@defgroup eps_io Postscript Exporting |
625 | 638 |
@ingroup io_group |
626 | 639 |
\brief General \c EPS drawer and graph exporter |
627 | 640 |
|
628 | 641 |
This group contains general \c EPS drawing methods and special |
629 | 642 |
graph exporting tools. |
630 | 643 |
*/ |
631 | 644 |
|
632 | 645 |
/** |
633 | 646 |
@defgroup dimacs_group DIMACS Format |
634 | 647 |
@ingroup io_group |
635 | 648 |
\brief Read and write files in DIMACS format |
636 | 649 |
|
637 | 650 |
Tools to read a digraph from or write it to a file in DIMACS format data. |
638 | 651 |
*/ |
639 | 652 |
|
640 | 653 |
/** |
641 | 654 |
@defgroup nauty_group NAUTY Format |
642 | 655 |
@ingroup io_group |
643 | 656 |
\brief Read \e Nauty format |
644 | 657 |
|
645 | 658 |
Tool to read graphs from \e Nauty format data. |
646 | 659 |
*/ |
647 | 660 |
|
648 | 661 |
/** |
649 | 662 |
@defgroup concept Concepts |
650 | 663 |
\brief Skeleton classes and concept checking classes |
651 | 664 |
|
652 | 665 |
This group contains the data/algorithm skeletons and concept checking |
653 | 666 |
classes implemented in LEMON. |
654 | 667 |
|
655 | 668 |
The purpose of the classes in this group is fourfold. |
656 | 669 |
|
657 | 670 |
- These classes contain the documentations of the %concepts. In order |
658 | 671 |
to avoid document multiplications, an implementation of a concept |
659 | 672 |
simply refers to the corresponding concept class. |
660 | 673 |
|
661 | 674 |
- These classes declare every functions, <tt>typedef</tt>s etc. an |
662 | 675 |
implementation of the %concepts should provide, however completely |
663 | 676 |
without implementations and real data structures behind the |
664 | 677 |
interface. On the other hand they should provide nothing else. All |
665 | 678 |
the algorithms working on a data structure meeting a certain concept |
666 | 679 |
should compile with these classes. (Though it will not run properly, |
667 | 680 |
of course.) In this way it is easily to check if an algorithm |
668 | 681 |
doesn't use any extra feature of a certain implementation. |
669 | 682 |
|
670 | 683 |
- The concept descriptor classes also provide a <em>checker class</em> |
671 | 684 |
that makes it possible to check whether a certain implementation of a |
672 | 685 |
concept indeed provides all the required features. |
673 | 686 |
|
674 | 687 |
- Finally, They can serve as a skeleton of a new implementation of a concept. |
675 | 688 |
*/ |
676 | 689 |
|
677 | 690 |
/** |
678 | 691 |
@defgroup graph_concepts Graph Structure Concepts |
679 | 692 |
@ingroup concept |
680 | 693 |
\brief Skeleton and concept checking classes for graph structures |
681 | 694 |
|
682 | 695 |
This group contains the skeletons and concept checking classes of |
683 | 696 |
graph structures. |
684 | 697 |
*/ |
685 | 698 |
|
686 | 699 |
/** |
687 | 700 |
@defgroup map_concepts Map Concepts |
688 | 701 |
@ingroup concept |
689 | 702 |
\brief Skeleton and concept checking classes for maps |
690 | 703 |
|
691 | 704 |
This group contains the skeletons and concept checking classes of maps. |
692 | 705 |
*/ |
693 | 706 |
|
694 | 707 |
/** |
695 | 708 |
@defgroup tools Standalone Utility Applications |
696 | 709 |
|
697 | 710 |
Some utility applications are listed here. |
698 | 711 |
|
699 | 712 |
The standard compilation procedure (<tt>./configure;make</tt>) will compile |
700 | 713 |
them, as well. |
701 | 714 |
*/ |
702 | 715 |
|
703 | 716 |
/** |
704 | 717 |
\anchor demoprograms |
705 | 718 |
|
706 | 719 |
@defgroup demos Demo Programs |
707 | 720 |
|
708 | 721 |
Some demo programs are listed here. Their full source codes can be found in |
709 | 722 |
the \c demo subdirectory of the source tree. |
710 | 723 |
|
711 | 724 |
In order to compile them, use the <tt>make demo</tt> or the |
712 | 725 |
<tt>make check</tt> commands. |
713 | 726 |
*/ |
714 | 727 |
|
715 | 728 |
} |
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 |
/** |
20 | 20 |
\mainpage LEMON Documentation |
21 | 21 |
|
22 | 22 |
\section intro Introduction |
23 | 23 |
|
24 |
\subsection whatis What is LEMON |
|
25 |
|
|
26 |
LEMON stands for <b>L</b>ibrary for <b>E</b>fficient <b>M</b>odeling |
|
27 |
and <b>O</b>ptimization in <b>N</b>etworks. |
|
28 |
It is a C++ template |
|
29 |
library aimed at combinatorial optimization tasks which |
|
30 |
often involve in working |
|
31 |
with graphs. |
|
24 |
<b>LEMON</b> stands for <i><b>L</b>ibrary for <b>E</b>fficient <b>M</b>odeling |
|
25 |
and <b>O</b>ptimization in <b>N</b>etworks</i>. |
|
26 |
It is a C++ template library providing efficient implementation of common |
|
27 |
data structures and algorithms with focus on combinatorial optimization |
|
28 |
problems in graphs and networks. |
|
32 | 29 |
|
33 | 30 |
<b> |
34 | 31 |
LEMON is an <a class="el" href="http://opensource.org/">open source</a> |
35 | 32 |
project. |
36 | 33 |
You are free to use it in your commercial or |
37 | 34 |
non-commercial applications under very permissive |
38 | 35 |
\ref license "license terms". |
39 | 36 |
</b> |
40 | 37 |
|
41 |
|
|
38 |
The project is maintained by the |
|
39 |
<a href="http://www.cs.elte.hu/egres/">Egerváry Research Group on |
|
40 |
Combinatorial Optimization</a> \ref egres |
|
41 |
at the Operations Research Department of the |
|
42 |
<a href="http://www.elte.hu/">Eötvös Loránd University, |
|
43 |
Budapest</a>, Hungary. |
|
44 |
LEMON is also a member of the <a href="http://www.coin-or.org/">COIN-OR</a> |
|
45 |
initiative \ref coinor. |
|
46 |
|
|
47 |
\section howtoread How to Read the Documentation |
|
42 | 48 |
|
43 | 49 |
If you would like to get to know the library, see |
44 | 50 |
<a class="el" href="http://lemon.cs.elte.hu/pub/tutorial/">LEMON Tutorial</a>. |
45 | 51 |
|
46 | 52 |
If you know what you are looking for, then try to find it under the |
47 | 53 |
<a class="el" href="modules.html">Modules</a> section. |
48 | 54 |
|
49 | 55 |
If you are a user of the old (0.x) series of LEMON, please check out the |
50 | 56 |
\ref migration "Migration Guide" for the backward incompatibilities. |
51 | 57 |
*/ |
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 |
\page min_cost_flow Minimum Cost Flow Problem |
23 | 23 |
|
24 | 24 |
\section mcf_def Definition (GEQ form) |
25 | 25 |
|
26 | 26 |
The \e minimum \e cost \e flow \e problem is to find a feasible flow of |
27 | 27 |
minimum total cost from a set of supply nodes to a set of demand nodes |
28 | 28 |
in a network with capacity constraints (lower and upper bounds) |
29 |
and arc costs. |
|
29 |
and arc costs \ref amo93networkflows. |
|
30 | 30 |
|
31 | 31 |
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$, |
32 | 32 |
\f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and |
33 | 33 |
upper bounds for the flow values on the arcs, for which |
34 | 34 |
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$, |
35 | 35 |
\f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow |
36 | 36 |
on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the |
37 | 37 |
signed supply values of the nodes. |
38 | 38 |
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ |
39 | 39 |
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with |
40 | 40 |
\f$-sup(u)\f$ demand. |
41 | 41 |
A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution |
42 | 42 |
of the following optimization problem. |
43 | 43 |
|
44 | 44 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
45 | 45 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq |
46 | 46 |
sup(u) \quad \forall u\in V \f] |
47 | 47 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
48 | 48 |
|
49 | 49 |
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be |
50 | 50 |
zero or negative in order to have a feasible solution (since the sum |
51 | 51 |
of the expressions on the left-hand side of the inequalities is zero). |
52 | 52 |
It means that the total demand must be greater or equal to the total |
53 | 53 |
supply and all the supplies have to be carried out from the supply nodes, |
54 | 54 |
but there could be demands that are not satisfied. |
55 | 55 |
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand |
56 | 56 |
constraints have to be satisfied with equality, i.e. all demands |
57 | 57 |
have to be satisfied and all supplies have to be used. |
58 | 58 |
|
59 | 59 |
|
60 | 60 |
\section mcf_algs Algorithms |
61 | 61 |
|
62 | 62 |
LEMON contains several algorithms for solving this problem, for more |
63 | 63 |
information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms". |
64 | 64 |
|
65 | 65 |
A feasible solution for this problem can be found using \ref Circulation. |
66 | 66 |
|
67 | 67 |
|
68 | 68 |
\section mcf_dual Dual Solution |
69 | 69 |
|
70 | 70 |
The dual solution of the minimum cost flow problem is represented by |
71 | 71 |
node potentials \f$\pi: V\rightarrow\mathbf{R}\f$. |
72 | 72 |
An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal |
73 | 73 |
if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials |
74 | 74 |
the following \e complementary \e slackness optimality conditions hold. |
75 | 75 |
|
76 | 76 |
- For all \f$uv\in A\f$ arcs: |
77 | 77 |
- if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; |
78 | 78 |
- if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; |
79 | 79 |
- if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. |
80 | 80 |
- For all \f$u\in V\f$ nodes: |
81 | 81 |
- \f$\pi(u)<=0\f$; |
82 | 82 |
- if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$, |
83 | 83 |
then \f$\pi(u)=0\f$. |
84 | 84 |
|
85 | 85 |
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc |
86 | 86 |
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e. |
87 | 87 |
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f] |
88 | 88 |
|
89 | 89 |
All algorithms provide dual solution (node potentials), as well, |
90 | 90 |
if an optimal flow is found. |
91 | 91 |
|
92 | 92 |
|
93 | 93 |
\section mcf_eq Equality Form |
94 | 94 |
|
95 | 95 |
The above \ref mcf_def "definition" is actually more general than the |
96 | 96 |
usual formulation of the minimum cost flow problem, in which strict |
97 | 97 |
equalities are required in the supply/demand contraints. |
98 | 98 |
|
99 | 99 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
100 | 100 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) = |
101 | 101 |
sup(u) \quad \forall u\in V \f] |
102 | 102 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
103 | 103 |
|
104 | 104 |
However if the sum of the supply values is zero, then these two problems |
105 | 105 |
are equivalent. |
106 | 106 |
The \ref min_cost_flow_algs "algorithms" in LEMON support the general |
107 | 107 |
form, so if you need the equality form, you have to ensure this additional |
108 | 108 |
contraint manually. |
109 | 109 |
|
110 | 110 |
|
111 | 111 |
\section mcf_leq Opposite Inequalites (LEQ Form) |
112 | 112 |
|
113 | 113 |
Another possible definition of the minimum cost flow problem is |
114 | 114 |
when there are <em>"less or equal"</em> (LEQ) supply/demand constraints, |
115 | 115 |
instead of the <em>"greater or equal"</em> (GEQ) constraints. |
116 | 116 |
|
117 | 117 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
118 | 118 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq |
119 | 119 |
sup(u) \quad \forall u\in V \f] |
120 | 120 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
121 | 121 |
|
122 | 122 |
It means that the total demand must be less or equal to the |
123 | 123 |
total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or |
124 | 124 |
positive) and all the demands have to be satisfied, but there |
125 | 125 |
could be supplies that are not carried out from the supply |
126 | 126 |
nodes. |
127 | 127 |
The equality form is also a special case of this form, of course. |
128 | 128 |
|
129 | 129 |
You could easily transform this case to the \ref mcf_def "GEQ form" |
130 | 130 |
of the problem by reversing the direction of the arcs and taking the |
131 | 131 |
negative of the supply values (e.g. using \ref ReverseDigraph and |
132 | 132 |
\ref NegMap adaptors). |
133 | 133 |
However \ref NetworkSimplex algorithm also supports this form directly |
134 | 134 |
for the sake of convenience. |
135 | 135 |
|
136 | 136 |
Note that the optimality conditions for this supply constraint type are |
137 | 137 |
slightly differ from the conditions that are discussed for the GEQ form, |
138 | 138 |
namely the potentials have to be non-negative instead of non-positive. |
139 | 139 |
An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem |
140 | 140 |
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ |
141 | 141 |
node potentials the following conditions hold. |
142 | 142 |
|
143 | 143 |
- For all \f$uv\in A\f$ arcs: |
144 | 144 |
- if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$; |
145 | 145 |
- if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$; |
146 | 146 |
- if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$. |
147 | 147 |
- For all \f$u\in V\f$ nodes: |
148 | 148 |
- \f$\pi(u)>=0\f$; |
149 | 149 |
- if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$, |
150 | 150 |
then \f$\pi(u)=0\f$. |
151 | 151 |
|
152 | 152 |
*/ |
153 | 153 |
} |
1 | 1 |
%%%%% Defining LEMON %%%%% |
2 | 2 |
|
3 | 3 |
@misc{lemon, |
4 | 4 |
key = {LEMON}, |
5 | 5 |
title = {{LEMON} -- {L}ibrary for {E}fficient {M}odeling and |
6 | 6 |
{O}ptimization in {N}etworks}, |
7 | 7 |
howpublished = {\url{http://lemon.cs.elte.hu/}}, |
8 | 8 |
year = 2009 |
9 | 9 |
} |
10 | 10 |
|
11 | 11 |
@misc{egres, |
12 | 12 |
key = {EGRES}, |
13 | 13 |
title = {{EGRES} -- {E}gerv{\'a}ry {R}esearch {G}roup on |
14 | 14 |
{C}ombinatorial {O}ptimization}, |
15 | 15 |
url = {http://www.cs.elte.hu/egres/} |
16 | 16 |
} |
17 | 17 |
|
18 | 18 |
@misc{coinor, |
19 | 19 |
key = {COIN-OR}, |
20 | 20 |
title = {{COIN-OR} -- {C}omputational {I}nfrastructure for |
21 | 21 |
{O}perations {R}esearch}, |
22 | 22 |
url = {http://www.coin-or.org/} |
23 | 23 |
} |
24 | 24 |
|
25 | 25 |
|
26 | 26 |
%%%%% Other libraries %%%%%% |
27 | 27 |
|
28 | 28 |
@misc{boost, |
29 | 29 |
key = {Boost}, |
30 | 30 |
title = {{B}oost {C++} {L}ibraries}, |
31 | 31 |
url = {http://www.boost.org/} |
32 | 32 |
} |
33 | 33 |
|
34 | 34 |
@book{bglbook, |
35 | 35 |
author = {Jeremy G. Siek and Lee-Quan Lee and Andrew |
36 | 36 |
Lumsdaine}, |
37 | 37 |
title = {The Boost Graph Library: User Guide and Reference |
38 | 38 |
Manual}, |
39 | 39 |
publisher = {Addison-Wesley}, |
40 | 40 |
year = 2002 |
41 | 41 |
} |
42 | 42 |
|
43 | 43 |
@misc{leda, |
44 | 44 |
key = {LEDA}, |
45 | 45 |
title = {{LEDA} -- {L}ibrary of {E}fficient {D}ata {T}ypes and |
46 | 46 |
{A}lgorithms}, |
47 | 47 |
url = {http://www.algorithmic-solutions.com/} |
48 | 48 |
} |
49 | 49 |
|
50 | 50 |
@book{ledabook, |
51 | 51 |
author = {Kurt Mehlhorn and Stefan N{\"a}her}, |
52 | 52 |
title = {{LEDA}: {A} platform for combinatorial and geometric |
53 | 53 |
computing}, |
54 | 54 |
isbn = {0-521-56329-1}, |
55 | 55 |
publisher = {Cambridge University Press}, |
56 | 56 |
address = {New York, NY, USA}, |
57 | 57 |
year = 1999 |
58 | 58 |
} |
59 | 59 |
|
60 | 60 |
|
61 | 61 |
%%%%% Tools that LEMON depends on %%%%% |
62 | 62 |
|
63 | 63 |
@misc{cmake, |
64 | 64 |
key = {CMake}, |
65 | 65 |
title = {{CMake} -- {C}ross {P}latform {M}ake}, |
66 | 66 |
url = {http://www.cmake.org/} |
67 | 67 |
} |
68 | 68 |
|
69 | 69 |
@misc{doxygen, |
70 | 70 |
key = {Doxygen}, |
71 | 71 |
title = {{Doxygen} -- {S}ource code documentation generator |
72 | 72 |
tool}, |
73 | 73 |
url = {http://www.doxygen.org/} |
74 | 74 |
} |
75 | 75 |
|
76 | 76 |
|
77 | 77 |
%%%%% LP/MIP libraries %%%%% |
78 | 78 |
|
79 | 79 |
@misc{glpk, |
80 | 80 |
key = {GLPK}, |
81 | 81 |
title = {{GLPK} -- {GNU} {L}inear {P}rogramming {K}it}, |
82 | 82 |
url = {http://www.gnu.org/software/glpk/} |
83 | 83 |
} |
84 | 84 |
|
85 | 85 |
@misc{clp, |
86 | 86 |
key = {Clp}, |
87 | 87 |
title = {{Clp} -- {Coin-Or} {L}inear {P}rogramming}, |
88 | 88 |
url = {http://projects.coin-or.org/Clp/} |
89 | 89 |
} |
90 | 90 |
|
91 | 91 |
@misc{cbc, |
92 | 92 |
key = {Cbc}, |
93 | 93 |
title = {{Cbc} -- {Coin-Or} {B}ranch and {C}ut}, |
94 | 94 |
url = {http://projects.coin-or.org/Cbc/} |
95 | 95 |
} |
96 | 96 |
|
97 | 97 |
@misc{cplex, |
98 | 98 |
key = {CPLEX}, |
99 | 99 |
title = {{ILOG} {CPLEX}}, |
100 | 100 |
url = {http://www.ilog.com/} |
101 | 101 |
} |
102 | 102 |
|
103 | 103 |
@misc{soplex, |
104 | 104 |
key = {SoPlex}, |
105 | 105 |
title = {{SoPlex} -- {T}he {S}equential {O}bject-{O}riented |
106 | 106 |
{S}implex}, |
107 | 107 |
url = {http://soplex.zib.de/} |
108 | 108 |
} |
109 | 109 |
|
110 | 110 |
|
111 | 111 |
%%%%% General books %%%%% |
112 | 112 |
|
113 | 113 |
@book{amo93networkflows, |
114 | 114 |
author = {Ravindra K. Ahuja and Thomas L. Magnanti and James |
115 | 115 |
B. Orlin}, |
116 | 116 |
title = {Network Flows: Theory, Algorithms, and Applications}, |
117 | 117 |
publisher = {Prentice-Hall, Inc.}, |
118 | 118 |
year = 1993, |
119 | 119 |
month = feb, |
120 | 120 |
isbn = {978-0136175490} |
121 | 121 |
} |
122 | 122 |
|
123 | 123 |
@book{schrijver03combinatorial, |
124 | 124 |
author = {Alexander Schrijver}, |
125 | 125 |
title = {Combinatorial Optimization: Polyhedra and Efficiency}, |
126 | 126 |
publisher = {Springer-Verlag}, |
127 | 127 |
year = 2003, |
128 | 128 |
isbn = {978-3540443896} |
129 | 129 |
} |
130 | 130 |
|
131 | 131 |
@book{clrs01algorithms, |
132 | 132 |
author = {Thomas H. Cormen and Charles E. Leiserson and Ronald |
133 | 133 |
L. Rivest and Clifford Stein}, |
134 | 134 |
title = {Introduction to Algorithms}, |
135 | 135 |
publisher = {The MIT Press}, |
136 | 136 |
year = 2001, |
137 | 137 |
edition = {2nd} |
138 | 138 |
} |
139 | 139 |
|
140 | 140 |
@book{stroustrup00cpp, |
141 | 141 |
author = {Bjarne Stroustrup}, |
142 | 142 |
title = {The C++ Programming Language}, |
143 | 143 |
edition = {3rd}, |
144 | 144 |
publisher = {Addison-Wesley Professional}, |
145 | 145 |
isbn = 0201700735, |
146 | 146 |
month = {February}, |
147 | 147 |
year = 2000 |
148 | 148 |
} |
149 | 149 |
|
150 | 150 |
|
151 | 151 |
%%%%% Maximum flow algorithms %%%%% |
152 | 152 |
|
153 |
@ |
|
153 |
@article{edmondskarp72theoretical, |
|
154 |
author = {Jack Edmonds and Richard M. Karp}, |
|
155 |
title = {Theoretical improvements in algorithmic efficiency |
|
156 |
for network flow problems}, |
|
157 |
journal = {Journal of the ACM}, |
|
158 |
year = 1972, |
|
159 |
volume = 19, |
|
160 |
number = 2, |
|
161 |
pages = {248-264} |
|
162 |
} |
|
163 |
|
|
164 |
@article{goldberg88newapproach, |
|
154 | 165 |
author = {Andrew V. Goldberg and Robert E. Tarjan}, |
155 | 166 |
title = {A new approach to the maximum flow problem}, |
156 |
booktitle = {STOC '86: Proceedings of the Eighteenth Annual ACM |
|
157 |
Symposium on Theory of Computing}, |
|
158 |
year = 1986, |
|
159 |
publisher = {ACM Press}, |
|
160 |
address = {New York, NY}, |
|
161 |
pages = {136-146} |
|
167 |
journal = {Journal of the ACM}, |
|
168 |
year = 1988, |
|
169 |
volume = 35, |
|
170 |
number = 4, |
|
171 |
pages = {921-940} |
|
162 | 172 |
} |
163 | 173 |
|
164 | 174 |
@article{dinic70algorithm, |
165 | 175 |
author = {E. A. Dinic}, |
166 | 176 |
title = {Algorithm for solution of a problem of maximum flow |
167 | 177 |
in a network with power estimation}, |
168 | 178 |
journal = {Soviet Math. Doklady}, |
169 | 179 |
year = 1970, |
170 | 180 |
volume = 11, |
171 | 181 |
pages = {1277-1280} |
172 | 182 |
} |
173 | 183 |
|
174 | 184 |
@article{goldberg08partial, |
175 | 185 |
author = {Andrew V. Goldberg}, |
176 | 186 |
title = {The Partial Augment-Relabel Algorithm for the |
177 | 187 |
Maximum Flow Problem}, |
178 | 188 |
journal = {16th Annual European Symposium on Algorithms}, |
179 | 189 |
year = 2008, |
180 | 190 |
pages = {466-477} |
181 | 191 |
} |
182 | 192 |
|
183 | 193 |
@article{sleator83dynamic, |
184 | 194 |
author = {Daniel D. Sleator and Robert E. Tarjan}, |
185 | 195 |
title = {A data structure for dynamic trees}, |
186 | 196 |
journal = {Journal of Computer and System Sciences}, |
187 | 197 |
year = 1983, |
188 | 198 |
volume = 26, |
189 | 199 |
number = 3, |
190 | 200 |
pages = {362-391} |
191 | 201 |
} |
192 | 202 |
|
193 | 203 |
|
194 | 204 |
%%%%% Minimum mean cycle algorithms %%%%% |
195 | 205 |
|
196 | 206 |
@article{karp78characterization, |
197 | 207 |
author = {Richard M. Karp}, |
198 | 208 |
title = {A characterization of the minimum cycle mean in a |
199 | 209 |
digraph}, |
200 | 210 |
journal = {Discrete Math.}, |
201 | 211 |
year = 1978, |
202 | 212 |
volume = 23, |
203 | 213 |
pages = {309-311} |
204 | 214 |
} |
205 | 215 |
|
206 | 216 |
@article{dasdan98minmeancycle, |
207 | 217 |
author = {Ali Dasdan and Rajesh K. Gupta}, |
208 | 218 |
title = {Faster Maximum and Minimum Mean Cycle Alogrithms for |
209 | 219 |
System Performance Analysis}, |
210 | 220 |
journal = {IEEE Transactions on Computer-Aided Design of |
211 | 221 |
Integrated Circuits and Systems}, |
212 | 222 |
year = 1998, |
213 | 223 |
volume = 17, |
214 | 224 |
number = 10, |
215 | 225 |
pages = {889-899} |
216 | 226 |
} |
217 | 227 |
|
218 | 228 |
|
219 | 229 |
%%%%% Minimum cost flow algorithms %%%%% |
220 | 230 |
|
221 | 231 |
@article{klein67primal, |
222 | 232 |
author = {Morton Klein}, |
223 | 233 |
title = {A primal method for minimal cost flows with |
224 | 234 |
applications to the assignment and transportation |
225 | 235 |
problems}, |
226 | 236 |
journal = {Management Science}, |
227 | 237 |
year = 1967, |
228 | 238 |
volume = 14, |
229 | 239 |
pages = {205-220} |
230 | 240 |
} |
231 | 241 |
|
232 |
@ |
|
242 |
@article{goldberg89cyclecanceling, |
|
233 | 243 |
author = {Andrew V. Goldberg and Robert E. Tarjan}, |
234 | 244 |
title = {Finding minimum-cost circulations by canceling |
235 | 245 |
negative cycles}, |
236 |
booktitle = {STOC '88: Proceedings of the Twentieth Annual ACM |
|
237 |
Symposium on Theory of Computing}, |
|
238 |
year = 1988, |
|
239 |
publisher = {ACM Press}, |
|
240 |
address = {New York, NY}, |
|
241 |
pages = {388-397} |
|
246 |
journal = {Journal of the ACM}, |
|
247 |
year = 1989, |
|
248 |
volume = 36, |
|
249 |
number = 4, |
|
250 |
pages = {873-886} |
|
242 | 251 |
} |
243 | 252 |
|
244 |
@article{edmondskarp72theoretical, |
|
245 |
author = {Jack Edmonds and Richard M. Karp}, |
|
246 |
title = {Theoretical improvements in algorithmic efficiency |
|
247 |
for network flow problems}, |
|
248 |
journal = {Journal of the ACM}, |
|
249 |
year = 1972, |
|
250 |
volume = 19, |
|
251 |
number = 2, |
|
252 |
pages = {248-264} |
|
253 |
} |
|
254 |
|
|
255 |
@inproceedings{goldberg87approximation, |
|
256 |
author = {Andrew V. Goldberg and Robert E. Tarjan}, |
|
257 |
title = {Solving minimum-cost flow problems by successive |
|
258 |
approximation}, |
|
259 |
booktitle = {STOC '87: Proceedings of the Nineteenth Annual ACM |
|
260 |
Symposium on Theory of Computing}, |
|
261 |
year = 1987, |
|
262 |
publisher = {ACM Press}, |
|
263 |
address = {New York, NY}, |
|
264 |
pages = {7-18} |
|
265 |
} |
|
266 |
|
|
267 |
@article{goldberg90finding, |
|
253 |
@article{goldberg90approximation, |
|
268 | 254 |
author = {Andrew V. Goldberg and Robert E. Tarjan}, |
269 | 255 |
title = {Finding Minimum-Cost Circulations by Successive |
270 | 256 |
Approximation}, |
271 | 257 |
journal = {Mathematics of Operations Research}, |
272 | 258 |
year = 1990, |
273 | 259 |
volume = 15, |
274 | 260 |
number = 3, |
275 | 261 |
pages = {430-466} |
276 | 262 |
} |
277 | 263 |
|
278 | 264 |
@article{goldberg97efficient, |
279 | 265 |
author = {Andrew V. Goldberg}, |
280 | 266 |
title = {An Efficient Implementation of a Scaling |
281 | 267 |
Minimum-Cost Flow Algorithm}, |
282 | 268 |
journal = {Journal of Algorithms}, |
283 | 269 |
year = 1997, |
284 | 270 |
volume = 22, |
285 | 271 |
number = 1, |
286 | 272 |
pages = {1-29} |
287 | 273 |
} |
288 | 274 |
|
289 | 275 |
@article{bunnagel98efficient, |
290 | 276 |
author = {Ursula B{\"u}nnagel and Bernhard Korte and Jens |
291 | 277 |
Vygen}, |
292 | 278 |
title = {Efficient implementation of the {G}oldberg-{T}arjan |
293 | 279 |
minimum-cost flow algorithm}, |
294 | 280 |
journal = {Optimization Methods and Software}, |
295 | 281 |
year = 1998, |
296 | 282 |
volume = 10, |
297 | 283 |
pages = {157-174} |
298 | 284 |
} |
299 | 285 |
|
286 |
@book{dantzig63linearprog, |
|
287 |
author = {George B. Dantzig}, |
|
288 |
title = {Linear Programming and Extensions}, |
|
289 |
publisher = {Princeton University Press}, |
|
290 |
year = 1963 |
|
291 |
} |
|
292 |
|
|
300 | 293 |
@mastersthesis{kellyoneill91netsimplex, |
301 | 294 |
author = {Damian J. Kelly and Garrett M. O'Neill}, |
302 | 295 |
title = {The Minimum Cost Flow Problem and The Network |
303 | 296 |
Simplex Method}, |
304 | 297 |
school = {University College}, |
305 | 298 |
address = {Dublin, Ireland}, |
306 | 299 |
year = 1991, |
307 | 300 |
month = sep, |
308 | 301 |
} |
309 |
|
|
310 |
@techreport{lobel96networksimplex, |
|
311 |
author = {Andreas L{\"o}bel}, |
|
312 |
title = {Solving large-scale real-world minimum-cost flow |
|
313 |
problems by a network simplex method}, |
|
314 |
institution = {Konrad-Zuse-Zentrum fur Informationstechnik Berlin |
|
315 |
({ZIB})}, |
|
316 |
address = {Berlin, Germany}, |
|
317 |
year = 1996, |
|
318 |
number = {SC 96-7} |
|
319 |
} |
|
320 |
|
|
321 |
@article{frangioni06computational, |
|
322 |
author = {Antonio Frangioni and Antonio Manca}, |
|
323 |
title = {A Computational Study of Cost Reoptimization for |
|
324 |
Min-Cost Flow Problems}, |
|
325 |
journal = {INFORMS Journal On Computing}, |
|
326 |
year = 2006, |
|
327 |
volume = 18, |
|
328 |
number = 1, |
|
329 |
pages = {61-70} |
|
330 |
} |
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 | 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 | 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 |
/// for finding a \ref min_cost_flow "minimum cost flow" |
|
43 |
/// for finding a \ref min_cost_flow "minimum cost flow" |
|
44 |
/// \ref amo93networkflows, \ref dantzig63linearprog, |
|
45 |
/// \ref kellyoneill91netsimplex. |
|
44 | 46 |
/// This algorithm is a specialized version of the linear programming |
45 | 47 |
/// simplex method directly for the minimum cost flow problem. |
46 | 48 |
/// It is one of the most efficient solution methods. |
47 | 49 |
/// |
48 | 50 |
/// In general this class is the fastest implementation available |
49 | 51 |
/// in LEMON for the minimum cost flow problem. |
50 | 52 |
/// Moreover it supports both directions of the supply/demand inequality |
51 | 53 |
/// constraints. For more information see \ref SupplyType. |
52 | 54 |
/// |
53 | 55 |
/// Most of the parameters of the problem (except for the digraph) |
54 | 56 |
/// can be given using separate functions, and the algorithm can be |
55 | 57 |
/// executed using the \ref run() function. If some parameters are not |
56 | 58 |
/// specified, then default values will be used. |
57 | 59 |
/// |
58 | 60 |
/// \tparam GR The digraph type the algorithm runs on. |
59 | 61 |
/// \tparam V The value type used for flow amounts, capacity bounds |
60 | 62 |
/// and supply values in the algorithm. By default it is \c int. |
61 | 63 |
/// \tparam C The value type used for costs and potentials in the |
62 | 64 |
/// algorithm. By default it is the same as \c V. |
63 | 65 |
/// |
64 | 66 |
/// \warning Both value types must be signed and all input data must |
65 | 67 |
/// be integer. |
66 | 68 |
/// |
67 | 69 |
/// \note %NetworkSimplex provides five different pivot rule |
68 | 70 |
/// implementations, from which the most efficient one is used |
69 | 71 |
/// by default. For more information see \ref PivotRule. |
70 | 72 |
template <typename GR, typename V = int, typename C = V> |
71 | 73 |
class NetworkSimplex |
72 | 74 |
{ |
73 | 75 |
public: |
74 | 76 |
|
75 | 77 |
/// The type of the flow amounts, capacity bounds and supply values |
76 | 78 |
typedef V Value; |
77 | 79 |
/// The type of the arc costs |
78 | 80 |
typedef C Cost; |
79 | 81 |
|
80 | 82 |
public: |
81 | 83 |
|
82 | 84 |
/// \brief Problem type constants for the \c run() function. |
83 | 85 |
/// |
84 | 86 |
/// Enum type containing the problem type constants that can be |
85 | 87 |
/// returned by the \ref run() function of the algorithm. |
86 | 88 |
enum ProblemType { |
87 | 89 |
/// The problem has no feasible solution (flow). |
88 | 90 |
INFEASIBLE, |
89 | 91 |
/// The problem has optimal solution (i.e. it is feasible and |
90 | 92 |
/// bounded), and the algorithm has found optimal flow and node |
91 | 93 |
/// potentials (primal and dual solutions). |
92 | 94 |
OPTIMAL, |
93 | 95 |
/// The objective function of the problem is unbounded, i.e. |
94 | 96 |
/// there is a directed cycle having negative total cost and |
95 | 97 |
/// infinite upper bound. |
96 | 98 |
UNBOUNDED |
97 | 99 |
}; |
98 | 100 |
|
99 | 101 |
/// \brief Constants for selecting the type of the supply constraints. |
100 | 102 |
/// |
101 | 103 |
/// Enum type containing constants for selecting the supply type, |
102 | 104 |
/// i.e. the direction of the inequalities in the supply/demand |
103 | 105 |
/// constraints of the \ref min_cost_flow "minimum cost flow problem". |
104 | 106 |
/// |
105 | 107 |
/// The default supply type is \c GEQ, the \c LEQ type can be |
106 | 108 |
/// selected using \ref supplyType(). |
107 | 109 |
/// The equality form is a special case of both supply types. |
108 | 110 |
enum SupplyType { |
109 | 111 |
/// This option means that there are <em>"greater or equal"</em> |
110 | 112 |
/// supply/demand constraints in the definition of the problem. |
111 | 113 |
GEQ, |
112 | 114 |
/// This option means that there are <em>"less or equal"</em> |
113 | 115 |
/// supply/demand constraints in the definition of the problem. |
114 | 116 |
LEQ |
115 | 117 |
}; |
116 | 118 |
|
117 | 119 |
/// \brief Constants for selecting the pivot rule. |
118 | 120 |
/// |
119 | 121 |
/// Enum type containing constants for selecting the pivot rule for |
120 | 122 |
/// the \ref run() function. |
121 | 123 |
/// |
122 | 124 |
/// \ref NetworkSimplex provides five different pivot rule |
123 | 125 |
/// implementations that significantly affect the running time |
124 | 126 |
/// of the algorithm. |
125 | 127 |
/// By default \ref BLOCK_SEARCH "Block Search" is used, which |
126 | 128 |
/// proved to be the most efficient and the most robust on various |
127 | 129 |
/// test inputs according to our benchmark tests. |
128 | 130 |
/// However another pivot rule can be selected using the \ref run() |
129 | 131 |
/// function with the proper parameter. |
130 | 132 |
enum PivotRule { |
131 | 133 |
|
132 | 134 |
/// The First Eligible pivot rule. |
133 | 135 |
/// The next eligible arc is selected in a wraparound fashion |
134 | 136 |
/// in every iteration. |
135 | 137 |
FIRST_ELIGIBLE, |
136 | 138 |
|
137 | 139 |
/// The Best Eligible pivot rule. |
138 | 140 |
/// The best eligible arc is selected in every iteration. |
139 | 141 |
BEST_ELIGIBLE, |
140 | 142 |
|
141 | 143 |
/// The Block Search pivot rule. |
142 | 144 |
/// A specified number of arcs are examined in every iteration |
143 | 145 |
/// in a wraparound fashion and the best eligible arc is selected |
144 | 146 |
/// from this block. |
145 | 147 |
BLOCK_SEARCH, |
146 | 148 |
|
147 | 149 |
/// The Candidate List pivot rule. |
148 | 150 |
/// In a major iteration a candidate list is built from eligible arcs |
149 | 151 |
/// in a wraparound fashion and in the following minor iterations |
150 | 152 |
/// the best eligible arc is selected from this list. |
151 | 153 |
CANDIDATE_LIST, |
152 | 154 |
|
153 | 155 |
/// The Altering Candidate List pivot rule. |
154 | 156 |
/// It is a modified version of the Candidate List method. |
155 | 157 |
/// It keeps only the several best eligible arcs from the former |
156 | 158 |
/// candidate list and extends this list in every iteration. |
157 | 159 |
ALTERING_LIST |
158 | 160 |
}; |
159 | 161 |
|
160 | 162 |
private: |
161 | 163 |
|
162 | 164 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
163 | 165 |
|
164 | 166 |
typedef std::vector<int> IntVector; |
165 | 167 |
typedef std::vector<bool> BoolVector; |
166 | 168 |
typedef std::vector<Value> ValueVector; |
167 | 169 |
typedef std::vector<Cost> CostVector; |
168 | 170 |
|
169 | 171 |
// State constants for arcs |
170 | 172 |
enum ArcStateEnum { |
171 | 173 |
STATE_UPPER = -1, |
172 | 174 |
STATE_TREE = 0, |
173 | 175 |
STATE_LOWER = 1 |
174 | 176 |
}; |
175 | 177 |
|
176 | 178 |
private: |
177 | 179 |
|
178 | 180 |
// Data related to the underlying digraph |
179 | 181 |
const GR &_graph; |
180 | 182 |
int _node_num; |
181 | 183 |
int _arc_num; |
182 | 184 |
int _all_arc_num; |
183 | 185 |
int _search_arc_num; |
184 | 186 |
|
185 | 187 |
// Parameters of the problem |
186 | 188 |
bool _have_lower; |
187 | 189 |
SupplyType _stype; |
188 | 190 |
Value _sum_supply; |
189 | 191 |
|
190 | 192 |
// Data structures for storing the digraph |
191 | 193 |
IntNodeMap _node_id; |
192 | 194 |
IntArcMap _arc_id; |
193 | 195 |
IntVector _source; |
194 | 196 |
IntVector _target; |
195 | 197 |
|
196 | 198 |
// Node and arc data |
197 | 199 |
ValueVector _lower; |
198 | 200 |
ValueVector _upper; |
199 | 201 |
ValueVector _cap; |
200 | 202 |
CostVector _cost; |
201 | 203 |
ValueVector _supply; |
202 | 204 |
ValueVector _flow; |
203 | 205 |
CostVector _pi; |
204 | 206 |
|
205 | 207 |
// Data for storing the spanning tree structure |
206 | 208 |
IntVector _parent; |
207 | 209 |
IntVector _pred; |
208 | 210 |
IntVector _thread; |
209 | 211 |
IntVector _rev_thread; |
210 | 212 |
IntVector _succ_num; |
211 | 213 |
IntVector _last_succ; |
212 | 214 |
IntVector _dirty_revs; |
213 | 215 |
BoolVector _forward; |
214 | 216 |
IntVector _state; |
215 | 217 |
int _root; |
216 | 218 |
|
217 | 219 |
// Temporary data used in the current pivot iteration |
218 | 220 |
int in_arc, join, u_in, v_in, u_out, v_out; |
219 | 221 |
int first, second, right, last; |
220 | 222 |
int stem, par_stem, new_stem; |
221 | 223 |
Value delta; |
222 | 224 |
|
223 | 225 |
public: |
224 | 226 |
|
225 | 227 |
/// \brief Constant for infinite upper bounds (capacities). |
226 | 228 |
/// |
227 | 229 |
/// Constant for infinite upper bounds (capacities). |
228 | 230 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
229 | 231 |
/// \c std::numeric_limits<Value>::max() otherwise. |
230 | 232 |
const Value INF; |
231 | 233 |
|
232 | 234 |
private: |
233 | 235 |
|
234 | 236 |
// Implementation of the First Eligible pivot rule |
235 | 237 |
class FirstEligiblePivotRule |
236 | 238 |
{ |
237 | 239 |
private: |
238 | 240 |
|
239 | 241 |
// References to the NetworkSimplex class |
240 | 242 |
const IntVector &_source; |
241 | 243 |
const IntVector &_target; |
242 | 244 |
const CostVector &_cost; |
243 | 245 |
const IntVector &_state; |
244 | 246 |
const CostVector &_pi; |
245 | 247 |
int &_in_arc; |
246 | 248 |
int _search_arc_num; |
247 | 249 |
|
248 | 250 |
// Pivot rule data |
249 | 251 |
int _next_arc; |
250 | 252 |
|
251 | 253 |
public: |
252 | 254 |
|
253 | 255 |
// Constructor |
254 | 256 |
FirstEligiblePivotRule(NetworkSimplex &ns) : |
255 | 257 |
_source(ns._source), _target(ns._target), |
256 | 258 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
257 | 259 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
258 | 260 |
_next_arc(0) |
259 | 261 |
{} |
260 | 262 |
|
261 | 263 |
// Find next entering arc |
262 | 264 |
bool findEnteringArc() { |
263 | 265 |
Cost c; |
264 | 266 |
for (int e = _next_arc; e < _search_arc_num; ++e) { |
265 | 267 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
266 | 268 |
if (c < 0) { |
267 | 269 |
_in_arc = e; |
268 | 270 |
_next_arc = e + 1; |
269 | 271 |
return true; |
270 | 272 |
} |
271 | 273 |
} |
272 | 274 |
for (int e = 0; e < _next_arc; ++e) { |
273 | 275 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
274 | 276 |
if (c < 0) { |
275 | 277 |
_in_arc = e; |
276 | 278 |
_next_arc = e + 1; |
277 | 279 |
return true; |
278 | 280 |
} |
279 | 281 |
} |
280 | 282 |
return false; |
281 | 283 |
} |
282 | 284 |
|
283 | 285 |
}; //class FirstEligiblePivotRule |
284 | 286 |
|
285 | 287 |
|
286 | 288 |
// Implementation of the Best Eligible pivot rule |
287 | 289 |
class BestEligiblePivotRule |
288 | 290 |
{ |
289 | 291 |
private: |
290 | 292 |
|
291 | 293 |
// References to the NetworkSimplex class |
292 | 294 |
const IntVector &_source; |
293 | 295 |
const IntVector &_target; |
294 | 296 |
const CostVector &_cost; |
295 | 297 |
const IntVector &_state; |
296 | 298 |
const CostVector &_pi; |
297 | 299 |
int &_in_arc; |
298 | 300 |
int _search_arc_num; |
299 | 301 |
|
300 | 302 |
public: |
301 | 303 |
|
302 | 304 |
// Constructor |
303 | 305 |
BestEligiblePivotRule(NetworkSimplex &ns) : |
304 | 306 |
_source(ns._source), _target(ns._target), |
305 | 307 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
306 | 308 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) |
307 | 309 |
{} |
308 | 310 |
|
309 | 311 |
// Find next entering arc |
310 | 312 |
bool findEnteringArc() { |
311 | 313 |
Cost c, min = 0; |
312 | 314 |
for (int e = 0; e < _search_arc_num; ++e) { |
313 | 315 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
314 | 316 |
if (c < min) { |
315 | 317 |
min = c; |
316 | 318 |
_in_arc = e; |
317 | 319 |
} |
318 | 320 |
} |
319 | 321 |
return min < 0; |
320 | 322 |
} |
321 | 323 |
|
322 | 324 |
}; //class BestEligiblePivotRule |
323 | 325 |
|
324 | 326 |
|
325 | 327 |
// Implementation of the Block Search pivot rule |
326 | 328 |
class BlockSearchPivotRule |
327 | 329 |
{ |
328 | 330 |
private: |
329 | 331 |
|
330 | 332 |
// References to the NetworkSimplex class |
331 | 333 |
const IntVector &_source; |
332 | 334 |
const IntVector &_target; |
333 | 335 |
const CostVector &_cost; |
334 | 336 |
const IntVector &_state; |
335 | 337 |
const CostVector &_pi; |
336 | 338 |
int &_in_arc; |
337 | 339 |
int _search_arc_num; |
338 | 340 |
|
339 | 341 |
// Pivot rule data |
340 | 342 |
int _block_size; |
341 | 343 |
int _next_arc; |
342 | 344 |
|
343 | 345 |
public: |
344 | 346 |
|
345 | 347 |
// Constructor |
346 | 348 |
BlockSearchPivotRule(NetworkSimplex &ns) : |
347 | 349 |
_source(ns._source), _target(ns._target), |
348 | 350 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
349 | 351 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
350 | 352 |
_next_arc(0) |
351 | 353 |
{ |
352 | 354 |
// The main parameters of the pivot rule |
353 | 355 |
const double BLOCK_SIZE_FACTOR = 0.5; |
354 | 356 |
const int MIN_BLOCK_SIZE = 10; |
355 | 357 |
|
356 | 358 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * |
357 | 359 |
std::sqrt(double(_search_arc_num))), |
358 | 360 |
MIN_BLOCK_SIZE ); |
359 | 361 |
} |
360 | 362 |
|
361 | 363 |
// Find next entering arc |
362 | 364 |
bool findEnteringArc() { |
363 | 365 |
Cost c, min = 0; |
364 | 366 |
int cnt = _block_size; |
365 | 367 |
int e; |
366 | 368 |
for (e = _next_arc; e < _search_arc_num; ++e) { |
367 | 369 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
368 | 370 |
if (c < min) { |
369 | 371 |
min = c; |
370 | 372 |
_in_arc = e; |
371 | 373 |
} |
372 | 374 |
if (--cnt == 0) { |
373 | 375 |
if (min < 0) goto search_end; |
374 | 376 |
cnt = _block_size; |
375 | 377 |
} |
376 | 378 |
} |
377 | 379 |
for (e = 0; e < _next_arc; ++e) { |
378 | 380 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
379 | 381 |
if (c < min) { |
380 | 382 |
min = c; |
381 | 383 |
_in_arc = e; |
382 | 384 |
} |
383 | 385 |
if (--cnt == 0) { |
384 | 386 |
if (min < 0) goto search_end; |
385 | 387 |
cnt = _block_size; |
386 | 388 |
} |
387 | 389 |
} |
388 | 390 |
if (min >= 0) return false; |
389 | 391 |
|
390 | 392 |
search_end: |
391 | 393 |
_next_arc = e; |
392 | 394 |
return true; |
393 | 395 |
} |
394 | 396 |
|
395 | 397 |
}; //class BlockSearchPivotRule |
396 | 398 |
|
397 | 399 |
|
398 | 400 |
// Implementation of the Candidate List pivot rule |
399 | 401 |
class CandidateListPivotRule |
400 | 402 |
{ |
401 | 403 |
private: |
402 | 404 |
|
403 | 405 |
// References to the NetworkSimplex class |
404 | 406 |
const IntVector &_source; |
405 | 407 |
const IntVector &_target; |
406 | 408 |
const CostVector &_cost; |
407 | 409 |
const IntVector &_state; |
408 | 410 |
const CostVector &_pi; |
409 | 411 |
int &_in_arc; |
410 | 412 |
int _search_arc_num; |
411 | 413 |
|
412 | 414 |
// Pivot rule data |
413 | 415 |
IntVector _candidates; |
414 | 416 |
int _list_length, _minor_limit; |
415 | 417 |
int _curr_length, _minor_count; |
416 | 418 |
int _next_arc; |
417 | 419 |
|
418 | 420 |
public: |
419 | 421 |
|
420 | 422 |
/// Constructor |
421 | 423 |
CandidateListPivotRule(NetworkSimplex &ns) : |
422 | 424 |
_source(ns._source), _target(ns._target), |
423 | 425 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
424 | 426 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
425 | 427 |
_next_arc(0) |
426 | 428 |
{ |
427 | 429 |
// The main parameters of the pivot rule |
428 | 430 |
const double LIST_LENGTH_FACTOR = 0.25; |
429 | 431 |
const int MIN_LIST_LENGTH = 10; |
430 | 432 |
const double MINOR_LIMIT_FACTOR = 0.1; |
431 | 433 |
const int MIN_MINOR_LIMIT = 3; |
432 | 434 |
|
433 | 435 |
_list_length = std::max( int(LIST_LENGTH_FACTOR * |
434 | 436 |
std::sqrt(double(_search_arc_num))), |
435 | 437 |
MIN_LIST_LENGTH ); |
436 | 438 |
_minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), |
437 | 439 |
MIN_MINOR_LIMIT ); |
438 | 440 |
_curr_length = _minor_count = 0; |
439 | 441 |
_candidates.resize(_list_length); |
440 | 442 |
} |
441 | 443 |
|
442 | 444 |
/// Find next entering arc |
443 | 445 |
bool findEnteringArc() { |
444 | 446 |
Cost min, c; |
445 | 447 |
int e; |
446 | 448 |
if (_curr_length > 0 && _minor_count < _minor_limit) { |
447 | 449 |
// Minor iteration: select the best eligible arc from the |
448 | 450 |
// current candidate list |
449 | 451 |
++_minor_count; |
450 | 452 |
min = 0; |
451 | 453 |
for (int i = 0; i < _curr_length; ++i) { |
452 | 454 |
e = _candidates[i]; |
453 | 455 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
454 | 456 |
if (c < min) { |
455 | 457 |
min = c; |
456 | 458 |
_in_arc = e; |
457 | 459 |
} |
458 | 460 |
else if (c >= 0) { |
459 | 461 |
_candidates[i--] = _candidates[--_curr_length]; |
460 | 462 |
} |
461 | 463 |
} |
462 | 464 |
if (min < 0) return true; |
463 | 465 |
} |
464 | 466 |
|
465 | 467 |
// Major iteration: build a new candidate list |
466 | 468 |
min = 0; |
467 | 469 |
_curr_length = 0; |
468 | 470 |
for (e = _next_arc; e < _search_arc_num; ++e) { |
469 | 471 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
470 | 472 |
if (c < 0) { |
471 | 473 |
_candidates[_curr_length++] = e; |
472 | 474 |
if (c < min) { |
473 | 475 |
min = c; |
474 | 476 |
_in_arc = e; |
475 | 477 |
} |
476 | 478 |
if (_curr_length == _list_length) goto search_end; |
477 | 479 |
} |
478 | 480 |
} |
479 | 481 |
for (e = 0; e < _next_arc; ++e) { |
480 | 482 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
481 | 483 |
if (c < 0) { |
482 | 484 |
_candidates[_curr_length++] = e; |
483 | 485 |
if (c < min) { |
484 | 486 |
min = c; |
485 | 487 |
_in_arc = e; |
486 | 488 |
} |
487 | 489 |
if (_curr_length == _list_length) goto search_end; |
488 | 490 |
} |
489 | 491 |
} |
490 | 492 |
if (_curr_length == 0) return false; |
491 | 493 |
|
492 | 494 |
search_end: |
493 | 495 |
_minor_count = 1; |
494 | 496 |
_next_arc = e; |
495 | 497 |
return true; |
496 | 498 |
} |
497 | 499 |
|
498 | 500 |
}; //class CandidateListPivotRule |
499 | 501 |
|
500 | 502 |
|
501 | 503 |
// Implementation of the Altering Candidate List pivot rule |
502 | 504 |
class AlteringListPivotRule |
503 | 505 |
{ |
504 | 506 |
private: |
505 | 507 |
|
506 | 508 |
// References to the NetworkSimplex class |
507 | 509 |
const IntVector &_source; |
508 | 510 |
const IntVector &_target; |
509 | 511 |
const CostVector &_cost; |
510 | 512 |
const IntVector &_state; |
511 | 513 |
const CostVector &_pi; |
512 | 514 |
int &_in_arc; |
513 | 515 |
int _search_arc_num; |
514 | 516 |
|
515 | 517 |
// Pivot rule data |
516 | 518 |
int _block_size, _head_length, _curr_length; |
517 | 519 |
int _next_arc; |
518 | 520 |
IntVector _candidates; |
519 | 521 |
CostVector _cand_cost; |
520 | 522 |
|
521 | 523 |
// Functor class to compare arcs during sort of the candidate list |
522 | 524 |
class SortFunc |
523 | 525 |
{ |
524 | 526 |
private: |
525 | 527 |
const CostVector &_map; |
526 | 528 |
public: |
527 | 529 |
SortFunc(const CostVector &map) : _map(map) {} |
528 | 530 |
bool operator()(int left, int right) { |
529 | 531 |
return _map[left] > _map[right]; |
530 | 532 |
} |
531 | 533 |
}; |
532 | 534 |
|
533 | 535 |
SortFunc _sort_func; |
534 | 536 |
|
535 | 537 |
public: |
536 | 538 |
|
537 | 539 |
// Constructor |
538 | 540 |
AlteringListPivotRule(NetworkSimplex &ns) : |
539 | 541 |
_source(ns._source), _target(ns._target), |
540 | 542 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
541 | 543 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
542 | 544 |
_next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) |
543 | 545 |
{ |
544 | 546 |
// The main parameters of the pivot rule |
545 | 547 |
const double BLOCK_SIZE_FACTOR = 1.0; |
546 | 548 |
const int MIN_BLOCK_SIZE = 10; |
547 | 549 |
const double HEAD_LENGTH_FACTOR = 0.1; |
548 | 550 |
const int MIN_HEAD_LENGTH = 3; |
549 | 551 |
|
550 | 552 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * |
551 | 553 |
std::sqrt(double(_search_arc_num))), |
552 | 554 |
MIN_BLOCK_SIZE ); |
553 | 555 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), |
554 | 556 |
MIN_HEAD_LENGTH ); |
555 | 557 |
_candidates.resize(_head_length + _block_size); |
556 | 558 |
_curr_length = 0; |
557 | 559 |
} |
558 | 560 |
|
559 | 561 |
// Find next entering arc |
560 | 562 |
bool findEnteringArc() { |
561 | 563 |
// Check the current candidate list |
562 | 564 |
int e; |
563 | 565 |
for (int i = 0; i < _curr_length; ++i) { |
564 | 566 |
e = _candidates[i]; |
565 | 567 |
_cand_cost[e] = _state[e] * |
566 | 568 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
567 | 569 |
if (_cand_cost[e] >= 0) { |
568 | 570 |
_candidates[i--] = _candidates[--_curr_length]; |
569 | 571 |
} |
570 | 572 |
} |
571 | 573 |
|
572 | 574 |
// Extend the list |
573 | 575 |
int cnt = _block_size; |
574 | 576 |
int limit = _head_length; |
575 | 577 |
|
576 | 578 |
for (e = _next_arc; e < _search_arc_num; ++e) { |
577 | 579 |
_cand_cost[e] = _state[e] * |
578 | 580 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
579 | 581 |
if (_cand_cost[e] < 0) { |
580 | 582 |
_candidates[_curr_length++] = e; |
581 | 583 |
} |
582 | 584 |
if (--cnt == 0) { |
583 | 585 |
if (_curr_length > limit) goto search_end; |
584 | 586 |
limit = 0; |
585 | 587 |
cnt = _block_size; |
586 | 588 |
} |
587 | 589 |
} |
588 | 590 |
for (e = 0; e < _next_arc; ++e) { |
589 | 591 |
_cand_cost[e] = _state[e] * |
590 | 592 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
591 | 593 |
if (_cand_cost[e] < 0) { |
592 | 594 |
_candidates[_curr_length++] = e; |
593 | 595 |
} |
594 | 596 |
if (--cnt == 0) { |
595 | 597 |
if (_curr_length > limit) goto search_end; |
596 | 598 |
limit = 0; |
597 | 599 |
cnt = _block_size; |
598 | 600 |
} |
599 | 601 |
} |
600 | 602 |
if (_curr_length == 0) return false; |
601 | 603 |
|
602 | 604 |
search_end: |
603 | 605 |
|
604 | 606 |
// Make heap of the candidate list (approximating a partial sort) |
605 | 607 |
make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
606 | 608 |
_sort_func ); |
607 | 609 |
|
608 | 610 |
// Pop the first element of the heap |
609 | 611 |
_in_arc = _candidates[0]; |
610 | 612 |
_next_arc = e; |
611 | 613 |
pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
612 | 614 |
_sort_func ); |
613 | 615 |
_curr_length = std::min(_head_length, _curr_length - 1); |
614 | 616 |
return true; |
615 | 617 |
} |
616 | 618 |
|
617 | 619 |
}; //class AlteringListPivotRule |
618 | 620 |
|
619 | 621 |
public: |
620 | 622 |
|
621 | 623 |
/// \brief Constructor. |
622 | 624 |
/// |
623 | 625 |
/// The constructor of the class. |
624 | 626 |
/// |
625 | 627 |
/// \param graph The digraph the algorithm runs on. |
626 | 628 |
/// \param arc_mixing Indicate if the arcs have to be stored in a |
627 | 629 |
/// mixed order in the internal data structure. |
628 | 630 |
/// In special cases, it could lead to better overall performance, |
629 | 631 |
/// but it is usually slower. Therefore it is disabled by default. |
630 | 632 |
NetworkSimplex(const GR& graph, bool arc_mixing = false) : |
631 | 633 |
_graph(graph), _node_id(graph), _arc_id(graph), |
632 | 634 |
INF(std::numeric_limits<Value>::has_infinity ? |
633 | 635 |
std::numeric_limits<Value>::infinity() : |
634 | 636 |
std::numeric_limits<Value>::max()) |
635 | 637 |
{ |
636 | 638 |
// Check the value types |
637 | 639 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
638 | 640 |
"The flow type of NetworkSimplex must be signed"); |
639 | 641 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
640 | 642 |
"The cost type of NetworkSimplex must be signed"); |
641 | 643 |
|
642 | 644 |
// Resize vectors |
643 | 645 |
_node_num = countNodes(_graph); |
644 | 646 |
_arc_num = countArcs(_graph); |
645 | 647 |
int all_node_num = _node_num + 1; |
646 | 648 |
int max_arc_num = _arc_num + 2 * _node_num; |
647 | 649 |
|
648 | 650 |
_source.resize(max_arc_num); |
649 | 651 |
_target.resize(max_arc_num); |
650 | 652 |
|
651 | 653 |
_lower.resize(_arc_num); |
652 | 654 |
_upper.resize(_arc_num); |
653 | 655 |
_cap.resize(max_arc_num); |
654 | 656 |
_cost.resize(max_arc_num); |
655 | 657 |
_supply.resize(all_node_num); |
656 | 658 |
_flow.resize(max_arc_num); |
657 | 659 |
_pi.resize(all_node_num); |
658 | 660 |
|
659 | 661 |
_parent.resize(all_node_num); |
660 | 662 |
_pred.resize(all_node_num); |
661 | 663 |
_forward.resize(all_node_num); |
662 | 664 |
_thread.resize(all_node_num); |
663 | 665 |
_rev_thread.resize(all_node_num); |
664 | 666 |
_succ_num.resize(all_node_num); |
665 | 667 |
_last_succ.resize(all_node_num); |
666 | 668 |
_state.resize(max_arc_num); |
667 | 669 |
|
668 | 670 |
// Copy the graph |
669 | 671 |
int i = 0; |
670 | 672 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
671 | 673 |
_node_id[n] = i; |
672 | 674 |
} |
673 | 675 |
if (arc_mixing) { |
674 | 676 |
// Store the arcs in a mixed order |
675 | 677 |
int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
676 | 678 |
int i = 0, j = 0; |
677 | 679 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
678 | 680 |
_arc_id[a] = i; |
679 | 681 |
_source[i] = _node_id[_graph.source(a)]; |
680 | 682 |
_target[i] = _node_id[_graph.target(a)]; |
681 | 683 |
if ((i += k) >= _arc_num) i = ++j; |
682 | 684 |
} |
683 | 685 |
} else { |
684 | 686 |
// Store the arcs in the original order |
685 | 687 |
int i = 0; |
686 | 688 |
for (ArcIt a(_graph); a != INVALID; ++a, ++i) { |
687 | 689 |
_arc_id[a] = i; |
688 | 690 |
_source[i] = _node_id[_graph.source(a)]; |
689 | 691 |
_target[i] = _node_id[_graph.target(a)]; |
690 | 692 |
} |
691 | 693 |
} |
692 | 694 |
|
693 | 695 |
// Reset parameters |
694 | 696 |
reset(); |
695 | 697 |
} |
696 | 698 |
|
697 | 699 |
/// \name Parameters |
698 | 700 |
/// The parameters of the algorithm can be specified using these |
699 | 701 |
/// functions. |
700 | 702 |
|
701 | 703 |
/// @{ |
702 | 704 |
|
703 | 705 |
/// \brief Set the lower bounds on the arcs. |
704 | 706 |
/// |
705 | 707 |
/// This function sets the lower bounds on the arcs. |
706 | 708 |
/// If it is not used before calling \ref run(), the lower bounds |
707 | 709 |
/// will be set to zero on all arcs. |
708 | 710 |
/// |
709 | 711 |
/// \param map An arc map storing the lower bounds. |
710 | 712 |
/// Its \c Value type must be convertible to the \c Value type |
711 | 713 |
/// of the algorithm. |
712 | 714 |
/// |
713 | 715 |
/// \return <tt>(*this)</tt> |
714 | 716 |
template <typename LowerMap> |
715 | 717 |
NetworkSimplex& lowerMap(const LowerMap& map) { |
716 | 718 |
_have_lower = true; |
717 | 719 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
718 | 720 |
_lower[_arc_id[a]] = map[a]; |
719 | 721 |
} |
720 | 722 |
return *this; |
721 | 723 |
} |
722 | 724 |
|
723 | 725 |
/// \brief Set the upper bounds (capacities) on the arcs. |
724 | 726 |
/// |
725 | 727 |
/// This function sets the upper bounds (capacities) on the arcs. |
726 | 728 |
/// If it is not used before calling \ref run(), the upper bounds |
727 | 729 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
728 | 730 |
/// unbounded from above on each arc). |
729 | 731 |
/// |
730 | 732 |
/// \param map An arc map storing the upper bounds. |
731 | 733 |
/// Its \c Value type must be convertible to the \c Value type |
732 | 734 |
/// of the algorithm. |
733 | 735 |
/// |
734 | 736 |
/// \return <tt>(*this)</tt> |
735 | 737 |
template<typename UpperMap> |
736 | 738 |
NetworkSimplex& upperMap(const UpperMap& map) { |
737 | 739 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
738 | 740 |
_upper[_arc_id[a]] = map[a]; |
739 | 741 |
} |
740 | 742 |
return *this; |
741 | 743 |
} |
742 | 744 |
|
743 | 745 |
/// \brief Set the costs of the arcs. |
744 | 746 |
/// |
745 | 747 |
/// This function sets the costs of the arcs. |
746 | 748 |
/// If it is not used before calling \ref run(), the costs |
747 | 749 |
/// will be set to \c 1 on all arcs. |
748 | 750 |
/// |
749 | 751 |
/// \param map An arc map storing the costs. |
750 | 752 |
/// Its \c Value type must be convertible to the \c Cost type |
751 | 753 |
/// of the algorithm. |
752 | 754 |
/// |
753 | 755 |
/// \return <tt>(*this)</tt> |
754 | 756 |
template<typename CostMap> |
755 | 757 |
NetworkSimplex& costMap(const CostMap& map) { |
756 | 758 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
757 | 759 |
_cost[_arc_id[a]] = map[a]; |
758 | 760 |
} |
759 | 761 |
return *this; |
760 | 762 |
} |
761 | 763 |
|
762 | 764 |
/// \brief Set the supply values of the nodes. |
763 | 765 |
/// |
764 | 766 |
/// This function sets the supply values of the nodes. |
765 | 767 |
/// If neither this function nor \ref stSupply() is used before |
766 | 768 |
/// calling \ref run(), the supply of each node will be set to zero. |
767 | 769 |
/// |
768 | 770 |
/// \param map A node map storing the supply values. |
769 | 771 |
/// Its \c Value type must be convertible to the \c Value type |
770 | 772 |
/// of the algorithm. |
771 | 773 |
/// |
772 | 774 |
/// \return <tt>(*this)</tt> |
773 | 775 |
template<typename SupplyMap> |
774 | 776 |
NetworkSimplex& supplyMap(const SupplyMap& map) { |
775 | 777 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
776 | 778 |
_supply[_node_id[n]] = map[n]; |
777 | 779 |
} |
778 | 780 |
return *this; |
779 | 781 |
} |
780 | 782 |
|
781 | 783 |
/// \brief Set single source and target nodes and a supply value. |
782 | 784 |
/// |
783 | 785 |
/// This function sets a single source node and a single target node |
784 | 786 |
/// and the required flow value. |
785 | 787 |
/// If neither this function nor \ref supplyMap() is used before |
786 | 788 |
/// calling \ref run(), the supply of each node will be set to zero. |
787 | 789 |
/// |
788 | 790 |
/// Using this function has the same effect as using \ref supplyMap() |
789 | 791 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
790 | 792 |
/// assigned to \c t and all other nodes have zero supply value. |
791 | 793 |
/// |
792 | 794 |
/// \param s The source node. |
793 | 795 |
/// \param t The target node. |
794 | 796 |
/// \param k The required amount of flow from node \c s to node \c t |
795 | 797 |
/// (i.e. the supply of \c s and the demand of \c t). |
796 | 798 |
/// |
797 | 799 |
/// \return <tt>(*this)</tt> |
798 | 800 |
NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { |
799 | 801 |
for (int i = 0; i != _node_num; ++i) { |
800 | 802 |
_supply[i] = 0; |
801 | 803 |
} |
802 | 804 |
_supply[_node_id[s]] = k; |
803 | 805 |
_supply[_node_id[t]] = -k; |
804 | 806 |
return *this; |
805 | 807 |
} |
806 | 808 |
|
807 | 809 |
/// \brief Set the type of the supply constraints. |
808 | 810 |
/// |
809 | 811 |
/// This function sets the type of the supply/demand constraints. |
810 | 812 |
/// If it is not used before calling \ref run(), the \ref GEQ supply |
811 | 813 |
/// type will be used. |
812 | 814 |
/// |
813 | 815 |
/// For more information see \ref SupplyType. |
814 | 816 |
/// |
815 | 817 |
/// \return <tt>(*this)</tt> |
816 | 818 |
NetworkSimplex& supplyType(SupplyType supply_type) { |
817 | 819 |
_stype = supply_type; |
818 | 820 |
return *this; |
819 | 821 |
} |
820 | 822 |
|
821 | 823 |
/// @} |
822 | 824 |
|
823 | 825 |
/// \name Execution Control |
824 | 826 |
/// The algorithm can be executed using \ref run(). |
825 | 827 |
|
826 | 828 |
/// @{ |
827 | 829 |
|
828 | 830 |
/// \brief Run the algorithm. |
829 | 831 |
/// |
830 | 832 |
/// This function runs the algorithm. |
831 | 833 |
/// The paramters can be specified using functions \ref lowerMap(), |
832 | 834 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), |
833 | 835 |
/// \ref supplyType(). |
834 | 836 |
/// For example, |
835 | 837 |
/// \code |
836 | 838 |
/// NetworkSimplex<ListDigraph> ns(graph); |
837 | 839 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
838 | 840 |
/// .supplyMap(sup).run(); |
839 | 841 |
/// \endcode |
840 | 842 |
/// |
841 | 843 |
/// This function can be called more than once. All the parameters |
842 | 844 |
/// that have been given are kept for the next call, unless |
843 | 845 |
/// \ref reset() is called, thus only the modified parameters |
844 | 846 |
/// have to be set again. See \ref reset() for examples. |
845 | 847 |
/// However the underlying digraph must not be modified after this |
846 | 848 |
/// class have been constructed, since it copies and extends the graph. |
847 | 849 |
/// |
848 | 850 |
/// \param pivot_rule The pivot rule that will be used during the |
849 | 851 |
/// algorithm. For more information see \ref PivotRule. |
850 | 852 |
/// |
851 | 853 |
/// \return \c INFEASIBLE if no feasible flow exists, |
852 | 854 |
/// \n \c OPTIMAL if the problem has optimal solution |
853 | 855 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
854 | 856 |
/// optimal flow and node potentials (primal and dual solutions), |
855 | 857 |
/// \n \c UNBOUNDED if the objective function of the problem is |
856 | 858 |
/// unbounded, i.e. there is a directed cycle having negative total |
857 | 859 |
/// cost and infinite upper bound. |
858 | 860 |
/// |
859 | 861 |
/// \see ProblemType, PivotRule |
860 | 862 |
ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { |
861 | 863 |
if (!init()) return INFEASIBLE; |
862 | 864 |
return start(pivot_rule); |
863 | 865 |
} |
864 | 866 |
|
865 | 867 |
/// \brief Reset all the parameters that have been given before. |
866 | 868 |
/// |
867 | 869 |
/// This function resets all the paramaters that have been given |
868 | 870 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
869 | 871 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). |
870 | 872 |
/// |
871 | 873 |
/// It is useful for multiple run() calls. If this function is not |
872 | 874 |
/// used, all the parameters given before are kept for the next |
873 | 875 |
/// \ref run() call. |
874 | 876 |
/// However the underlying digraph must not be modified after this |
875 | 877 |
/// class have been constructed, since it copies and extends the graph. |
876 | 878 |
/// |
877 | 879 |
/// For example, |
878 | 880 |
/// \code |
879 | 881 |
/// NetworkSimplex<ListDigraph> ns(graph); |
880 | 882 |
/// |
881 | 883 |
/// // First run |
882 | 884 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
883 | 885 |
/// .supplyMap(sup).run(); |
884 | 886 |
/// |
885 | 887 |
/// // Run again with modified cost map (reset() is not called, |
886 | 888 |
/// // so only the cost map have to be set again) |
887 | 889 |
/// cost[e] += 100; |
888 | 890 |
/// ns.costMap(cost).run(); |
889 | 891 |
/// |
890 | 892 |
/// // Run again from scratch using reset() |
891 | 893 |
/// // (the lower bounds will be set to zero on all arcs) |
892 | 894 |
/// ns.reset(); |
893 | 895 |
/// ns.upperMap(capacity).costMap(cost) |
894 | 896 |
/// .supplyMap(sup).run(); |
895 | 897 |
/// \endcode |
896 | 898 |
/// |
897 | 899 |
/// \return <tt>(*this)</tt> |
898 | 900 |
NetworkSimplex& reset() { |
899 | 901 |
for (int i = 0; i != _node_num; ++i) { |
900 | 902 |
_supply[i] = 0; |
901 | 903 |
} |
902 | 904 |
for (int i = 0; i != _arc_num; ++i) { |
903 | 905 |
_lower[i] = 0; |
904 | 906 |
_upper[i] = INF; |
905 | 907 |
_cost[i] = 1; |
906 | 908 |
} |
907 | 909 |
_have_lower = false; |
908 | 910 |
_stype = GEQ; |
909 | 911 |
return *this; |
910 | 912 |
} |
911 | 913 |
|
912 | 914 |
/// @} |
913 | 915 |
|
914 | 916 |
/// \name Query Functions |
915 | 917 |
/// The results of the algorithm can be obtained using these |
916 | 918 |
/// functions.\n |
917 | 919 |
/// The \ref run() function must be called before using them. |
918 | 920 |
|
919 | 921 |
/// @{ |
920 | 922 |
|
921 | 923 |
/// \brief Return the total cost of the found flow. |
922 | 924 |
/// |
923 | 925 |
/// This function returns the total cost of the found flow. |
924 | 926 |
/// Its complexity is O(e). |
925 | 927 |
/// |
926 | 928 |
/// \note The return type of the function can be specified as a |
927 | 929 |
/// template parameter. For example, |
928 | 930 |
/// \code |
929 | 931 |
/// ns.totalCost<double>(); |
930 | 932 |
/// \endcode |
931 | 933 |
/// It is useful if the total cost cannot be stored in the \c Cost |
932 | 934 |
/// type of the algorithm, which is the default return type of the |
933 | 935 |
/// function. |
934 | 936 |
/// |
935 | 937 |
/// \pre \ref run() must be called before using this function. |
936 | 938 |
template <typename Number> |
937 | 939 |
Number totalCost() const { |
938 | 940 |
Number c = 0; |
939 | 941 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
940 | 942 |
int i = _arc_id[a]; |
941 | 943 |
c += Number(_flow[i]) * Number(_cost[i]); |
942 | 944 |
} |
943 | 945 |
return c; |
944 | 946 |
} |
945 | 947 |
|
946 | 948 |
#ifndef DOXYGEN |
947 | 949 |
Cost totalCost() const { |
948 | 950 |
return totalCost<Cost>(); |
949 | 951 |
} |
950 | 952 |
#endif |
951 | 953 |
|
952 | 954 |
/// \brief Return the flow on the given arc. |
953 | 955 |
/// |
954 | 956 |
/// This function returns the flow on the given arc. |
955 | 957 |
/// |
956 | 958 |
/// \pre \ref run() must be called before using this function. |
957 | 959 |
Value flow(const Arc& a) const { |
958 | 960 |
return _flow[_arc_id[a]]; |
959 | 961 |
} |
960 | 962 |
|
961 | 963 |
/// \brief Return the flow map (the primal solution). |
962 | 964 |
/// |
963 | 965 |
/// This function copies the flow value on each arc into the given |
964 | 966 |
/// map. The \c Value type of the algorithm must be convertible to |
965 | 967 |
/// the \c Value type of the map. |
966 | 968 |
/// |
967 | 969 |
/// \pre \ref run() must be called before using this function. |
968 | 970 |
template <typename FlowMap> |
969 | 971 |
void flowMap(FlowMap &map) const { |
970 | 972 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
971 | 973 |
map.set(a, _flow[_arc_id[a]]); |
972 | 974 |
} |
973 | 975 |
} |
974 | 976 |
|
975 | 977 |
/// \brief Return the potential (dual value) of the given node. |
976 | 978 |
/// |
977 | 979 |
/// This function returns the potential (dual value) of the |
978 | 980 |
/// given node. |
979 | 981 |
/// |
980 | 982 |
/// \pre \ref run() must be called before using this function. |
981 | 983 |
Cost potential(const Node& n) const { |
982 | 984 |
return _pi[_node_id[n]]; |
983 | 985 |
} |
984 | 986 |
|
985 | 987 |
/// \brief Return the potential map (the dual solution). |
986 | 988 |
/// |
987 | 989 |
/// This function copies the potential (dual value) of each node |
988 | 990 |
/// into the given map. |
989 | 991 |
/// The \c Cost type of the algorithm must be convertible to the |
990 | 992 |
/// \c Value type of the map. |
991 | 993 |
/// |
992 | 994 |
/// \pre \ref run() must be called before using this function. |
993 | 995 |
template <typename PotentialMap> |
994 | 996 |
void potentialMap(PotentialMap &map) const { |
995 | 997 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
996 | 998 |
map.set(n, _pi[_node_id[n]]); |
997 | 999 |
} |
998 | 1000 |
} |
999 | 1001 |
|
1000 | 1002 |
/// @} |
1001 | 1003 |
|
1002 | 1004 |
private: |
1003 | 1005 |
|
1004 | 1006 |
// Initialize internal data structures |
1005 | 1007 |
bool init() { |
1006 | 1008 |
if (_node_num == 0) return false; |
1007 | 1009 |
|
1008 | 1010 |
// Check the sum of supply values |
1009 | 1011 |
_sum_supply = 0; |
1010 | 1012 |
for (int i = 0; i != _node_num; ++i) { |
1011 | 1013 |
_sum_supply += _supply[i]; |
1012 | 1014 |
} |
1013 | 1015 |
if ( !((_stype == GEQ && _sum_supply <= 0) || |
1014 | 1016 |
(_stype == LEQ && _sum_supply >= 0)) ) return false; |
1015 | 1017 |
|
1016 | 1018 |
// Remove non-zero lower bounds |
1017 | 1019 |
if (_have_lower) { |
1018 | 1020 |
for (int i = 0; i != _arc_num; ++i) { |
1019 | 1021 |
Value c = _lower[i]; |
1020 | 1022 |
if (c >= 0) { |
1021 | 1023 |
_cap[i] = _upper[i] < INF ? _upper[i] - c : INF; |
1022 | 1024 |
} else { |
1023 | 1025 |
_cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF; |
1024 | 1026 |
} |
1025 | 1027 |
_supply[_source[i]] -= c; |
1026 | 1028 |
_supply[_target[i]] += c; |
1027 | 1029 |
} |
1028 | 1030 |
} else { |
1029 | 1031 |
for (int i = 0; i != _arc_num; ++i) { |
1030 | 1032 |
_cap[i] = _upper[i]; |
1031 | 1033 |
} |
1032 | 1034 |
} |
1033 | 1035 |
|
1034 | 1036 |
// Initialize artifical cost |
1035 | 1037 |
Cost ART_COST; |
1036 | 1038 |
if (std::numeric_limits<Cost>::is_exact) { |
1037 | 1039 |
ART_COST = std::numeric_limits<Cost>::max() / 2 + 1; |
1038 | 1040 |
} else { |
1039 | 1041 |
ART_COST = std::numeric_limits<Cost>::min(); |
1040 | 1042 |
for (int i = 0; i != _arc_num; ++i) { |
1041 | 1043 |
if (_cost[i] > ART_COST) ART_COST = _cost[i]; |
1042 | 1044 |
} |
1043 | 1045 |
ART_COST = (ART_COST + 1) * _node_num; |
1044 | 1046 |
} |
1045 | 1047 |
|
1046 | 1048 |
// Initialize arc maps |
1047 | 1049 |
for (int i = 0; i != _arc_num; ++i) { |
1048 | 1050 |
_flow[i] = 0; |
1049 | 1051 |
_state[i] = STATE_LOWER; |
1050 | 1052 |
} |
1051 | 1053 |
|
1052 | 1054 |
// Set data for the artificial root node |
1053 | 1055 |
_root = _node_num; |
1054 | 1056 |
_parent[_root] = -1; |
1055 | 1057 |
_pred[_root] = -1; |
1056 | 1058 |
_thread[_root] = 0; |
1057 | 1059 |
_rev_thread[0] = _root; |
1058 | 1060 |
_succ_num[_root] = _node_num + 1; |
1059 | 1061 |
_last_succ[_root] = _root - 1; |
1060 | 1062 |
_supply[_root] = -_sum_supply; |
1061 | 1063 |
_pi[_root] = 0; |
1062 | 1064 |
|
1063 | 1065 |
// Add artificial arcs and initialize the spanning tree data structure |
1064 | 1066 |
if (_sum_supply == 0) { |
1065 | 1067 |
// EQ supply constraints |
1066 | 1068 |
_search_arc_num = _arc_num; |
1067 | 1069 |
_all_arc_num = _arc_num + _node_num; |
1068 | 1070 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
1069 | 1071 |
_parent[u] = _root; |
1070 | 1072 |
_pred[u] = e; |
1071 | 1073 |
_thread[u] = u + 1; |
1072 | 1074 |
_rev_thread[u + 1] = u; |
1073 | 1075 |
_succ_num[u] = 1; |
1074 | 1076 |
_last_succ[u] = u; |
1075 | 1077 |
_cap[e] = INF; |
1076 | 1078 |
_state[e] = STATE_TREE; |
1077 | 1079 |
if (_supply[u] >= 0) { |
1078 | 1080 |
_forward[u] = true; |
1079 | 1081 |
_pi[u] = 0; |
1080 | 1082 |
_source[e] = u; |
1081 | 1083 |
_target[e] = _root; |
1082 | 1084 |
_flow[e] = _supply[u]; |
1083 | 1085 |
_cost[e] = 0; |
1084 | 1086 |
} else { |
1085 | 1087 |
_forward[u] = false; |
1086 | 1088 |
_pi[u] = ART_COST; |
1087 | 1089 |
_source[e] = _root; |
1088 | 1090 |
_target[e] = u; |
1089 | 1091 |
_flow[e] = -_supply[u]; |
1090 | 1092 |
_cost[e] = ART_COST; |
1091 | 1093 |
} |
1092 | 1094 |
} |
1093 | 1095 |
} |
1094 | 1096 |
else if (_sum_supply > 0) { |
1095 | 1097 |
// LEQ supply constraints |
1096 | 1098 |
_search_arc_num = _arc_num + _node_num; |
1097 | 1099 |
int f = _arc_num + _node_num; |
1098 | 1100 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
1099 | 1101 |
_parent[u] = _root; |
1100 | 1102 |
_thread[u] = u + 1; |
1101 | 1103 |
_rev_thread[u + 1] = u; |
1102 | 1104 |
_succ_num[u] = 1; |
1103 | 1105 |
_last_succ[u] = u; |
1104 | 1106 |
if (_supply[u] >= 0) { |
1105 | 1107 |
_forward[u] = true; |
1106 | 1108 |
_pi[u] = 0; |
1107 | 1109 |
_pred[u] = e; |
1108 | 1110 |
_source[e] = u; |
1109 | 1111 |
_target[e] = _root; |
1110 | 1112 |
_cap[e] = INF; |
1111 | 1113 |
_flow[e] = _supply[u]; |
1112 | 1114 |
_cost[e] = 0; |
1113 | 1115 |
_state[e] = STATE_TREE; |
1114 | 1116 |
} else { |
1115 | 1117 |
_forward[u] = false; |
1116 | 1118 |
_pi[u] = ART_COST; |
1117 | 1119 |
_pred[u] = f; |
1118 | 1120 |
_source[f] = _root; |
1119 | 1121 |
_target[f] = u; |
1120 | 1122 |
_cap[f] = INF; |
1121 | 1123 |
_flow[f] = -_supply[u]; |
1122 | 1124 |
_cost[f] = ART_COST; |
1123 | 1125 |
_state[f] = STATE_TREE; |
1124 | 1126 |
_source[e] = u; |
1125 | 1127 |
_target[e] = _root; |
1126 | 1128 |
_cap[e] = INF; |
1127 | 1129 |
_flow[e] = 0; |
1128 | 1130 |
_cost[e] = 0; |
1129 | 1131 |
_state[e] = STATE_LOWER; |
1130 | 1132 |
++f; |
1131 | 1133 |
} |
1132 | 1134 |
} |
1133 | 1135 |
_all_arc_num = f; |
1134 | 1136 |
} |
1135 | 1137 |
else { |
1136 | 1138 |
// GEQ supply constraints |
1137 | 1139 |
_search_arc_num = _arc_num + _node_num; |
1138 | 1140 |
int f = _arc_num + _node_num; |
1139 | 1141 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
1140 | 1142 |
_parent[u] = _root; |
1141 | 1143 |
_thread[u] = u + 1; |
1142 | 1144 |
_rev_thread[u + 1] = u; |
1143 | 1145 |
_succ_num[u] = 1; |
1144 | 1146 |
_last_succ[u] = u; |
1145 | 1147 |
if (_supply[u] <= 0) { |
1146 | 1148 |
_forward[u] = false; |
1147 | 1149 |
_pi[u] = 0; |
1148 | 1150 |
_pred[u] = e; |
1149 | 1151 |
_source[e] = _root; |
1150 | 1152 |
_target[e] = u; |
1151 | 1153 |
_cap[e] = INF; |
1152 | 1154 |
_flow[e] = -_supply[u]; |
1153 | 1155 |
_cost[e] = 0; |
1154 | 1156 |
_state[e] = STATE_TREE; |
1155 | 1157 |
} else { |
1156 | 1158 |
_forward[u] = true; |
1157 | 1159 |
_pi[u] = -ART_COST; |
1158 | 1160 |
_pred[u] = f; |
1159 | 1161 |
_source[f] = u; |
1160 | 1162 |
_target[f] = _root; |
1161 | 1163 |
_cap[f] = INF; |
1162 | 1164 |
_flow[f] = _supply[u]; |
1163 | 1165 |
_state[f] = STATE_TREE; |
1164 | 1166 |
_cost[f] = ART_COST; |
1165 | 1167 |
_source[e] = _root; |
1166 | 1168 |
_target[e] = u; |
1167 | 1169 |
_cap[e] = INF; |
1168 | 1170 |
_flow[e] = 0; |
1169 | 1171 |
_cost[e] = 0; |
1170 | 1172 |
_state[e] = STATE_LOWER; |
1171 | 1173 |
++f; |
1172 | 1174 |
} |
1173 | 1175 |
} |
1174 | 1176 |
_all_arc_num = f; |
1175 | 1177 |
} |
1176 | 1178 |
|
1177 | 1179 |
return true; |
1178 | 1180 |
} |
1179 | 1181 |
|
1180 | 1182 |
// Find the join node |
1181 | 1183 |
void findJoinNode() { |
1182 | 1184 |
int u = _source[in_arc]; |
1183 | 1185 |
int v = _target[in_arc]; |
1184 | 1186 |
while (u != v) { |
1185 | 1187 |
if (_succ_num[u] < _succ_num[v]) { |
1186 | 1188 |
u = _parent[u]; |
1187 | 1189 |
} else { |
1188 | 1190 |
v = _parent[v]; |
1189 | 1191 |
} |
1190 | 1192 |
} |
1191 | 1193 |
join = u; |
1192 | 1194 |
} |
1193 | 1195 |
|
1194 | 1196 |
// Find the leaving arc of the cycle and returns true if the |
1195 | 1197 |
// leaving arc is not the same as the entering arc |
1196 | 1198 |
bool findLeavingArc() { |
1197 | 1199 |
// Initialize first and second nodes according to the direction |
1198 | 1200 |
// of the cycle |
1199 | 1201 |
if (_state[in_arc] == STATE_LOWER) { |
1200 | 1202 |
first = _source[in_arc]; |
1201 | 1203 |
second = _target[in_arc]; |
1202 | 1204 |
} else { |
1203 | 1205 |
first = _target[in_arc]; |
1204 | 1206 |
second = _source[in_arc]; |
1205 | 1207 |
} |
1206 | 1208 |
delta = _cap[in_arc]; |
1207 | 1209 |
int result = 0; |
1208 | 1210 |
Value d; |
1209 | 1211 |
int e; |
1210 | 1212 |
|
1211 | 1213 |
// Search the cycle along the path form the first node to the root |
1212 | 1214 |
for (int u = first; u != join; u = _parent[u]) { |
1213 | 1215 |
e = _pred[u]; |
1214 | 1216 |
d = _forward[u] ? |
1215 | 1217 |
_flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]); |
1216 | 1218 |
if (d < delta) { |
1217 | 1219 |
delta = d; |
1218 | 1220 |
u_out = u; |
1219 | 1221 |
result = 1; |
1220 | 1222 |
} |
1221 | 1223 |
} |
1222 | 1224 |
// Search the cycle along the path form the second node to the root |
1223 | 1225 |
for (int u = second; u != join; u = _parent[u]) { |
1224 | 1226 |
e = _pred[u]; |
1225 | 1227 |
d = _forward[u] ? |
1226 | 1228 |
(_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e]; |
1227 | 1229 |
if (d <= delta) { |
1228 | 1230 |
delta = d; |
1229 | 1231 |
u_out = u; |
1230 | 1232 |
result = 2; |
1231 | 1233 |
} |
1232 | 1234 |
} |
1233 | 1235 |
|
1234 | 1236 |
if (result == 1) { |
1235 | 1237 |
u_in = first; |
1236 | 1238 |
v_in = second; |
1237 | 1239 |
} else { |
1238 | 1240 |
u_in = second; |
1239 | 1241 |
v_in = first; |
1240 | 1242 |
} |
1241 | 1243 |
return result != 0; |
1242 | 1244 |
} |
1243 | 1245 |
|
1244 | 1246 |
// Change _flow and _state vectors |
1245 | 1247 |
void changeFlow(bool change) { |
1246 | 1248 |
// Augment along the cycle |
1247 | 1249 |
if (delta > 0) { |
1248 | 1250 |
Value val = _state[in_arc] * delta; |
1249 | 1251 |
_flow[in_arc] += val; |
1250 | 1252 |
for (int u = _source[in_arc]; u != join; u = _parent[u]) { |
1251 | 1253 |
_flow[_pred[u]] += _forward[u] ? -val : val; |
1252 | 1254 |
} |
1253 | 1255 |
for (int u = _target[in_arc]; u != join; u = _parent[u]) { |
1254 | 1256 |
_flow[_pred[u]] += _forward[u] ? val : -val; |
1255 | 1257 |
} |
1256 | 1258 |
} |
1257 | 1259 |
// Update the state of the entering and leaving arcs |
1258 | 1260 |
if (change) { |
1259 | 1261 |
_state[in_arc] = STATE_TREE; |
1260 | 1262 |
_state[_pred[u_out]] = |
1261 | 1263 |
(_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; |
1262 | 1264 |
} else { |
1263 | 1265 |
_state[in_arc] = -_state[in_arc]; |
1264 | 1266 |
} |
1265 | 1267 |
} |
1266 | 1268 |
|
1267 | 1269 |
// Update the tree structure |
1268 | 1270 |
void updateTreeStructure() { |
1269 | 1271 |
int u, w; |
1270 | 1272 |
int old_rev_thread = _rev_thread[u_out]; |
1271 | 1273 |
int old_succ_num = _succ_num[u_out]; |
1272 | 1274 |
int old_last_succ = _last_succ[u_out]; |
1273 | 1275 |
v_out = _parent[u_out]; |
1274 | 1276 |
|
1275 | 1277 |
u = _last_succ[u_in]; // the last successor of u_in |
1276 | 1278 |
right = _thread[u]; // the node after it |
1277 | 1279 |
|
1278 | 1280 |
// Handle the case when old_rev_thread equals to v_in |
1279 | 1281 |
// (it also means that join and v_out coincide) |
1280 | 1282 |
if (old_rev_thread == v_in) { |
1281 | 1283 |
last = _thread[_last_succ[u_out]]; |
1282 | 1284 |
} else { |
1283 | 1285 |
last = _thread[v_in]; |
1284 | 1286 |
} |
1285 | 1287 |
|
1286 | 1288 |
// Update _thread and _parent along the stem nodes (i.e. the nodes |
1287 | 1289 |
// between u_in and u_out, whose parent have to be changed) |
1288 | 1290 |
_thread[v_in] = stem = u_in; |
1289 | 1291 |
_dirty_revs.clear(); |
1290 | 1292 |
_dirty_revs.push_back(v_in); |
1291 | 1293 |
par_stem = v_in; |
1292 | 1294 |
while (stem != u_out) { |
1293 | 1295 |
// Insert the next stem node into the thread list |
1294 | 1296 |
new_stem = _parent[stem]; |
1295 | 1297 |
_thread[u] = new_stem; |
1296 | 1298 |
_dirty_revs.push_back(u); |
1297 | 1299 |
|
1298 | 1300 |
// Remove the subtree of stem from the thread list |
1299 | 1301 |
w = _rev_thread[stem]; |
1300 | 1302 |
_thread[w] = right; |
1301 | 1303 |
_rev_thread[right] = w; |
1302 | 1304 |
|
1303 | 1305 |
// Change the parent node and shift stem nodes |
1304 | 1306 |
_parent[stem] = par_stem; |
1305 | 1307 |
par_stem = stem; |
1306 | 1308 |
stem = new_stem; |
1307 | 1309 |
|
1308 | 1310 |
// Update u and right |
1309 | 1311 |
u = _last_succ[stem] == _last_succ[par_stem] ? |
1310 | 1312 |
_rev_thread[par_stem] : _last_succ[stem]; |
1311 | 1313 |
right = _thread[u]; |
1312 | 1314 |
} |
1313 | 1315 |
_parent[u_out] = par_stem; |
1314 | 1316 |
_thread[u] = last; |
1315 | 1317 |
_rev_thread[last] = u; |
1316 | 1318 |
_last_succ[u_out] = u; |
1317 | 1319 |
|
1318 | 1320 |
// Remove the subtree of u_out from the thread list except for |
1319 | 1321 |
// the case when old_rev_thread equals to v_in |
1320 | 1322 |
// (it also means that join and v_out coincide) |
1321 | 1323 |
if (old_rev_thread != v_in) { |
1322 | 1324 |
_thread[old_rev_thread] = right; |
1323 | 1325 |
_rev_thread[right] = old_rev_thread; |
1324 | 1326 |
} |
1325 | 1327 |
|
1326 | 1328 |
// Update _rev_thread using the new _thread values |
1327 | 1329 |
for (int i = 0; i < int(_dirty_revs.size()); ++i) { |
1328 | 1330 |
u = _dirty_revs[i]; |
1329 | 1331 |
_rev_thread[_thread[u]] = u; |
1330 | 1332 |
} |
1331 | 1333 |
|
1332 | 1334 |
// Update _pred, _forward, _last_succ and _succ_num for the |
1333 | 1335 |
// stem nodes from u_out to u_in |
1334 | 1336 |
int tmp_sc = 0, tmp_ls = _last_succ[u_out]; |
1335 | 1337 |
u = u_out; |
1336 | 1338 |
while (u != u_in) { |
1337 | 1339 |
w = _parent[u]; |
1338 | 1340 |
_pred[u] = _pred[w]; |
1339 | 1341 |
_forward[u] = !_forward[w]; |
1340 | 1342 |
tmp_sc += _succ_num[u] - _succ_num[w]; |
1341 | 1343 |
_succ_num[u] = tmp_sc; |
1342 | 1344 |
_last_succ[w] = tmp_ls; |
1343 | 1345 |
u = w; |
1344 | 1346 |
} |
1345 | 1347 |
_pred[u_in] = in_arc; |
1346 | 1348 |
_forward[u_in] = (u_in == _source[in_arc]); |
1347 | 1349 |
_succ_num[u_in] = old_succ_num; |
1348 | 1350 |
|
1349 | 1351 |
// Set limits for updating _last_succ form v_in and v_out |
1350 | 1352 |
// towards the root |
1351 | 1353 |
int up_limit_in = -1; |
1352 | 1354 |
int up_limit_out = -1; |
1353 | 1355 |
if (_last_succ[join] == v_in) { |
1354 | 1356 |
up_limit_out = join; |
1355 | 1357 |
} else { |
1356 | 1358 |
up_limit_in = join; |
1357 | 1359 |
} |
1358 | 1360 |
|
1359 | 1361 |
// Update _last_succ from v_in towards the root |
1360 | 1362 |
for (u = v_in; u != up_limit_in && _last_succ[u] == v_in; |
1361 | 1363 |
u = _parent[u]) { |
1362 | 1364 |
_last_succ[u] = _last_succ[u_out]; |
1363 | 1365 |
} |
1364 | 1366 |
// Update _last_succ from v_out towards the root |
1365 | 1367 |
if (join != old_rev_thread && v_in != old_rev_thread) { |
1366 | 1368 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; |
1367 | 1369 |
u = _parent[u]) { |
1368 | 1370 |
_last_succ[u] = old_rev_thread; |
1369 | 1371 |
} |
1370 | 1372 |
} else { |
1371 | 1373 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ; |
1372 | 1374 |
u = _parent[u]) { |
1373 | 1375 |
_last_succ[u] = _last_succ[u_out]; |
1374 | 1376 |
} |
1375 | 1377 |
} |
1376 | 1378 |
|
1377 | 1379 |
// Update _succ_num from v_in to join |
1378 | 1380 |
for (u = v_in; u != join; u = _parent[u]) { |
1379 | 1381 |
_succ_num[u] += old_succ_num; |
1380 | 1382 |
} |
1381 | 1383 |
// Update _succ_num from v_out to join |
1382 | 1384 |
for (u = v_out; u != join; u = _parent[u]) { |
1383 | 1385 |
_succ_num[u] -= old_succ_num; |
1384 | 1386 |
} |
1385 | 1387 |
} |
1386 | 1388 |
|
1387 | 1389 |
// Update potentials |
1388 | 1390 |
void updatePotential() { |
1389 | 1391 |
Cost sigma = _forward[u_in] ? |
1390 | 1392 |
_pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] : |
1391 | 1393 |
_pi[v_in] - _pi[u_in] + _cost[_pred[u_in]]; |
1392 | 1394 |
// Update potentials in the subtree, which has been moved |
1393 | 1395 |
int end = _thread[_last_succ[u_in]]; |
1394 | 1396 |
for (int u = u_in; u != end; u = _thread[u]) { |
1395 | 1397 |
_pi[u] += sigma; |
1396 | 1398 |
} |
1397 | 1399 |
} |
1398 | 1400 |
|
1399 | 1401 |
// Execute the algorithm |
1400 | 1402 |
ProblemType start(PivotRule pivot_rule) { |
1401 | 1403 |
// Select the pivot rule implementation |
1402 | 1404 |
switch (pivot_rule) { |
1403 | 1405 |
case FIRST_ELIGIBLE: |
1404 | 1406 |
return start<FirstEligiblePivotRule>(); |
1405 | 1407 |
case BEST_ELIGIBLE: |
1406 | 1408 |
return start<BestEligiblePivotRule>(); |
1407 | 1409 |
case BLOCK_SEARCH: |
1408 | 1410 |
return start<BlockSearchPivotRule>(); |
1409 | 1411 |
case CANDIDATE_LIST: |
1410 | 1412 |
return start<CandidateListPivotRule>(); |
1411 | 1413 |
case ALTERING_LIST: |
1412 | 1414 |
return start<AlteringListPivotRule>(); |
1413 | 1415 |
} |
1414 | 1416 |
return INFEASIBLE; // avoid warning |
1415 | 1417 |
} |
1416 | 1418 |
|
1417 | 1419 |
template <typename PivotRuleImpl> |
1418 | 1420 |
ProblemType start() { |
1419 | 1421 |
PivotRuleImpl pivot(*this); |
1420 | 1422 |
|
1421 | 1423 |
// Execute the Network Simplex algorithm |
1422 | 1424 |
while (pivot.findEnteringArc()) { |
1423 | 1425 |
findJoinNode(); |
1424 | 1426 |
bool change = findLeavingArc(); |
1425 | 1427 |
if (delta >= INF) return UNBOUNDED; |
1426 | 1428 |
changeFlow(change); |
1427 | 1429 |
if (change) { |
1428 | 1430 |
updateTreeStructure(); |
1429 | 1431 |
updatePotential(); |
1430 | 1432 |
} |
1431 | 1433 |
} |
1432 | 1434 |
|
1433 | 1435 |
// Check feasibility |
1434 | 1436 |
for (int e = _search_arc_num; e != _all_arc_num; ++e) { |
1435 | 1437 |
if (_flow[e] != 0) return INFEASIBLE; |
1436 | 1438 |
} |
1437 | 1439 |
|
1438 | 1440 |
// Transform the solution and the supply map to the original form |
1439 | 1441 |
if (_have_lower) { |
1440 | 1442 |
for (int i = 0; i != _arc_num; ++i) { |
1441 | 1443 |
Value c = _lower[i]; |
1442 | 1444 |
if (c != 0) { |
1443 | 1445 |
_flow[i] += c; |
1444 | 1446 |
_supply[_source[i]] += c; |
1445 | 1447 |
_supply[_target[i]] -= c; |
1446 | 1448 |
} |
1447 | 1449 |
} |
1448 | 1450 |
} |
1449 | 1451 |
|
1450 | 1452 |
// Shift potentials to meet the requirements of the GEQ/LEQ type |
1451 | 1453 |
// optimality conditions |
1452 | 1454 |
if (_sum_supply == 0) { |
1453 | 1455 |
if (_stype == GEQ) { |
1454 | 1456 |
Cost max_pot = std::numeric_limits<Cost>::min(); |
1455 | 1457 |
for (int i = 0; i != _node_num; ++i) { |
1456 | 1458 |
if (_pi[i] > max_pot) max_pot = _pi[i]; |
1457 | 1459 |
} |
1458 | 1460 |
if (max_pot > 0) { |
1459 | 1461 |
for (int i = 0; i != _node_num; ++i) |
1460 | 1462 |
_pi[i] -= max_pot; |
1461 | 1463 |
} |
1462 | 1464 |
} else { |
1463 | 1465 |
Cost min_pot = std::numeric_limits<Cost>::max(); |
1464 | 1466 |
for (int i = 0; i != _node_num; ++i) { |
1465 | 1467 |
if (_pi[i] < min_pot) min_pot = _pi[i]; |
1466 | 1468 |
} |
1467 | 1469 |
if (min_pot < 0) { |
1468 | 1470 |
for (int i = 0; i != _node_num; ++i) |
1469 | 1471 |
_pi[i] -= min_pot; |
1470 | 1472 |
} |
1471 | 1473 |
} |
1472 | 1474 |
} |
1473 | 1475 |
|
1474 | 1476 |
return OPTIMAL; |
1475 | 1477 |
} |
1476 | 1478 |
|
1477 | 1479 |
}; //class NetworkSimplex |
1478 | 1480 |
|
1479 | 1481 |
///@} |
1480 | 1482 |
|
1481 | 1483 |
} //namespace lemon |
1482 | 1484 |
|
1483 | 1485 |
#endif //LEMON_NETWORK_SIMPLEX_H |
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_PREFLOW_H |
20 | 20 |
#define LEMON_PREFLOW_H |
21 | 21 |
|
22 | 22 |
#include <lemon/tolerance.h> |
23 | 23 |
#include <lemon/elevator.h> |
24 | 24 |
|
25 | 25 |
/// \file |
26 | 26 |
/// \ingroup max_flow |
27 | 27 |
/// \brief Implementation of the preflow algorithm. |
28 | 28 |
|
29 | 29 |
namespace lemon { |
30 | 30 |
|
31 | 31 |
/// \brief Default traits class of Preflow class. |
32 | 32 |
/// |
33 | 33 |
/// Default traits class of Preflow class. |
34 | 34 |
/// \tparam GR Digraph type. |
35 | 35 |
/// \tparam CAP Capacity map type. |
36 | 36 |
template <typename GR, typename CAP> |
37 | 37 |
struct PreflowDefaultTraits { |
38 | 38 |
|
39 | 39 |
/// \brief The type of the digraph the algorithm runs on. |
40 | 40 |
typedef GR Digraph; |
41 | 41 |
|
42 | 42 |
/// \brief The type of the map that stores the arc capacities. |
43 | 43 |
/// |
44 | 44 |
/// The type of the map that stores the arc capacities. |
45 | 45 |
/// It must meet the \ref concepts::ReadMap "ReadMap" concept. |
46 | 46 |
typedef CAP CapacityMap; |
47 | 47 |
|
48 | 48 |
/// \brief The type of the flow values. |
49 | 49 |
typedef typename CapacityMap::Value Value; |
50 | 50 |
|
51 | 51 |
/// \brief The type of the map that stores the flow values. |
52 | 52 |
/// |
53 | 53 |
/// The type of the map that stores the flow values. |
54 | 54 |
/// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. |
55 | 55 |
#ifdef DOXYGEN |
56 | 56 |
typedef GR::ArcMap<Value> FlowMap; |
57 | 57 |
#else |
58 | 58 |
typedef typename Digraph::template ArcMap<Value> FlowMap; |
59 | 59 |
#endif |
60 | 60 |
|
61 | 61 |
/// \brief Instantiates a FlowMap. |
62 | 62 |
/// |
63 | 63 |
/// This function instantiates a \ref FlowMap. |
64 | 64 |
/// \param digraph The digraph for which we would like to define |
65 | 65 |
/// the flow map. |
66 | 66 |
static FlowMap* createFlowMap(const Digraph& digraph) { |
67 | 67 |
return new FlowMap(digraph); |
68 | 68 |
} |
69 | 69 |
|
70 | 70 |
/// \brief The elevator type used by Preflow algorithm. |
71 | 71 |
/// |
72 | 72 |
/// The elevator type used by Preflow algorithm. |
73 | 73 |
/// |
74 | 74 |
/// \sa Elevator, LinkedElevator |
75 | 75 |
#ifdef DOXYGEN |
76 | 76 |
typedef lemon::Elevator<GR, GR::Node> Elevator; |
77 | 77 |
#else |
78 | 78 |
typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator; |
79 | 79 |
#endif |
80 | 80 |
|
81 | 81 |
/// \brief Instantiates an Elevator. |
82 | 82 |
/// |
83 | 83 |
/// This function instantiates an \ref Elevator. |
84 | 84 |
/// \param digraph The digraph for which we would like to define |
85 | 85 |
/// the elevator. |
86 | 86 |
/// \param max_level The maximum level of the elevator. |
87 | 87 |
static Elevator* createElevator(const Digraph& digraph, int max_level) { |
88 | 88 |
return new Elevator(digraph, max_level); |
89 | 89 |
} |
90 | 90 |
|
91 | 91 |
/// \brief The tolerance used by the algorithm |
92 | 92 |
/// |
93 | 93 |
/// The tolerance used by the algorithm to handle inexact computation. |
94 | 94 |
typedef lemon::Tolerance<Value> Tolerance; |
95 | 95 |
|
96 | 96 |
}; |
97 | 97 |
|
98 | 98 |
|
99 | 99 |
/// \ingroup max_flow |
100 | 100 |
/// |
101 | 101 |
/// \brief %Preflow algorithm class. |
102 | 102 |
/// |
103 | 103 |
/// This class provides an implementation of Goldberg-Tarjan's \e preflow |
104 | 104 |
/// \e push-relabel algorithm producing a \ref max_flow |
105 |
/// "flow of maximum value" in a digraph |
|
105 |
/// "flow of maximum value" in a digraph \ref clrs01algorithms, |
|
106 |
/// \ref amo93networkflows, \ref goldberg88newapproach. |
|
106 | 107 |
/// The preflow algorithms are the fastest known maximum |
107 | 108 |
/// flow algorithms. The current implementation uses a mixture of the |
108 | 109 |
/// \e "highest label" and the \e "bound decrease" heuristics. |
109 | 110 |
/// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$. |
110 | 111 |
/// |
111 | 112 |
/// The algorithm consists of two phases. After the first phase |
112 | 113 |
/// the maximum flow value and the minimum cut is obtained. The |
113 | 114 |
/// second phase constructs a feasible maximum flow on each arc. |
114 | 115 |
/// |
115 | 116 |
/// \tparam GR The type of the digraph the algorithm runs on. |
116 | 117 |
/// \tparam CAP The type of the capacity map. The default map |
117 | 118 |
/// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>". |
118 | 119 |
#ifdef DOXYGEN |
119 | 120 |
template <typename GR, typename CAP, typename TR> |
120 | 121 |
#else |
121 | 122 |
template <typename GR, |
122 | 123 |
typename CAP = typename GR::template ArcMap<int>, |
123 | 124 |
typename TR = PreflowDefaultTraits<GR, CAP> > |
124 | 125 |
#endif |
125 | 126 |
class Preflow { |
126 | 127 |
public: |
127 | 128 |
|
128 | 129 |
///The \ref PreflowDefaultTraits "traits class" of the algorithm. |
129 | 130 |
typedef TR Traits; |
130 | 131 |
///The type of the digraph the algorithm runs on. |
131 | 132 |
typedef typename Traits::Digraph Digraph; |
132 | 133 |
///The type of the capacity map. |
133 | 134 |
typedef typename Traits::CapacityMap CapacityMap; |
134 | 135 |
///The type of the flow values. |
135 | 136 |
typedef typename Traits::Value Value; |
136 | 137 |
|
137 | 138 |
///The type of the flow map. |
138 | 139 |
typedef typename Traits::FlowMap FlowMap; |
139 | 140 |
///The type of the elevator. |
140 | 141 |
typedef typename Traits::Elevator Elevator; |
141 | 142 |
///The type of the tolerance. |
142 | 143 |
typedef typename Traits::Tolerance Tolerance; |
143 | 144 |
|
144 | 145 |
private: |
145 | 146 |
|
146 | 147 |
TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
147 | 148 |
|
148 | 149 |
const Digraph& _graph; |
149 | 150 |
const CapacityMap* _capacity; |
150 | 151 |
|
151 | 152 |
int _node_num; |
152 | 153 |
|
153 | 154 |
Node _source, _target; |
154 | 155 |
|
155 | 156 |
FlowMap* _flow; |
156 | 157 |
bool _local_flow; |
157 | 158 |
|
158 | 159 |
Elevator* _level; |
159 | 160 |
bool _local_level; |
160 | 161 |
|
161 | 162 |
typedef typename Digraph::template NodeMap<Value> ExcessMap; |
162 | 163 |
ExcessMap* _excess; |
163 | 164 |
|
164 | 165 |
Tolerance _tolerance; |
165 | 166 |
|
166 | 167 |
bool _phase; |
167 | 168 |
|
168 | 169 |
|
169 | 170 |
void createStructures() { |
170 | 171 |
_node_num = countNodes(_graph); |
171 | 172 |
|
172 | 173 |
if (!_flow) { |
173 | 174 |
_flow = Traits::createFlowMap(_graph); |
174 | 175 |
_local_flow = true; |
175 | 176 |
} |
176 | 177 |
if (!_level) { |
177 | 178 |
_level = Traits::createElevator(_graph, _node_num); |
178 | 179 |
_local_level = true; |
179 | 180 |
} |
180 | 181 |
if (!_excess) { |
181 | 182 |
_excess = new ExcessMap(_graph); |
182 | 183 |
} |
183 | 184 |
} |
184 | 185 |
|
185 | 186 |
void destroyStructures() { |
186 | 187 |
if (_local_flow) { |
187 | 188 |
delete _flow; |
188 | 189 |
} |
189 | 190 |
if (_local_level) { |
190 | 191 |
delete _level; |
191 | 192 |
} |
192 | 193 |
if (_excess) { |
193 | 194 |
delete _excess; |
194 | 195 |
} |
195 | 196 |
} |
196 | 197 |
|
197 | 198 |
public: |
198 | 199 |
|
199 | 200 |
typedef Preflow Create; |
200 | 201 |
|
201 | 202 |
///\name Named Template Parameters |
202 | 203 |
|
203 | 204 |
///@{ |
204 | 205 |
|
205 | 206 |
template <typename T> |
206 | 207 |
struct SetFlowMapTraits : public Traits { |
207 | 208 |
typedef T FlowMap; |
208 | 209 |
static FlowMap *createFlowMap(const Digraph&) { |
209 | 210 |
LEMON_ASSERT(false, "FlowMap is not initialized"); |
210 | 211 |
return 0; // ignore warnings |
211 | 212 |
} |
212 | 213 |
}; |
213 | 214 |
|
214 | 215 |
/// \brief \ref named-templ-param "Named parameter" for setting |
215 | 216 |
/// FlowMap type |
216 | 217 |
/// |
217 | 218 |
/// \ref named-templ-param "Named parameter" for setting FlowMap |
218 | 219 |
/// type. |
219 | 220 |
template <typename T> |
220 | 221 |
struct SetFlowMap |
221 | 222 |
: public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > { |
222 | 223 |
typedef Preflow<Digraph, CapacityMap, |
223 | 224 |
SetFlowMapTraits<T> > Create; |
224 | 225 |
}; |
225 | 226 |
|
226 | 227 |
template <typename T> |
227 | 228 |
struct SetElevatorTraits : public Traits { |
228 | 229 |
typedef T Elevator; |
229 | 230 |
static Elevator *createElevator(const Digraph&, int) { |
230 | 231 |
LEMON_ASSERT(false, "Elevator is not initialized"); |
231 | 232 |
return 0; // ignore warnings |
232 | 233 |
} |
233 | 234 |
}; |
234 | 235 |
|
235 | 236 |
/// \brief \ref named-templ-param "Named parameter" for setting |
236 | 237 |
/// Elevator type |
237 | 238 |
/// |
238 | 239 |
/// \ref named-templ-param "Named parameter" for setting Elevator |
239 | 240 |
/// type. If this named parameter is used, then an external |
240 | 241 |
/// elevator object must be passed to the algorithm using the |
241 | 242 |
/// \ref elevator(Elevator&) "elevator()" function before calling |
242 | 243 |
/// \ref run() or \ref init(). |
243 | 244 |
/// \sa SetStandardElevator |
244 | 245 |
template <typename T> |
245 | 246 |
struct SetElevator |
246 | 247 |
: public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > { |
247 | 248 |
typedef Preflow<Digraph, CapacityMap, |
248 | 249 |
SetElevatorTraits<T> > Create; |
249 | 250 |
}; |
250 | 251 |
|
251 | 252 |
template <typename T> |
252 | 253 |
struct SetStandardElevatorTraits : public Traits { |
253 | 254 |
typedef T Elevator; |
254 | 255 |
static Elevator *createElevator(const Digraph& digraph, int max_level) { |
255 | 256 |
return new Elevator(digraph, max_level); |
256 | 257 |
} |
257 | 258 |
}; |
258 | 259 |
|
259 | 260 |
/// \brief \ref named-templ-param "Named parameter" for setting |
260 | 261 |
/// Elevator type with automatic allocation |
261 | 262 |
/// |
262 | 263 |
/// \ref named-templ-param "Named parameter" for setting Elevator |
263 | 264 |
/// type with automatic allocation. |
264 | 265 |
/// The Elevator should have standard constructor interface to be |
265 | 266 |
/// able to automatically created by the algorithm (i.e. the |
266 | 267 |
/// digraph and the maximum level should be passed to it). |
267 | 268 |
/// However an external elevator object could also be passed to the |
268 | 269 |
/// algorithm with the \ref elevator(Elevator&) "elevator()" function |
269 | 270 |
/// before calling \ref run() or \ref init(). |
270 | 271 |
/// \sa SetElevator |
271 | 272 |
template <typename T> |
272 | 273 |
struct SetStandardElevator |
273 | 274 |
: public Preflow<Digraph, CapacityMap, |
274 | 275 |
SetStandardElevatorTraits<T> > { |
275 | 276 |
typedef Preflow<Digraph, CapacityMap, |
276 | 277 |
SetStandardElevatorTraits<T> > Create; |
277 | 278 |
}; |
278 | 279 |
|
279 | 280 |
/// @} |
280 | 281 |
|
281 | 282 |
protected: |
282 | 283 |
|
283 | 284 |
Preflow() {} |
284 | 285 |
|
285 | 286 |
public: |
286 | 287 |
|
287 | 288 |
|
288 | 289 |
/// \brief The constructor of the class. |
289 | 290 |
/// |
290 | 291 |
/// The constructor of the class. |
291 | 292 |
/// \param digraph The digraph the algorithm runs on. |
292 | 293 |
/// \param capacity The capacity of the arcs. |
293 | 294 |
/// \param source The source node. |
294 | 295 |
/// \param target The target node. |
295 | 296 |
Preflow(const Digraph& digraph, const CapacityMap& capacity, |
296 | 297 |
Node source, Node target) |
297 | 298 |
: _graph(digraph), _capacity(&capacity), |
298 | 299 |
_node_num(0), _source(source), _target(target), |
299 | 300 |
_flow(0), _local_flow(false), |
300 | 301 |
_level(0), _local_level(false), |
301 | 302 |
_excess(0), _tolerance(), _phase() {} |
302 | 303 |
|
303 | 304 |
/// \brief Destructor. |
304 | 305 |
/// |
305 | 306 |
/// Destructor. |
306 | 307 |
~Preflow() { |
307 | 308 |
destroyStructures(); |
308 | 309 |
} |
309 | 310 |
|
310 | 311 |
/// \brief Sets the capacity map. |
311 | 312 |
/// |
312 | 313 |
/// Sets the capacity map. |
313 | 314 |
/// \return <tt>(*this)</tt> |
314 | 315 |
Preflow& capacityMap(const CapacityMap& map) { |
315 | 316 |
_capacity = ↦ |
316 | 317 |
return *this; |
317 | 318 |
} |
318 | 319 |
|
319 | 320 |
/// \brief Sets the flow map. |
320 | 321 |
/// |
321 | 322 |
/// Sets the flow map. |
322 | 323 |
/// If you don't use this function before calling \ref run() or |
323 | 324 |
/// \ref init(), an instance will be allocated automatically. |
324 | 325 |
/// The destructor deallocates this automatically allocated map, |
325 | 326 |
/// of course. |
326 | 327 |
/// \return <tt>(*this)</tt> |
327 | 328 |
Preflow& flowMap(FlowMap& map) { |
328 | 329 |
if (_local_flow) { |
329 | 330 |
delete _flow; |
330 | 331 |
_local_flow = false; |
331 | 332 |
} |
332 | 333 |
_flow = ↦ |
333 | 334 |
return *this; |
334 | 335 |
} |
335 | 336 |
|
336 | 337 |
/// \brief Sets the source node. |
337 | 338 |
/// |
338 | 339 |
/// Sets the source node. |
339 | 340 |
/// \return <tt>(*this)</tt> |
340 | 341 |
Preflow& source(const Node& node) { |
341 | 342 |
_source = node; |
342 | 343 |
return *this; |
343 | 344 |
} |
344 | 345 |
|
345 | 346 |
/// \brief Sets the target node. |
346 | 347 |
/// |
347 | 348 |
/// Sets the target node. |
348 | 349 |
/// \return <tt>(*this)</tt> |
349 | 350 |
Preflow& target(const Node& node) { |
350 | 351 |
_target = node; |
351 | 352 |
return *this; |
352 | 353 |
} |
353 | 354 |
|
354 | 355 |
/// \brief Sets the elevator used by algorithm. |
355 | 356 |
/// |
356 | 357 |
/// Sets the elevator used by algorithm. |
357 | 358 |
/// If you don't use this function before calling \ref run() or |
358 | 359 |
/// \ref init(), an instance will be allocated automatically. |
359 | 360 |
/// The destructor deallocates this automatically allocated elevator, |
360 | 361 |
/// of course. |
361 | 362 |
/// \return <tt>(*this)</tt> |
362 | 363 |
Preflow& elevator(Elevator& elevator) { |
363 | 364 |
if (_local_level) { |
364 | 365 |
delete _level; |
365 | 366 |
_local_level = false; |
366 | 367 |
} |
367 | 368 |
_level = &elevator; |
368 | 369 |
return *this; |
369 | 370 |
} |
370 | 371 |
|
371 | 372 |
/// \brief Returns a const reference to the elevator. |
372 | 373 |
/// |
373 | 374 |
/// Returns a const reference to the elevator. |
374 | 375 |
/// |
375 | 376 |
/// \pre Either \ref run() or \ref init() must be called before |
376 | 377 |
/// using this function. |
377 | 378 |
const Elevator& elevator() const { |
378 | 379 |
return *_level; |
379 | 380 |
} |
380 | 381 |
|
381 | 382 |
/// \brief Sets the tolerance used by the algorithm. |
382 | 383 |
/// |
383 | 384 |
/// Sets the tolerance object used by the algorithm. |
384 | 385 |
/// \return <tt>(*this)</tt> |
385 | 386 |
Preflow& tolerance(const Tolerance& tolerance) { |
386 | 387 |
_tolerance = tolerance; |
387 | 388 |
return *this; |
388 | 389 |
} |
389 | 390 |
|
390 | 391 |
/// \brief Returns a const reference to the tolerance. |
391 | 392 |
/// |
392 | 393 |
/// Returns a const reference to the tolerance object used by |
393 | 394 |
/// the algorithm. |
394 | 395 |
const Tolerance& tolerance() const { |
395 | 396 |
return _tolerance; |
396 | 397 |
} |
397 | 398 |
|
398 | 399 |
/// \name Execution Control |
399 | 400 |
/// The simplest way to execute the preflow algorithm is to use |
400 | 401 |
/// \ref run() or \ref runMinCut().\n |
401 | 402 |
/// If you need better control on the initial solution or the execution, |
402 | 403 |
/// you have to call one of the \ref init() functions first, then |
403 | 404 |
/// \ref startFirstPhase() and if you need it \ref startSecondPhase(). |
404 | 405 |
|
405 | 406 |
///@{ |
406 | 407 |
|
407 | 408 |
/// \brief Initializes the internal data structures. |
408 | 409 |
/// |
409 | 410 |
/// Initializes the internal data structures and sets the initial |
410 | 411 |
/// flow to zero on each arc. |
411 | 412 |
void init() { |
412 | 413 |
createStructures(); |
413 | 414 |
|
414 | 415 |
_phase = true; |
415 | 416 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
416 | 417 |
(*_excess)[n] = 0; |
417 | 418 |
} |
418 | 419 |
|
419 | 420 |
for (ArcIt e(_graph); e != INVALID; ++e) { |
420 | 421 |
_flow->set(e, 0); |
421 | 422 |
} |
422 | 423 |
|
423 | 424 |
typename Digraph::template NodeMap<bool> reached(_graph, false); |
424 | 425 |
|
425 | 426 |
_level->initStart(); |
426 | 427 |
_level->initAddItem(_target); |
427 | 428 |
|
428 | 429 |
std::vector<Node> queue; |
429 | 430 |
reached[_source] = true; |
430 | 431 |
|
431 | 432 |
queue.push_back(_target); |
432 | 433 |
reached[_target] = true; |
433 | 434 |
while (!queue.empty()) { |
434 | 435 |
_level->initNewLevel(); |
435 | 436 |
std::vector<Node> nqueue; |
436 | 437 |
for (int i = 0; i < int(queue.size()); ++i) { |
437 | 438 |
Node n = queue[i]; |
438 | 439 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
439 | 440 |
Node u = _graph.source(e); |
440 | 441 |
if (!reached[u] && _tolerance.positive((*_capacity)[e])) { |
441 | 442 |
reached[u] = true; |
442 | 443 |
_level->initAddItem(u); |
443 | 444 |
nqueue.push_back(u); |
444 | 445 |
} |
445 | 446 |
} |
446 | 447 |
} |
447 | 448 |
queue.swap(nqueue); |
448 | 449 |
} |
449 | 450 |
_level->initFinish(); |
450 | 451 |
|
451 | 452 |
for (OutArcIt e(_graph, _source); e != INVALID; ++e) { |
452 | 453 |
if (_tolerance.positive((*_capacity)[e])) { |
453 | 454 |
Node u = _graph.target(e); |
454 | 455 |
if ((*_level)[u] == _level->maxLevel()) continue; |
455 | 456 |
_flow->set(e, (*_capacity)[e]); |
456 | 457 |
(*_excess)[u] += (*_capacity)[e]; |
457 | 458 |
if (u != _target && !_level->active(u)) { |
458 | 459 |
_level->activate(u); |
459 | 460 |
} |
460 | 461 |
} |
461 | 462 |
} |
462 | 463 |
} |
463 | 464 |
|
464 | 465 |
/// \brief Initializes the internal data structures using the |
465 | 466 |
/// given flow map. |
466 | 467 |
/// |
467 | 468 |
/// Initializes the internal data structures and sets the initial |
468 | 469 |
/// flow to the given \c flowMap. The \c flowMap should contain a |
469 | 470 |
/// flow or at least a preflow, i.e. at each node excluding the |
470 | 471 |
/// source node the incoming flow should greater or equal to the |
471 | 472 |
/// outgoing flow. |
472 | 473 |
/// \return \c false if the given \c flowMap is not a preflow. |
473 | 474 |
template <typename FlowMap> |
474 | 475 |
bool init(const FlowMap& flowMap) { |
475 | 476 |
createStructures(); |
476 | 477 |
|
477 | 478 |
for (ArcIt e(_graph); e != INVALID; ++e) { |
478 | 479 |
_flow->set(e, flowMap[e]); |
479 | 480 |
} |
480 | 481 |
|
481 | 482 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
482 | 483 |
Value excess = 0; |
483 | 484 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
484 | 485 |
excess += (*_flow)[e]; |
485 | 486 |
} |
486 | 487 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
487 | 488 |
excess -= (*_flow)[e]; |
488 | 489 |
} |
489 | 490 |
if (excess < 0 && n != _source) return false; |
490 | 491 |
(*_excess)[n] = excess; |
491 | 492 |
} |
492 | 493 |
|
493 | 494 |
typename Digraph::template NodeMap<bool> reached(_graph, false); |
494 | 495 |
|
495 | 496 |
_level->initStart(); |
496 | 497 |
_level->initAddItem(_target); |
497 | 498 |
|
498 | 499 |
std::vector<Node> queue; |
499 | 500 |
reached[_source] = true; |
500 | 501 |
|
501 | 502 |
queue.push_back(_target); |
502 | 503 |
reached[_target] = true; |
503 | 504 |
while (!queue.empty()) { |
504 | 505 |
_level->initNewLevel(); |
505 | 506 |
std::vector<Node> nqueue; |
506 | 507 |
for (int i = 0; i < int(queue.size()); ++i) { |
507 | 508 |
Node n = queue[i]; |
508 | 509 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
509 | 510 |
Node u = _graph.source(e); |
510 | 511 |
if (!reached[u] && |
511 | 512 |
_tolerance.positive((*_capacity)[e] - (*_flow)[e])) { |
512 | 513 |
reached[u] = true; |
513 | 514 |
_level->initAddItem(u); |
514 | 515 |
nqueue.push_back(u); |
515 | 516 |
} |
516 | 517 |
} |
517 | 518 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
518 | 519 |
Node v = _graph.target(e); |
519 | 520 |
if (!reached[v] && _tolerance.positive((*_flow)[e])) { |
520 | 521 |
reached[v] = true; |
521 | 522 |
_level->initAddItem(v); |
522 | 523 |
nqueue.push_back(v); |
523 | 524 |
} |
524 | 525 |
} |
525 | 526 |
} |
526 | 527 |
queue.swap(nqueue); |
527 | 528 |
} |
528 | 529 |
_level->initFinish(); |
529 | 530 |
|
530 | 531 |
for (OutArcIt e(_graph, _source); e != INVALID; ++e) { |
531 | 532 |
Value rem = (*_capacity)[e] - (*_flow)[e]; |
532 | 533 |
if (_tolerance.positive(rem)) { |
533 | 534 |
Node u = _graph.target(e); |
534 | 535 |
if ((*_level)[u] == _level->maxLevel()) continue; |
535 | 536 |
_flow->set(e, (*_capacity)[e]); |
536 | 537 |
(*_excess)[u] += rem; |
537 | 538 |
if (u != _target && !_level->active(u)) { |
538 | 539 |
_level->activate(u); |
539 | 540 |
} |
540 | 541 |
} |
541 | 542 |
} |
542 | 543 |
for (InArcIt e(_graph, _source); e != INVALID; ++e) { |
543 | 544 |
Value rem = (*_flow)[e]; |
544 | 545 |
if (_tolerance.positive(rem)) { |
545 | 546 |
Node v = _graph.source(e); |
546 | 547 |
if ((*_level)[v] == _level->maxLevel()) continue; |
547 | 548 |
_flow->set(e, 0); |
548 | 549 |
(*_excess)[v] += rem; |
549 | 550 |
if (v != _target && !_level->active(v)) { |
550 | 551 |
_level->activate(v); |
551 | 552 |
} |
552 | 553 |
} |
553 | 554 |
} |
554 | 555 |
return true; |
555 | 556 |
} |
556 | 557 |
|
557 | 558 |
/// \brief Starts the first phase of the preflow algorithm. |
558 | 559 |
/// |
559 | 560 |
/// The preflow algorithm consists of two phases, this method runs |
560 | 561 |
/// the first phase. After the first phase the maximum flow value |
561 | 562 |
/// and a minimum value cut can already be computed, although a |
562 | 563 |
/// maximum flow is not yet obtained. So after calling this method |
563 | 564 |
/// \ref flowValue() returns the value of a maximum flow and \ref |
564 | 565 |
/// minCut() returns a minimum cut. |
565 | 566 |
/// \pre One of the \ref init() functions must be called before |
566 | 567 |
/// using this function. |
567 | 568 |
void startFirstPhase() { |
568 | 569 |
_phase = true; |
569 | 570 |
|
570 | 571 |
Node n = _level->highestActive(); |
571 | 572 |
int level = _level->highestActiveLevel(); |
572 | 573 |
while (n != INVALID) { |
573 | 574 |
int num = _node_num; |
574 | 575 |
|
575 | 576 |
while (num > 0 && n != INVALID) { |
576 | 577 |
Value excess = (*_excess)[n]; |
577 | 578 |
int new_level = _level->maxLevel(); |
578 | 579 |
|
579 | 580 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
580 | 581 |
Value rem = (*_capacity)[e] - (*_flow)[e]; |
581 | 582 |
if (!_tolerance.positive(rem)) continue; |
582 | 583 |
Node v = _graph.target(e); |
583 | 584 |
if ((*_level)[v] < level) { |
584 | 585 |
if (!_level->active(v) && v != _target) { |
585 | 586 |
_level->activate(v); |
586 | 587 |
} |
587 | 588 |
if (!_tolerance.less(rem, excess)) { |
588 | 589 |
_flow->set(e, (*_flow)[e] + excess); |
589 | 590 |
(*_excess)[v] += excess; |
590 | 591 |
excess = 0; |
591 | 592 |
goto no_more_push_1; |
592 | 593 |
} else { |
593 | 594 |
excess -= rem; |
594 | 595 |
(*_excess)[v] += rem; |
595 | 596 |
_flow->set(e, (*_capacity)[e]); |
596 | 597 |
} |
597 | 598 |
} else if (new_level > (*_level)[v]) { |
598 | 599 |
new_level = (*_level)[v]; |
599 | 600 |
} |
600 | 601 |
} |
601 | 602 |
|
602 | 603 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
603 | 604 |
Value rem = (*_flow)[e]; |
604 | 605 |
if (!_tolerance.positive(rem)) continue; |
605 | 606 |
Node v = _graph.source(e); |
606 | 607 |
if ((*_level)[v] < level) { |
607 | 608 |
if (!_level->active(v) && v != _target) { |
608 | 609 |
_level->activate(v); |
609 | 610 |
} |
610 | 611 |
if (!_tolerance.less(rem, excess)) { |
611 | 612 |
_flow->set(e, (*_flow)[e] - excess); |
612 | 613 |
(*_excess)[v] += excess; |
613 | 614 |
excess = 0; |
614 | 615 |
goto no_more_push_1; |
615 | 616 |
} else { |
616 | 617 |
excess -= rem; |
617 | 618 |
(*_excess)[v] += rem; |
618 | 619 |
_flow->set(e, 0); |
619 | 620 |
} |
620 | 621 |
} else if (new_level > (*_level)[v]) { |
621 | 622 |
new_level = (*_level)[v]; |
622 | 623 |
} |
623 | 624 |
} |
624 | 625 |
|
625 | 626 |
no_more_push_1: |
626 | 627 |
|
627 | 628 |
(*_excess)[n] = excess; |
628 | 629 |
|
629 | 630 |
if (excess != 0) { |
630 | 631 |
if (new_level + 1 < _level->maxLevel()) { |
631 | 632 |
_level->liftHighestActive(new_level + 1); |
632 | 633 |
} else { |
633 | 634 |
_level->liftHighestActiveToTop(); |
634 | 635 |
} |
635 | 636 |
if (_level->emptyLevel(level)) { |
636 | 637 |
_level->liftToTop(level); |
637 | 638 |
} |
638 | 639 |
} else { |
639 | 640 |
_level->deactivate(n); |
640 | 641 |
} |
641 | 642 |
|
642 | 643 |
n = _level->highestActive(); |
643 | 644 |
level = _level->highestActiveLevel(); |
644 | 645 |
--num; |
645 | 646 |
} |
646 | 647 |
|
647 | 648 |
num = _node_num * 20; |
648 | 649 |
while (num > 0 && n != INVALID) { |
649 | 650 |
Value excess = (*_excess)[n]; |
650 | 651 |
int new_level = _level->maxLevel(); |
651 | 652 |
|
652 | 653 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
653 | 654 |
Value rem = (*_capacity)[e] - (*_flow)[e]; |
654 | 655 |
if (!_tolerance.positive(rem)) continue; |
655 | 656 |
Node v = _graph.target(e); |
656 | 657 |
if ((*_level)[v] < level) { |
657 | 658 |
if (!_level->active(v) && v != _target) { |
658 | 659 |
_level->activate(v); |
659 | 660 |
} |
660 | 661 |
if (!_tolerance.less(rem, excess)) { |
661 | 662 |
_flow->set(e, (*_flow)[e] + excess); |
662 | 663 |
(*_excess)[v] += excess; |
663 | 664 |
excess = 0; |
664 | 665 |
goto no_more_push_2; |
665 | 666 |
} else { |
666 | 667 |
excess -= rem; |
667 | 668 |
(*_excess)[v] += rem; |
668 | 669 |
_flow->set(e, (*_capacity)[e]); |
669 | 670 |
} |
670 | 671 |
} else if (new_level > (*_level)[v]) { |
671 | 672 |
new_level = (*_level)[v]; |
672 | 673 |
} |
673 | 674 |
} |
674 | 675 |
|
675 | 676 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
676 | 677 |
Value rem = (*_flow)[e]; |
677 | 678 |
if (!_tolerance.positive(rem)) continue; |
678 | 679 |
Node v = _graph.source(e); |
679 | 680 |
if ((*_level)[v] < level) { |
680 | 681 |
if (!_level->active(v) && v != _target) { |
681 | 682 |
_level->activate(v); |
682 | 683 |
} |
683 | 684 |
if (!_tolerance.less(rem, excess)) { |
684 | 685 |
_flow->set(e, (*_flow)[e] - excess); |
685 | 686 |
(*_excess)[v] += excess; |
686 | 687 |
excess = 0; |
687 | 688 |
goto no_more_push_2; |
688 | 689 |
} else { |
689 | 690 |
excess -= rem; |
690 | 691 |
(*_excess)[v] += rem; |
691 | 692 |
_flow->set(e, 0); |
692 | 693 |
} |
693 | 694 |
} else if (new_level > (*_level)[v]) { |
694 | 695 |
new_level = (*_level)[v]; |
695 | 696 |
} |
696 | 697 |
} |
697 | 698 |
|
698 | 699 |
no_more_push_2: |
699 | 700 |
|
700 | 701 |
(*_excess)[n] = excess; |
701 | 702 |
|
702 | 703 |
if (excess != 0) { |
703 | 704 |
if (new_level + 1 < _level->maxLevel()) { |
704 | 705 |
_level->liftActiveOn(level, new_level + 1); |
705 | 706 |
} else { |
706 | 707 |
_level->liftActiveToTop(level); |
707 | 708 |
} |
708 | 709 |
if (_level->emptyLevel(level)) { |
709 | 710 |
_level->liftToTop(level); |
710 | 711 |
} |
711 | 712 |
} else { |
712 | 713 |
_level->deactivate(n); |
713 | 714 |
} |
714 | 715 |
|
715 | 716 |
while (level >= 0 && _level->activeFree(level)) { |
716 | 717 |
--level; |
717 | 718 |
} |
718 | 719 |
if (level == -1) { |
719 | 720 |
n = _level->highestActive(); |
720 | 721 |
level = _level->highestActiveLevel(); |
721 | 722 |
} else { |
722 | 723 |
n = _level->activeOn(level); |
723 | 724 |
} |
724 | 725 |
--num; |
725 | 726 |
} |
726 | 727 |
} |
727 | 728 |
} |
728 | 729 |
|
729 | 730 |
/// \brief Starts the second phase of the preflow algorithm. |
730 | 731 |
/// |
731 | 732 |
/// The preflow algorithm consists of two phases, this method runs |
732 | 733 |
/// the second phase. After calling one of the \ref init() functions |
733 | 734 |
/// and \ref startFirstPhase() and then \ref startSecondPhase(), |
734 | 735 |
/// \ref flowMap() returns a maximum flow, \ref flowValue() returns the |
735 | 736 |
/// value of a maximum flow, \ref minCut() returns a minimum cut |
736 | 737 |
/// \pre One of the \ref init() functions and \ref startFirstPhase() |
737 | 738 |
/// must be called before using this function. |
738 | 739 |
void startSecondPhase() { |
739 | 740 |
_phase = false; |
740 | 741 |
|
741 | 742 |
typename Digraph::template NodeMap<bool> reached(_graph); |
742 | 743 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
743 | 744 |
reached[n] = (*_level)[n] < _level->maxLevel(); |
744 | 745 |
} |
745 | 746 |
|
746 | 747 |
_level->initStart(); |
747 | 748 |
_level->initAddItem(_source); |
748 | 749 |
|
749 | 750 |
std::vector<Node> queue; |
750 | 751 |
queue.push_back(_source); |
751 | 752 |
reached[_source] = true; |
752 | 753 |
|
753 | 754 |
while (!queue.empty()) { |
754 | 755 |
_level->initNewLevel(); |
755 | 756 |
std::vector<Node> nqueue; |
756 | 757 |
for (int i = 0; i < int(queue.size()); ++i) { |
757 | 758 |
Node n = queue[i]; |
758 | 759 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
759 | 760 |
Node v = _graph.target(e); |
760 | 761 |
if (!reached[v] && _tolerance.positive((*_flow)[e])) { |
761 | 762 |
reached[v] = true; |
762 | 763 |
_level->initAddItem(v); |
763 | 764 |
nqueue.push_back(v); |
764 | 765 |
} |
765 | 766 |
} |
766 | 767 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
767 | 768 |
Node u = _graph.source(e); |
768 | 769 |
if (!reached[u] && |
769 | 770 |
_tolerance.positive((*_capacity)[e] - (*_flow)[e])) { |
770 | 771 |
reached[u] = true; |
771 | 772 |
_level->initAddItem(u); |
772 | 773 |
nqueue.push_back(u); |
773 | 774 |
} |
774 | 775 |
} |
775 | 776 |
} |
776 | 777 |
queue.swap(nqueue); |
777 | 778 |
} |
778 | 779 |
_level->initFinish(); |
779 | 780 |
|
780 | 781 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
781 | 782 |
if (!reached[n]) { |
782 | 783 |
_level->dirtyTopButOne(n); |
783 | 784 |
} else if ((*_excess)[n] > 0 && _target != n) { |
784 | 785 |
_level->activate(n); |
785 | 786 |
} |
786 | 787 |
} |
787 | 788 |
|
788 | 789 |
Node n; |
789 | 790 |
while ((n = _level->highestActive()) != INVALID) { |
790 | 791 |
Value excess = (*_excess)[n]; |
791 | 792 |
int level = _level->highestActiveLevel(); |
792 | 793 |
int new_level = _level->maxLevel(); |
793 | 794 |
|
794 | 795 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
795 | 796 |
Value rem = (*_capacity)[e] - (*_flow)[e]; |
796 | 797 |
if (!_tolerance.positive(rem)) continue; |
797 | 798 |
Node v = _graph.target(e); |
798 | 799 |
if ((*_level)[v] < level) { |
799 | 800 |
if (!_level->active(v) && v != _source) { |
800 | 801 |
_level->activate(v); |
801 | 802 |
} |
802 | 803 |
if (!_tolerance.less(rem, excess)) { |
803 | 804 |
_flow->set(e, (*_flow)[e] + excess); |
804 | 805 |
(*_excess)[v] += excess; |
805 | 806 |
excess = 0; |
806 | 807 |
goto no_more_push; |
807 | 808 |
} else { |
808 | 809 |
excess -= rem; |
809 | 810 |
(*_excess)[v] += rem; |
810 | 811 |
_flow->set(e, (*_capacity)[e]); |
811 | 812 |
} |
812 | 813 |
} else if (new_level > (*_level)[v]) { |
813 | 814 |
new_level = (*_level)[v]; |
814 | 815 |
} |
815 | 816 |
} |
816 | 817 |
|
817 | 818 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
818 | 819 |
Value rem = (*_flow)[e]; |
819 | 820 |
if (!_tolerance.positive(rem)) continue; |
820 | 821 |
Node v = _graph.source(e); |
821 | 822 |
if ((*_level)[v] < level) { |
822 | 823 |
if (!_level->active(v) && v != _source) { |
823 | 824 |
_level->activate(v); |
824 | 825 |
} |
825 | 826 |
if (!_tolerance.less(rem, excess)) { |
826 | 827 |
_flow->set(e, (*_flow)[e] - excess); |
827 | 828 |
(*_excess)[v] += excess; |
828 | 829 |
excess = 0; |
829 | 830 |
goto no_more_push; |
830 | 831 |
} else { |
831 | 832 |
excess -= rem; |
832 | 833 |
(*_excess)[v] += rem; |
833 | 834 |
_flow->set(e, 0); |
834 | 835 |
} |
835 | 836 |
} else if (new_level > (*_level)[v]) { |
836 | 837 |
new_level = (*_level)[v]; |
837 | 838 |
} |
838 | 839 |
} |
839 | 840 |
|
840 | 841 |
no_more_push: |
841 | 842 |
|
842 | 843 |
(*_excess)[n] = excess; |
843 | 844 |
|
844 | 845 |
if (excess != 0) { |
845 | 846 |
if (new_level + 1 < _level->maxLevel()) { |
846 | 847 |
_level->liftHighestActive(new_level + 1); |
847 | 848 |
} else { |
848 | 849 |
// Calculation error |
849 | 850 |
_level->liftHighestActiveToTop(); |
850 | 851 |
} |
851 | 852 |
if (_level->emptyLevel(level)) { |
852 | 853 |
// Calculation error |
853 | 854 |
_level->liftToTop(level); |
854 | 855 |
} |
855 | 856 |
} else { |
856 | 857 |
_level->deactivate(n); |
857 | 858 |
} |
858 | 859 |
|
859 | 860 |
} |
860 | 861 |
} |
861 | 862 |
|
862 | 863 |
/// \brief Runs the preflow algorithm. |
863 | 864 |
/// |
864 | 865 |
/// Runs the preflow algorithm. |
865 | 866 |
/// \note pf.run() is just a shortcut of the following code. |
866 | 867 |
/// \code |
867 | 868 |
/// pf.init(); |
868 | 869 |
/// pf.startFirstPhase(); |
869 | 870 |
/// pf.startSecondPhase(); |
870 | 871 |
/// \endcode |
871 | 872 |
void run() { |
872 | 873 |
init(); |
873 | 874 |
startFirstPhase(); |
874 | 875 |
startSecondPhase(); |
875 | 876 |
} |
876 | 877 |
|
877 | 878 |
/// \brief Runs the preflow algorithm to compute the minimum cut. |
878 | 879 |
/// |
879 | 880 |
/// Runs the preflow algorithm to compute the minimum cut. |
880 | 881 |
/// \note pf.runMinCut() is just a shortcut of the following code. |
881 | 882 |
/// \code |
882 | 883 |
/// pf.init(); |
883 | 884 |
/// pf.startFirstPhase(); |
884 | 885 |
/// \endcode |
885 | 886 |
void runMinCut() { |
886 | 887 |
init(); |
887 | 888 |
startFirstPhase(); |
888 | 889 |
} |
889 | 890 |
|
890 | 891 |
/// @} |
891 | 892 |
|
892 | 893 |
/// \name Query Functions |
893 | 894 |
/// The results of the preflow algorithm can be obtained using these |
894 | 895 |
/// functions.\n |
895 | 896 |
/// Either one of the \ref run() "run*()" functions or one of the |
896 | 897 |
/// \ref startFirstPhase() "start*()" functions should be called |
897 | 898 |
/// before using them. |
898 | 899 |
|
899 | 900 |
///@{ |
900 | 901 |
|
901 | 902 |
/// \brief Returns the value of the maximum flow. |
902 | 903 |
/// |
903 | 904 |
/// Returns the value of the maximum flow by returning the excess |
904 | 905 |
/// of the target node. This value equals to the value of |
905 | 906 |
/// the maximum flow already after the first phase of the algorithm. |
906 | 907 |
/// |
907 | 908 |
/// \pre Either \ref run() or \ref init() must be called before |
908 | 909 |
/// using this function. |
909 | 910 |
Value flowValue() const { |
910 | 911 |
return (*_excess)[_target]; |
911 | 912 |
} |
912 | 913 |
|
913 | 914 |
/// \brief Returns the flow value on the given arc. |
914 | 915 |
/// |
915 | 916 |
/// Returns the flow value on the given arc. This method can |
916 | 917 |
/// be called after the second phase of the algorithm. |
917 | 918 |
/// |
918 | 919 |
/// \pre Either \ref run() or \ref init() must be called before |
919 | 920 |
/// using this function. |
920 | 921 |
Value flow(const Arc& arc) const { |
921 | 922 |
return (*_flow)[arc]; |
922 | 923 |
} |
923 | 924 |
|
924 | 925 |
/// \brief Returns a const reference to the flow map. |
925 | 926 |
/// |
926 | 927 |
/// Returns a const reference to the arc map storing the found flow. |
927 | 928 |
/// This method can be called after the second phase of the algorithm. |
928 | 929 |
/// |
929 | 930 |
/// \pre Either \ref run() or \ref init() must be called before |
930 | 931 |
/// using this function. |
931 | 932 |
const FlowMap& flowMap() const { |
932 | 933 |
return *_flow; |
933 | 934 |
} |
934 | 935 |
|
935 | 936 |
/// \brief Returns \c true when the node is on the source side of the |
936 | 937 |
/// minimum cut. |
937 | 938 |
/// |
938 | 939 |
/// Returns true when the node is on the source side of the found |
939 | 940 |
/// minimum cut. This method can be called both after running \ref |
940 | 941 |
/// startFirstPhase() and \ref startSecondPhase(). |
941 | 942 |
/// |
942 | 943 |
/// \pre Either \ref run() or \ref init() must be called before |
943 | 944 |
/// using this function. |
944 | 945 |
bool minCut(const Node& node) const { |
945 | 946 |
return ((*_level)[node] == _level->maxLevel()) == _phase; |
946 | 947 |
} |
947 | 948 |
|
948 | 949 |
/// \brief Gives back a minimum value cut. |
949 | 950 |
/// |
950 | 951 |
/// Sets \c cutMap to the characteristic vector of a minimum value |
951 | 952 |
/// cut. \c cutMap should be a \ref concepts::WriteMap "writable" |
952 | 953 |
/// node map with \c bool (or convertible) value type. |
953 | 954 |
/// |
954 | 955 |
/// This method can be called both after running \ref startFirstPhase() |
955 | 956 |
/// and \ref startSecondPhase(). The result after the second phase |
956 | 957 |
/// could be slightly different if inexact computation is used. |
957 | 958 |
/// |
958 | 959 |
/// \note This function calls \ref minCut() for each node, so it runs in |
959 | 960 |
/// O(n) time. |
960 | 961 |
/// |
961 | 962 |
/// \pre Either \ref run() or \ref init() must be called before |
962 | 963 |
/// using this function. |
963 | 964 |
template <typename CutMap> |
964 | 965 |
void minCutMap(CutMap& cutMap) const { |
965 | 966 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
966 | 967 |
cutMap.set(n, minCut(n)); |
967 | 968 |
} |
968 | 969 |
} |
969 | 970 |
|
970 | 971 |
/// @} |
971 | 972 |
}; |
972 | 973 |
} |
973 | 974 |
|
974 | 975 |
#endif |
0 comments (0 inline)