0
3
0
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 | 319 |
\e breadth-first \e search (BFS) and \e depth-first \e search (DFS) |
320 | 320 |
\ref clrs01algorithms. |
321 | 321 |
*/ |
322 | 322 |
|
323 | 323 |
/** |
324 | 324 |
@defgroup shortest_path Shortest Path Algorithms |
325 | 325 |
@ingroup algs |
326 | 326 |
\brief Algorithms for finding shortest paths. |
327 | 327 |
|
328 | 328 |
This group contains the algorithms for finding shortest paths in digraphs |
329 | 329 |
\ref clrs01algorithms. |
330 | 330 |
|
331 | 331 |
- \ref Dijkstra algorithm for finding shortest paths from a source node |
332 | 332 |
when all arc lengths are non-negative. |
333 | 333 |
- \ref BellmanFord "Bellman-Ford" algorithm for finding shortest paths |
334 | 334 |
from a source node when arc lenghts can be either positive or negative, |
335 | 335 |
but the digraph should not contain directed cycles with negative total |
336 | 336 |
length. |
337 | 337 |
- \ref FloydWarshall "Floyd-Warshall" and \ref Johnson "Johnson" algorithms |
338 | 338 |
for solving the \e all-pairs \e shortest \e paths \e problem when arc |
339 | 339 |
lenghts can be either positive or negative, but the digraph should |
340 | 340 |
not contain directed cycles with negative total length. |
341 | 341 |
- \ref Suurballe A successive shortest path algorithm for finding |
342 | 342 |
arc-disjoint paths between two nodes having minimum total length. |
343 | 343 |
*/ |
344 | 344 |
|
345 | 345 |
/** |
346 | 346 |
@defgroup spantree Minimum Spanning Tree Algorithms |
347 | 347 |
@ingroup algs |
348 | 348 |
\brief Algorithms for finding minimum cost spanning trees and arborescences. |
349 | 349 |
|
350 | 350 |
This group contains the algorithms for finding minimum cost spanning |
351 | 351 |
trees and arborescences \ref clrs01algorithms. |
352 | 352 |
*/ |
353 | 353 |
|
354 | 354 |
/** |
355 | 355 |
@defgroup max_flow Maximum Flow Algorithms |
356 | 356 |
@ingroup algs |
357 | 357 |
\brief Algorithms for finding maximum flows. |
358 | 358 |
|
359 | 359 |
This group contains the algorithms for finding maximum flows and |
360 | 360 |
feasible circulations \ref clrs01algorithms, \ref amo93networkflows. |
361 | 361 |
|
362 | 362 |
The \e maximum \e flow \e problem is to find a flow of maximum value between |
363 | 363 |
a single source and a single target. Formally, there is a \f$G=(V,A)\f$ |
364 | 364 |
digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and |
365 | 365 |
\f$s, t \in V\f$ source and target nodes. |
366 | 366 |
A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the |
367 | 367 |
following optimization problem. |
368 | 368 |
|
369 | 369 |
\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f] |
370 | 370 |
\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu) |
371 | 371 |
\quad \forall u\in V\setminus\{s,t\} \f] |
372 | 372 |
\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f] |
373 | 373 |
|
374 | 374 |
LEMON contains several algorithms for solving maximum flow problems: |
375 | 375 |
- \ref EdmondsKarp Edmonds-Karp algorithm |
376 | 376 |
\ref edmondskarp72theoretical. |
377 | 377 |
- \ref Preflow Goldberg-Tarjan's preflow push-relabel algorithm |
378 | 378 |
\ref goldberg88newapproach. |
379 | 379 |
- \ref DinitzSleatorTarjan Dinitz's blocking flow algorithm with dynamic trees |
380 | 380 |
\ref dinic70algorithm, \ref sleator83dynamic. |
381 | 381 |
- \ref GoldbergTarjan !Preflow push-relabel algorithm with dynamic trees |
382 | 382 |
\ref goldberg88newapproach, \ref sleator83dynamic. |
383 | 383 |
|
384 | 384 |
In most cases the \ref Preflow algorithm provides the |
385 | 385 |
fastest method for computing a maximum flow. All implementations |
386 | 386 |
also provide functions to query the minimum cut, which is the dual |
387 | 387 |
problem of maximum flow. |
388 | 388 |
|
389 | 389 |
\ref Circulation is a preflow push-relabel algorithm implemented directly |
390 | 390 |
for finding feasible circulations, which is a somewhat different problem, |
391 | 391 |
but it is strongly related to maximum flow. |
392 | 392 |
For more information, see \ref Circulation. |
393 | 393 |
*/ |
394 | 394 |
|
395 | 395 |
/** |
396 | 396 |
@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms |
397 | 397 |
@ingroup algs |
398 | 398 |
|
399 | 399 |
\brief Algorithms for finding minimum cost flows and circulations. |
400 | 400 |
|
401 | 401 |
This group contains the algorithms for finding minimum cost flows and |
402 | 402 |
circulations \ref amo93networkflows. For more information about this |
403 | 403 |
problem and its dual solution, see \ref min_cost_flow |
404 | 404 |
"Minimum Cost Flow Problem". |
405 | 405 |
|
406 | 406 |
LEMON contains several algorithms for this problem. |
407 | 407 |
- \ref NetworkSimplex Primal Network Simplex algorithm with various |
408 | 408 |
pivot strategies \ref dantzig63linearprog, \ref kellyoneill91netsimplex. |
409 |
- \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on |
|
410 |
cost scaling \ref goldberg90approximation, \ref goldberg97efficient, |
|
409 |
- \ref CostScaling Cost Scaling algorithm based on push/augment and |
|
410 |
relabel operations \ref goldberg90approximation, \ref goldberg97efficient, |
|
411 | 411 |
\ref bunnagel98efficient. |
412 |
- \ref CapacityScaling Successive Shortest %Path algorithm with optional |
|
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. |
|
412 |
- \ref CapacityScaling Capacity Scaling algorithm based on the successive |
|
413 |
shortest path method \ref edmondskarp72theoretical. |
|
414 |
- \ref CycleCanceling Cycle-Canceling algorithms, two of which are |
|
415 |
strongly polynomial \ref klein67primal, \ref goldberg89cyclecanceling. |
|
418 | 416 |
|
419 | 417 |
In general NetworkSimplex is the most efficient implementation, |
420 | 418 |
but in special cases other algorithms could be faster. |
421 | 419 |
For example, if the total supply and/or capacities are rather small, |
422 | 420 |
CapacityScaling is usually the fastest algorithm (without effective scaling). |
423 | 421 |
*/ |
424 | 422 |
|
425 | 423 |
/** |
426 | 424 |
@defgroup min_cut Minimum Cut Algorithms |
427 | 425 |
@ingroup algs |
428 | 426 |
|
429 | 427 |
\brief Algorithms for finding minimum cut in graphs. |
430 | 428 |
|
431 | 429 |
This group contains the algorithms for finding minimum cut in graphs. |
432 | 430 |
|
433 | 431 |
The \e minimum \e cut \e problem is to find a non-empty and non-complete |
434 | 432 |
\f$X\f$ subset of the nodes with minimum overall capacity on |
435 | 433 |
outgoing arcs. Formally, there is a \f$G=(V,A)\f$ digraph, a |
436 | 434 |
\f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function. The minimum |
437 | 435 |
cut is the \f$X\f$ solution of the next optimization problem: |
438 | 436 |
|
439 | 437 |
\f[ \min_{X \subset V, X\not\in \{\emptyset, V\}} |
440 | 438 |
\sum_{uv\in A: u\in X, v\not\in X}cap(uv) \f] |
441 | 439 |
|
442 | 440 |
LEMON contains several algorithms related to minimum cut problems: |
443 | 441 |
|
444 | 442 |
- \ref HaoOrlin "Hao-Orlin algorithm" for calculating minimum cut |
445 | 443 |
in directed graphs. |
446 | 444 |
- \ref NagamochiIbaraki "Nagamochi-Ibaraki algorithm" for |
447 | 445 |
calculating minimum cut in undirected graphs. |
448 | 446 |
- \ref GomoryHu "Gomory-Hu tree computation" for calculating |
449 | 447 |
all-pairs minimum cut in undirected graphs. |
450 | 448 |
|
451 | 449 |
If you want to find minimum cut just between two distinict nodes, |
452 | 450 |
see the \ref max_flow "maximum flow problem". |
453 | 451 |
*/ |
454 | 452 |
|
455 | 453 |
/** |
456 | 454 |
@defgroup min_mean_cycle Minimum Mean Cycle Algorithms |
457 | 455 |
@ingroup algs |
458 | 456 |
\brief Algorithms for finding minimum mean cycles. |
459 | 457 |
|
460 | 458 |
This group contains the algorithms for finding minimum mean cycles |
461 | 459 |
\ref clrs01algorithms, \ref amo93networkflows. |
462 | 460 |
|
463 | 461 |
The \e minimum \e mean \e cycle \e problem is to find a directed cycle |
464 | 462 |
of minimum mean length (cost) in a digraph. |
465 | 463 |
The mean length of a cycle is the average length of its arcs, i.e. the |
466 | 464 |
ratio between the total length of the cycle and the number of arcs on it. |
467 | 465 |
|
468 | 466 |
This problem has an important connection to \e conservative \e length |
469 | 467 |
\e functions, too. A length function on the arcs of a digraph is called |
470 | 468 |
conservative if and only if there is no directed cycle of negative total |
471 | 469 |
length. For an arbitrary length function, the negative of the minimum |
472 | 470 |
cycle mean is the smallest \f$\epsilon\f$ value so that increasing the |
473 | 471 |
arc lengths uniformly by \f$\epsilon\f$ results in a conservative length |
474 | 472 |
function. |
475 | 473 |
|
476 | 474 |
LEMON contains three algorithms for solving the minimum mean cycle problem: |
477 | 475 |
- \ref Karp "Karp"'s original algorithm \ref amo93networkflows, |
478 | 476 |
\ref dasdan98minmeancycle. |
479 | 477 |
- \ref HartmannOrlin "Hartmann-Orlin"'s algorithm, which is an improved |
480 | 478 |
version of Karp's algorithm \ref dasdan98minmeancycle. |
481 | 479 |
- \ref Howard "Howard"'s policy iteration algorithm |
482 | 480 |
\ref dasdan98minmeancycle. |
483 | 481 |
|
484 | 482 |
In practice, the Howard algorithm proved to be by far the most efficient |
485 | 483 |
one, though the best known theoretical bound on its running time is |
486 | 484 |
exponential. |
487 | 485 |
Both Karp and HartmannOrlin algorithms run in time O(ne) and use space |
488 | 486 |
O(n<sup>2</sup>+e), but the latter one is typically faster due to the |
489 | 487 |
applied early termination scheme. |
490 | 488 |
*/ |
491 | 489 |
|
492 | 490 |
/** |
493 | 491 |
@defgroup matching Matching Algorithms |
494 | 492 |
@ingroup algs |
495 | 493 |
\brief Algorithms for finding matchings in graphs and bipartite graphs. |
496 | 494 |
|
497 | 495 |
This group contains the algorithms for calculating |
498 | 496 |
matchings in graphs and bipartite graphs. The general matching problem is |
499 | 497 |
finding a subset of the edges for which each node has at most one incident |
500 | 498 |
edge. |
501 | 499 |
|
502 | 500 |
There are several different algorithms for calculate matchings in |
503 | 501 |
graphs. The matching problems in bipartite graphs are generally |
504 | 502 |
easier than in general graphs. The goal of the matching optimization |
505 | 503 |
can be finding maximum cardinality, maximum weight or minimum cost |
506 | 504 |
matching. The search can be constrained to find perfect or |
507 | 505 |
maximum cardinality matching. |
508 | 506 |
|
509 | 507 |
The matching algorithms implemented in LEMON: |
510 | 508 |
- \ref MaxBipartiteMatching Hopcroft-Karp augmenting path algorithm |
511 | 509 |
for calculating maximum cardinality matching in bipartite graphs. |
512 | 510 |
- \ref PrBipartiteMatching Push-relabel algorithm |
513 | 511 |
for calculating maximum cardinality matching in bipartite graphs. |
514 | 512 |
- \ref MaxWeightedBipartiteMatching |
515 | 513 |
Successive shortest path algorithm for calculating maximum weighted |
516 | 514 |
matching and maximum weighted bipartite matching in bipartite graphs. |
517 | 515 |
- \ref MinCostMaxBipartiteMatching |
518 | 516 |
Successive shortest path algorithm for calculating minimum cost maximum |
519 | 517 |
matching in bipartite graphs. |
520 | 518 |
- \ref MaxMatching Edmond's blossom shrinking algorithm for calculating |
521 | 519 |
maximum cardinality matching in general graphs. |
522 | 520 |
- \ref MaxWeightedMatching Edmond's blossom shrinking algorithm for calculating |
523 | 521 |
maximum weighted matching in general graphs. |
524 | 522 |
- \ref MaxWeightedPerfectMatching |
525 | 523 |
Edmond's blossom shrinking algorithm for calculating maximum weighted |
526 | 524 |
perfect matching in general graphs. |
527 | 525 |
|
528 | 526 |
\image html bipartite_matching.png |
529 | 527 |
\image latex bipartite_matching.eps "Bipartite Matching" width=\textwidth |
530 | 528 |
*/ |
531 | 529 |
|
532 | 530 |
/** |
533 | 531 |
@defgroup graph_properties Connectivity and Other Graph Properties |
534 | 532 |
@ingroup algs |
535 | 533 |
\brief Algorithms for discovering the graph properties |
536 | 534 |
|
537 | 535 |
This group contains the algorithms for discovering the graph properties |
538 | 536 |
like connectivity, bipartiteness, euler property, simplicity etc. |
539 | 537 |
|
540 | 538 |
\image html connected_components.png |
541 | 539 |
\image latex connected_components.eps "Connected components" width=\textwidth |
542 | 540 |
*/ |
543 | 541 |
|
544 | 542 |
/** |
545 | 543 |
@defgroup planar Planarity Embedding and Drawing |
546 | 544 |
@ingroup algs |
547 | 545 |
\brief Algorithms for planarity checking, embedding and drawing |
548 | 546 |
|
549 | 547 |
This group contains the algorithms for planarity checking, |
550 | 548 |
embedding and drawing. |
551 | 549 |
|
552 | 550 |
\image html planar.png |
553 | 551 |
\image latex planar.eps "Plane graph" width=\textwidth |
554 | 552 |
*/ |
555 | 553 |
|
556 | 554 |
/** |
557 | 555 |
@defgroup approx Approximation Algorithms |
558 | 556 |
@ingroup algs |
559 | 557 |
\brief Approximation algorithms. |
560 | 558 |
|
561 | 559 |
This group contains the approximation and heuristic algorithms |
562 | 560 |
implemented in LEMON. |
563 | 561 |
*/ |
564 | 562 |
|
565 | 563 |
/** |
566 | 564 |
@defgroup auxalg Auxiliary Algorithms |
567 | 565 |
@ingroup algs |
568 | 566 |
\brief Auxiliary algorithms implemented in LEMON. |
569 | 567 |
|
570 | 568 |
This group contains some algorithms implemented in LEMON |
571 | 569 |
in order to make it easier to implement complex algorithms. |
572 | 570 |
*/ |
573 | 571 |
|
574 | 572 |
/** |
575 | 573 |
@defgroup gen_opt_group General Optimization Tools |
576 | 574 |
\brief This group contains some general optimization frameworks |
577 | 575 |
implemented in LEMON. |
578 | 576 |
|
579 | 577 |
This group contains some general optimization frameworks |
580 | 578 |
implemented in LEMON. |
581 | 579 |
*/ |
582 | 580 |
|
583 | 581 |
/** |
584 | 582 |
@defgroup lp_group LP and MIP Solvers |
585 | 583 |
@ingroup gen_opt_group |
586 | 584 |
\brief LP and MIP solver interfaces for LEMON. |
587 | 585 |
|
588 | 586 |
This group contains LP and MIP solver interfaces for LEMON. |
589 | 587 |
Various LP solvers could be used in the same manner with this |
590 | 588 |
high-level interface. |
591 | 589 |
|
592 | 590 |
The currently supported solvers are \ref glpk, \ref clp, \ref cbc, |
593 | 591 |
\ref cplex, \ref soplex. |
594 | 592 |
*/ |
595 | 593 |
|
596 | 594 |
/** |
597 | 595 |
@defgroup lp_utils Tools for Lp and Mip Solvers |
598 | 596 |
@ingroup lp_group |
599 | 597 |
\brief Helper tools to the Lp and Mip solvers. |
600 | 598 |
|
601 | 599 |
This group adds some helper tools to general optimization framework |
602 | 600 |
implemented in LEMON. |
603 | 601 |
*/ |
604 | 602 |
|
605 | 603 |
/** |
606 | 604 |
@defgroup metah Metaheuristics |
607 | 605 |
@ingroup gen_opt_group |
608 | 606 |
\brief Metaheuristics for LEMON library. |
609 | 607 |
|
610 | 608 |
This group contains some metaheuristic optimization tools. |
611 | 609 |
*/ |
612 | 610 |
|
613 | 611 |
/** |
614 | 612 |
@defgroup utils Tools and Utilities |
615 | 613 |
\brief Tools and utilities for programming in LEMON |
616 | 614 |
|
617 | 615 |
Tools and utilities for programming in LEMON. |
618 | 616 |
*/ |
619 | 617 |
|
620 | 618 |
/** |
621 | 619 |
@defgroup gutils Basic Graph Utilities |
622 | 620 |
@ingroup utils |
623 | 621 |
\brief Simple basic graph utilities. |
624 | 622 |
|
625 | 623 |
This group contains some simple basic graph utilities. |
626 | 624 |
*/ |
627 | 625 |
|
628 | 626 |
/** |
629 | 627 |
@defgroup misc Miscellaneous Tools |
630 | 628 |
@ingroup utils |
631 | 629 |
\brief Tools for development, debugging and testing. |
632 | 630 |
|
633 | 631 |
This group contains several useful tools for development, |
634 | 632 |
debugging and testing. |
635 | 633 |
*/ |
636 | 634 |
|
637 | 635 |
/** |
638 | 636 |
@defgroup timecount Time Measuring and Counting |
639 | 637 |
@ingroup misc |
640 | 638 |
\brief Simple tools for measuring the performance of algorithms. |
641 | 639 |
|
642 | 640 |
This group contains simple tools for measuring the performance |
643 | 641 |
of algorithms. |
644 | 642 |
*/ |
645 | 643 |
|
646 | 644 |
/** |
647 | 645 |
@defgroup exceptions Exceptions |
648 | 646 |
@ingroup utils |
649 | 647 |
\brief Exceptions defined in LEMON. |
650 | 648 |
|
651 | 649 |
This group contains the exceptions defined in LEMON. |
652 | 650 |
*/ |
653 | 651 |
|
654 | 652 |
/** |
655 | 653 |
@defgroup io_group Input-Output |
656 | 654 |
\brief Graph Input-Output methods |
657 | 655 |
|
658 | 656 |
This group contains the tools for importing and exporting graphs |
659 | 657 |
and graph related data. Now it supports the \ref lgf-format |
660 | 658 |
"LEMON Graph Format", the \c DIMACS format and the encapsulated |
661 | 659 |
postscript (EPS) format. |
662 | 660 |
*/ |
663 | 661 |
|
664 | 662 |
/** |
665 | 663 |
@defgroup lemon_io LEMON Graph Format |
666 | 664 |
@ingroup io_group |
667 | 665 |
\brief Reading and writing LEMON Graph Format. |
668 | 666 |
|
669 | 667 |
This group contains methods for reading and writing |
670 | 668 |
\ref lgf-format "LEMON Graph Format". |
671 | 669 |
*/ |
672 | 670 |
|
673 | 671 |
/** |
674 | 672 |
@defgroup eps_io Postscript Exporting |
675 | 673 |
@ingroup io_group |
676 | 674 |
\brief General \c EPS drawer and graph exporter |
677 | 675 |
|
678 | 676 |
This group contains general \c EPS drawing methods and special |
679 | 677 |
graph exporting tools. |
680 | 678 |
*/ |
681 | 679 |
|
682 | 680 |
/** |
683 | 681 |
@defgroup dimacs_group DIMACS Format |
684 | 682 |
@ingroup io_group |
685 | 683 |
\brief Read and write files in DIMACS format |
686 | 684 |
|
687 | 685 |
Tools to read a digraph from or write it to a file in DIMACS format data. |
688 | 686 |
*/ |
689 | 687 |
|
690 | 688 |
/** |
691 | 689 |
@defgroup nauty_group NAUTY Format |
692 | 690 |
@ingroup io_group |
693 | 691 |
\brief Read \e Nauty format |
694 | 692 |
|
695 | 693 |
Tool to read graphs from \e Nauty format data. |
696 | 694 |
*/ |
697 | 695 |
|
698 | 696 |
/** |
699 | 697 |
@defgroup concept Concepts |
700 | 698 |
\brief Skeleton classes and concept checking classes |
701 | 699 |
|
702 | 700 |
This group contains the data/algorithm skeletons and concept checking |
703 | 701 |
classes implemented in LEMON. |
704 | 702 |
|
705 | 703 |
The purpose of the classes in this group is fourfold. |
706 | 704 |
|
707 | 705 |
- These classes contain the documentations of the %concepts. In order |
708 | 706 |
to avoid document multiplications, an implementation of a concept |
709 | 707 |
simply refers to the corresponding concept class. |
710 | 708 |
|
711 | 709 |
- These classes declare every functions, <tt>typedef</tt>s etc. an |
712 | 710 |
implementation of the %concepts should provide, however completely |
713 | 711 |
without implementations and real data structures behind the |
714 | 712 |
interface. On the other hand they should provide nothing else. All |
715 | 713 |
the algorithms working on a data structure meeting a certain concept |
716 | 714 |
should compile with these classes. (Though it will not run properly, |
717 | 715 |
of course.) In this way it is easily to check if an algorithm |
718 | 716 |
doesn't use any extra feature of a certain implementation. |
719 | 717 |
|
720 | 718 |
- The concept descriptor classes also provide a <em>checker class</em> |
721 | 719 |
that makes it possible to check whether a certain implementation of a |
722 | 720 |
concept indeed provides all the required features. |
723 | 721 |
|
724 | 722 |
- Finally, They can serve as a skeleton of a new implementation of a concept. |
725 | 723 |
*/ |
726 | 724 |
|
727 | 725 |
/** |
728 | 726 |
@defgroup graph_concepts Graph Structure Concepts |
729 | 727 |
@ingroup concept |
730 | 728 |
\brief Skeleton and concept checking classes for graph structures |
731 | 729 |
|
732 | 730 |
This group contains the skeletons and concept checking classes of |
733 | 731 |
graph structures. |
734 | 732 |
*/ |
735 | 733 |
|
736 | 734 |
/** |
737 | 735 |
@defgroup map_concepts Map Concepts |
738 | 736 |
@ingroup concept |
739 | 737 |
\brief Skeleton and concept checking classes for maps |
740 | 738 |
|
741 | 739 |
This group contains the skeletons and concept checking classes of maps. |
742 | 740 |
*/ |
743 | 741 |
|
744 | 742 |
/** |
745 | 743 |
@defgroup tools Standalone Utility Applications |
746 | 744 |
|
747 | 745 |
Some utility applications are listed here. |
748 | 746 |
|
749 | 747 |
The standard compilation procedure (<tt>./configure;make</tt>) will compile |
750 | 748 |
them, as well. |
751 | 749 |
*/ |
752 | 750 |
|
753 | 751 |
/** |
754 | 752 |
\anchor demoprograms |
755 | 753 |
|
756 | 754 |
@defgroup demos Demo Programs |
757 | 755 |
|
758 | 756 |
Some demo programs are listed here. Their full source codes can be found in |
759 | 757 |
the \c demo subdirectory of the source tree. |
760 | 758 |
|
761 | 759 |
In order to compile them, use the <tt>make demo</tt> or the |
762 | 760 |
<tt>make check</tt> commands. |
763 | 761 |
*/ |
764 | 762 |
|
765 | 763 |
} |
1 | 1 |
/* -*- C++ -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2008 |
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_CAPACITY_SCALING_H |
20 | 20 |
#define LEMON_CAPACITY_SCALING_H |
21 | 21 |
|
22 | 22 |
/// \ingroup min_cost_flow_algs |
23 | 23 |
/// |
24 | 24 |
/// \file |
25 | 25 |
/// \brief Capacity Scaling algorithm for finding a minimum cost flow. |
26 | 26 |
|
27 | 27 |
#include <vector> |
28 | 28 |
#include <limits> |
29 | 29 |
#include <lemon/core.h> |
30 | 30 |
#include <lemon/bin_heap.h> |
31 | 31 |
|
32 | 32 |
namespace lemon { |
33 | 33 |
|
34 | 34 |
/// \brief Default traits class of CapacityScaling algorithm. |
35 | 35 |
/// |
36 | 36 |
/// Default traits class of CapacityScaling algorithm. |
37 | 37 |
/// \tparam GR Digraph type. |
38 | 38 |
/// \tparam V The number type used for flow amounts, capacity bounds |
39 | 39 |
/// and supply values. By default it is \c int. |
40 | 40 |
/// \tparam C The number type used for costs and potentials. |
41 | 41 |
/// By default it is the same as \c V. |
42 | 42 |
template <typename GR, typename V = int, typename C = V> |
43 | 43 |
struct CapacityScalingDefaultTraits |
44 | 44 |
{ |
45 | 45 |
/// The type of the digraph |
46 | 46 |
typedef GR Digraph; |
47 | 47 |
/// The type of the flow amounts, capacity bounds and supply values |
48 | 48 |
typedef V Value; |
49 | 49 |
/// The type of the arc costs |
50 | 50 |
typedef C Cost; |
51 | 51 |
|
52 | 52 |
/// \brief The type of the heap used for internal Dijkstra computations. |
53 | 53 |
/// |
54 | 54 |
/// The type of the heap used for internal Dijkstra computations. |
55 | 55 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept, |
56 | 56 |
/// its priority type must be \c Cost and its cross reference type |
57 | 57 |
/// must be \ref RangeMap "RangeMap<int>". |
58 | 58 |
typedef BinHeap<Cost, RangeMap<int> > Heap; |
59 | 59 |
}; |
60 | 60 |
|
61 | 61 |
/// \addtogroup min_cost_flow_algs |
62 | 62 |
/// @{ |
63 | 63 |
|
64 | 64 |
/// \brief Implementation of the Capacity Scaling algorithm for |
65 | 65 |
/// finding a \ref min_cost_flow "minimum cost flow". |
66 | 66 |
/// |
67 | 67 |
/// \ref CapacityScaling implements the capacity scaling version |
68 | 68 |
/// of the successive shortest path algorithm for finding a |
69 |
/// \ref min_cost_flow "minimum cost flow" |
|
69 |
/// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows, |
|
70 |
/// \ref edmondskarp72theoretical. It is an efficient dual |
|
70 | 71 |
/// solution method. |
71 | 72 |
/// |
72 | 73 |
/// Most of the parameters of the problem (except for the digraph) |
73 | 74 |
/// can be given using separate functions, and the algorithm can be |
74 | 75 |
/// executed using the \ref run() function. If some parameters are not |
75 | 76 |
/// specified, then default values will be used. |
76 | 77 |
/// |
77 | 78 |
/// \tparam GR The digraph type the algorithm runs on. |
78 | 79 |
/// \tparam V The number type used for flow amounts, capacity bounds |
79 | 80 |
/// and supply values in the algorithm. By default it is \c int. |
80 | 81 |
/// \tparam C The number type used for costs and potentials in the |
81 | 82 |
/// algorithm. By default it is the same as \c V. |
82 | 83 |
/// |
83 | 84 |
/// \warning Both number types must be signed and all input data must |
84 | 85 |
/// be integer. |
85 | 86 |
/// \warning This algorithm does not support negative costs for such |
86 | 87 |
/// arcs that have infinite upper bound. |
87 | 88 |
#ifdef DOXYGEN |
88 | 89 |
template <typename GR, typename V, typename C, typename TR> |
89 | 90 |
#else |
90 | 91 |
template < typename GR, typename V = int, typename C = V, |
91 | 92 |
typename TR = CapacityScalingDefaultTraits<GR, V, C> > |
92 | 93 |
#endif |
93 | 94 |
class CapacityScaling |
94 | 95 |
{ |
95 | 96 |
public: |
96 | 97 |
|
97 | 98 |
/// The type of the digraph |
98 | 99 |
typedef typename TR::Digraph Digraph; |
99 | 100 |
/// The type of the flow amounts, capacity bounds and supply values |
100 | 101 |
typedef typename TR::Value Value; |
101 | 102 |
/// The type of the arc costs |
102 | 103 |
typedef typename TR::Cost Cost; |
103 | 104 |
|
104 | 105 |
/// The type of the heap used for internal Dijkstra computations |
105 | 106 |
typedef typename TR::Heap Heap; |
106 | 107 |
|
107 | 108 |
/// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm |
108 | 109 |
typedef TR Traits; |
109 | 110 |
|
110 | 111 |
public: |
111 | 112 |
|
112 | 113 |
/// \brief Problem type constants for the \c run() function. |
113 | 114 |
/// |
114 | 115 |
/// Enum type containing the problem type constants that can be |
115 | 116 |
/// returned by the \ref run() function of the algorithm. |
116 | 117 |
enum ProblemType { |
117 | 118 |
/// The problem has no feasible solution (flow). |
118 | 119 |
INFEASIBLE, |
119 | 120 |
/// The problem has optimal solution (i.e. it is feasible and |
120 | 121 |
/// bounded), and the algorithm has found optimal flow and node |
121 | 122 |
/// potentials (primal and dual solutions). |
122 | 123 |
OPTIMAL, |
123 | 124 |
/// The digraph contains an arc of negative cost and infinite |
124 | 125 |
/// upper bound. It means that the objective function is unbounded |
125 | 126 |
/// on that arc, however, note that it could actually be bounded |
126 | 127 |
/// over the feasible flows, but this algroithm cannot handle |
127 | 128 |
/// these cases. |
128 | 129 |
UNBOUNDED |
129 | 130 |
}; |
130 | 131 |
|
131 | 132 |
private: |
132 | 133 |
|
133 | 134 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
134 | 135 |
|
135 | 136 |
typedef std::vector<int> IntVector; |
136 | 137 |
typedef std::vector<char> BoolVector; |
137 | 138 |
typedef std::vector<Value> ValueVector; |
138 | 139 |
typedef std::vector<Cost> CostVector; |
139 | 140 |
|
140 | 141 |
private: |
141 | 142 |
|
142 | 143 |
// Data related to the underlying digraph |
143 | 144 |
const GR &_graph; |
144 | 145 |
int _node_num; |
145 | 146 |
int _arc_num; |
146 | 147 |
int _res_arc_num; |
147 | 148 |
int _root; |
148 | 149 |
|
149 | 150 |
// Parameters of the problem |
150 | 151 |
bool _have_lower; |
151 | 152 |
Value _sum_supply; |
152 | 153 |
|
153 | 154 |
// Data structures for storing the digraph |
154 | 155 |
IntNodeMap _node_id; |
155 | 156 |
IntArcMap _arc_idf; |
156 | 157 |
IntArcMap _arc_idb; |
157 | 158 |
IntVector _first_out; |
158 | 159 |
BoolVector _forward; |
159 | 160 |
IntVector _source; |
160 | 161 |
IntVector _target; |
161 | 162 |
IntVector _reverse; |
162 | 163 |
|
163 | 164 |
// Node and arc data |
164 | 165 |
ValueVector _lower; |
165 | 166 |
ValueVector _upper; |
166 | 167 |
CostVector _cost; |
167 | 168 |
ValueVector _supply; |
168 | 169 |
|
169 | 170 |
ValueVector _res_cap; |
170 | 171 |
CostVector _pi; |
171 | 172 |
ValueVector _excess; |
172 | 173 |
IntVector _excess_nodes; |
173 | 174 |
IntVector _deficit_nodes; |
174 | 175 |
|
175 | 176 |
Value _delta; |
176 | 177 |
int _factor; |
177 | 178 |
IntVector _pred; |
178 | 179 |
|
179 | 180 |
public: |
180 | 181 |
|
181 | 182 |
/// \brief Constant for infinite upper bounds (capacities). |
182 | 183 |
/// |
183 | 184 |
/// Constant for infinite upper bounds (capacities). |
184 | 185 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
185 | 186 |
/// \c std::numeric_limits<Value>::max() otherwise. |
186 | 187 |
const Value INF; |
187 | 188 |
|
188 | 189 |
private: |
189 | 190 |
|
190 | 191 |
// Special implementation of the Dijkstra algorithm for finding |
191 | 192 |
// shortest paths in the residual network of the digraph with |
192 | 193 |
// respect to the reduced arc costs and modifying the node |
193 | 194 |
// potentials according to the found distance labels. |
194 | 195 |
class ResidualDijkstra |
195 | 196 |
{ |
196 | 197 |
private: |
197 | 198 |
|
198 | 199 |
int _node_num; |
199 | 200 |
bool _geq; |
200 | 201 |
const IntVector &_first_out; |
201 | 202 |
const IntVector &_target; |
202 | 203 |
const CostVector &_cost; |
203 | 204 |
const ValueVector &_res_cap; |
204 | 205 |
const ValueVector &_excess; |
205 | 206 |
CostVector &_pi; |
206 | 207 |
IntVector &_pred; |
207 | 208 |
|
208 | 209 |
IntVector _proc_nodes; |
209 | 210 |
CostVector _dist; |
210 | 211 |
|
211 | 212 |
public: |
212 | 213 |
|
213 | 214 |
ResidualDijkstra(CapacityScaling& cs) : |
214 | 215 |
_node_num(cs._node_num), _geq(cs._sum_supply < 0), |
215 | 216 |
_first_out(cs._first_out), _target(cs._target), _cost(cs._cost), |
216 | 217 |
_res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi), |
217 | 218 |
_pred(cs._pred), _dist(cs._node_num) |
218 | 219 |
{} |
219 | 220 |
|
220 | 221 |
int run(int s, Value delta = 1) { |
221 | 222 |
RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP); |
222 | 223 |
Heap heap(heap_cross_ref); |
223 | 224 |
heap.push(s, 0); |
224 | 225 |
_pred[s] = -1; |
225 | 226 |
_proc_nodes.clear(); |
226 | 227 |
|
227 | 228 |
// Process nodes |
228 | 229 |
while (!heap.empty() && _excess[heap.top()] > -delta) { |
229 | 230 |
int u = heap.top(), v; |
230 | 231 |
Cost d = heap.prio() + _pi[u], dn; |
231 | 232 |
_dist[u] = heap.prio(); |
232 | 233 |
_proc_nodes.push_back(u); |
233 | 234 |
heap.pop(); |
234 | 235 |
|
235 | 236 |
// Traverse outgoing residual arcs |
236 | 237 |
int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1; |
237 | 238 |
for (int a = _first_out[u]; a != last_out; ++a) { |
238 | 239 |
if (_res_cap[a] < delta) continue; |
239 | 240 |
v = _target[a]; |
240 | 241 |
switch (heap.state(v)) { |
241 | 242 |
case Heap::PRE_HEAP: |
242 | 243 |
heap.push(v, d + _cost[a] - _pi[v]); |
243 | 244 |
_pred[v] = a; |
244 | 245 |
break; |
245 | 246 |
case Heap::IN_HEAP: |
246 | 247 |
dn = d + _cost[a] - _pi[v]; |
247 | 248 |
if (dn < heap[v]) { |
248 | 249 |
heap.decrease(v, dn); |
249 | 250 |
_pred[v] = a; |
250 | 251 |
} |
251 | 252 |
break; |
252 | 253 |
case Heap::POST_HEAP: |
253 | 254 |
break; |
254 | 255 |
} |
255 | 256 |
} |
256 | 257 |
} |
257 | 258 |
if (heap.empty()) return -1; |
258 | 259 |
|
259 | 260 |
// Update potentials of processed nodes |
260 | 261 |
int t = heap.top(); |
261 | 262 |
Cost dt = heap.prio(); |
262 | 263 |
for (int i = 0; i < int(_proc_nodes.size()); ++i) { |
263 | 264 |
_pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt; |
264 | 265 |
} |
265 | 266 |
|
266 | 267 |
return t; |
267 | 268 |
} |
268 | 269 |
|
269 | 270 |
}; //class ResidualDijkstra |
270 | 271 |
|
271 | 272 |
public: |
272 | 273 |
|
273 | 274 |
/// \name Named Template Parameters |
274 | 275 |
/// @{ |
275 | 276 |
|
276 | 277 |
template <typename T> |
277 | 278 |
struct SetHeapTraits : public Traits { |
278 | 279 |
typedef T Heap; |
279 | 280 |
}; |
280 | 281 |
|
281 | 282 |
/// \brief \ref named-templ-param "Named parameter" for setting |
282 | 283 |
/// \c Heap type. |
283 | 284 |
/// |
284 | 285 |
/// \ref named-templ-param "Named parameter" for setting \c Heap |
285 | 286 |
/// type, which is used for internal Dijkstra computations. |
286 | 287 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept, |
287 | 288 |
/// its priority type must be \c Cost and its cross reference type |
288 | 289 |
/// must be \ref RangeMap "RangeMap<int>". |
289 | 290 |
template <typename T> |
290 | 291 |
struct SetHeap |
291 | 292 |
: public CapacityScaling<GR, V, C, SetHeapTraits<T> > { |
292 | 293 |
typedef CapacityScaling<GR, V, C, SetHeapTraits<T> > Create; |
293 | 294 |
}; |
294 | 295 |
|
295 | 296 |
/// @} |
296 | 297 |
|
297 | 298 |
public: |
298 | 299 |
|
299 | 300 |
/// \brief Constructor. |
300 | 301 |
/// |
301 | 302 |
/// The constructor of the class. |
302 | 303 |
/// |
303 | 304 |
/// \param graph The digraph the algorithm runs on. |
304 | 305 |
CapacityScaling(const GR& graph) : |
305 | 306 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
306 | 307 |
INF(std::numeric_limits<Value>::has_infinity ? |
307 | 308 |
std::numeric_limits<Value>::infinity() : |
308 | 309 |
std::numeric_limits<Value>::max()) |
309 | 310 |
{ |
310 | 311 |
// Check the number types |
311 | 312 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
312 | 313 |
"The flow type of CapacityScaling must be signed"); |
313 | 314 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
314 | 315 |
"The cost type of CapacityScaling must be signed"); |
315 | 316 |
|
316 | 317 |
// Resize vectors |
317 | 318 |
_node_num = countNodes(_graph); |
318 | 319 |
_arc_num = countArcs(_graph); |
319 | 320 |
_res_arc_num = 2 * (_arc_num + _node_num); |
320 | 321 |
_root = _node_num; |
321 | 322 |
++_node_num; |
322 | 323 |
|
323 | 324 |
_first_out.resize(_node_num + 1); |
324 | 325 |
_forward.resize(_res_arc_num); |
325 | 326 |
_source.resize(_res_arc_num); |
326 | 327 |
_target.resize(_res_arc_num); |
327 | 328 |
_reverse.resize(_res_arc_num); |
328 | 329 |
|
329 | 330 |
_lower.resize(_res_arc_num); |
330 | 331 |
_upper.resize(_res_arc_num); |
331 | 332 |
_cost.resize(_res_arc_num); |
332 | 333 |
_supply.resize(_node_num); |
333 | 334 |
|
334 | 335 |
_res_cap.resize(_res_arc_num); |
335 | 336 |
_pi.resize(_node_num); |
336 | 337 |
_excess.resize(_node_num); |
337 | 338 |
_pred.resize(_node_num); |
338 | 339 |
|
339 | 340 |
// Copy the graph |
340 | 341 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1; |
341 | 342 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
342 | 343 |
_node_id[n] = i; |
343 | 344 |
} |
344 | 345 |
i = 0; |
345 | 346 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
346 | 347 |
_first_out[i] = j; |
347 | 348 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
348 | 349 |
_arc_idf[a] = j; |
349 | 350 |
_forward[j] = true; |
350 | 351 |
_source[j] = i; |
351 | 352 |
_target[j] = _node_id[_graph.runningNode(a)]; |
352 | 353 |
} |
353 | 354 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
354 | 355 |
_arc_idb[a] = j; |
355 | 356 |
_forward[j] = false; |
356 | 357 |
_source[j] = i; |
357 | 358 |
_target[j] = _node_id[_graph.runningNode(a)]; |
358 | 359 |
} |
359 | 360 |
_forward[j] = false; |
360 | 361 |
_source[j] = i; |
361 | 362 |
_target[j] = _root; |
362 | 363 |
_reverse[j] = k; |
363 | 364 |
_forward[k] = true; |
364 | 365 |
_source[k] = _root; |
365 | 366 |
_target[k] = i; |
366 | 367 |
_reverse[k] = j; |
367 | 368 |
++j; ++k; |
368 | 369 |
} |
369 | 370 |
_first_out[i] = j; |
370 | 371 |
_first_out[_node_num] = k; |
371 | 372 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
372 | 373 |
int fi = _arc_idf[a]; |
373 | 374 |
int bi = _arc_idb[a]; |
374 | 375 |
_reverse[fi] = bi; |
375 | 376 |
_reverse[bi] = fi; |
376 | 377 |
} |
377 | 378 |
|
378 | 379 |
// Reset parameters |
379 | 380 |
reset(); |
380 | 381 |
} |
381 | 382 |
|
382 | 383 |
/// \name Parameters |
383 | 384 |
/// The parameters of the algorithm can be specified using these |
384 | 385 |
/// functions. |
385 | 386 |
|
386 | 387 |
/// @{ |
387 | 388 |
|
388 | 389 |
/// \brief Set the lower bounds on the arcs. |
389 | 390 |
/// |
390 | 391 |
/// This function sets the lower bounds on the arcs. |
391 | 392 |
/// If it is not used before calling \ref run(), the lower bounds |
392 | 393 |
/// will be set to zero on all arcs. |
393 | 394 |
/// |
394 | 395 |
/// \param map An arc map storing the lower bounds. |
395 | 396 |
/// Its \c Value type must be convertible to the \c Value type |
396 | 397 |
/// of the algorithm. |
397 | 398 |
/// |
398 | 399 |
/// \return <tt>(*this)</tt> |
399 | 400 |
template <typename LowerMap> |
400 | 401 |
CapacityScaling& lowerMap(const LowerMap& map) { |
401 | 402 |
_have_lower = true; |
402 | 403 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
403 | 404 |
_lower[_arc_idf[a]] = map[a]; |
404 | 405 |
_lower[_arc_idb[a]] = map[a]; |
405 | 406 |
} |
406 | 407 |
return *this; |
407 | 408 |
} |
408 | 409 |
|
409 | 410 |
/// \brief Set the upper bounds (capacities) on the arcs. |
410 | 411 |
/// |
411 | 412 |
/// This function sets the upper bounds (capacities) on the arcs. |
412 | 413 |
/// If it is not used before calling \ref run(), the upper bounds |
413 | 414 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
414 | 415 |
/// unbounded from above). |
415 | 416 |
/// |
416 | 417 |
/// \param map An arc map storing the upper bounds. |
417 | 418 |
/// Its \c Value type must be convertible to the \c Value type |
418 | 419 |
/// of the algorithm. |
419 | 420 |
/// |
420 | 421 |
/// \return <tt>(*this)</tt> |
421 | 422 |
template<typename UpperMap> |
422 | 423 |
CapacityScaling& upperMap(const UpperMap& map) { |
423 | 424 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
424 | 425 |
_upper[_arc_idf[a]] = map[a]; |
425 | 426 |
} |
426 | 427 |
return *this; |
427 | 428 |
} |
428 | 429 |
|
429 | 430 |
/// \brief Set the costs of the arcs. |
430 | 431 |
/// |
431 | 432 |
/// This function sets the costs of the arcs. |
432 | 433 |
/// If it is not used before calling \ref run(), the costs |
433 | 434 |
/// will be set to \c 1 on all arcs. |
434 | 435 |
/// |
435 | 436 |
/// \param map An arc map storing the costs. |
436 | 437 |
/// Its \c Value type must be convertible to the \c Cost type |
437 | 438 |
/// of the algorithm. |
438 | 439 |
/// |
439 | 440 |
/// \return <tt>(*this)</tt> |
440 | 441 |
template<typename CostMap> |
441 | 442 |
CapacityScaling& costMap(const CostMap& map) { |
442 | 443 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
443 | 444 |
_cost[_arc_idf[a]] = map[a]; |
444 | 445 |
_cost[_arc_idb[a]] = -map[a]; |
445 | 446 |
} |
446 | 447 |
return *this; |
447 | 448 |
} |
448 | 449 |
|
449 | 450 |
/// \brief Set the supply values of the nodes. |
450 | 451 |
/// |
451 | 452 |
/// This function sets the supply values of the nodes. |
452 | 453 |
/// If neither this function nor \ref stSupply() is used before |
453 | 454 |
/// calling \ref run(), the supply of each node will be set to zero. |
454 | 455 |
/// |
455 | 456 |
/// \param map A node map storing the supply values. |
456 | 457 |
/// Its \c Value type must be convertible to the \c Value type |
457 | 458 |
/// of the algorithm. |
458 | 459 |
/// |
459 | 460 |
/// \return <tt>(*this)</tt> |
460 | 461 |
template<typename SupplyMap> |
461 | 462 |
CapacityScaling& supplyMap(const SupplyMap& map) { |
462 | 463 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
463 | 464 |
_supply[_node_id[n]] = map[n]; |
464 | 465 |
} |
465 | 466 |
return *this; |
466 | 467 |
} |
467 | 468 |
|
468 | 469 |
/// \brief Set single source and target nodes and a supply value. |
469 | 470 |
/// |
470 | 471 |
/// This function sets a single source node and a single target node |
471 | 472 |
/// and the required flow value. |
472 | 473 |
/// If neither this function nor \ref supplyMap() is used before |
473 | 474 |
/// calling \ref run(), the supply of each node will be set to zero. |
474 | 475 |
/// |
475 | 476 |
/// Using this function has the same effect as using \ref supplyMap() |
476 | 477 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
477 | 478 |
/// assigned to \c t and all other nodes have zero supply value. |
478 | 479 |
/// |
479 | 480 |
/// \param s The source node. |
480 | 481 |
/// \param t The target node. |
481 | 482 |
/// \param k The required amount of flow from node \c s to node \c t |
482 | 483 |
/// (i.e. the supply of \c s and the demand of \c t). |
483 | 484 |
/// |
484 | 485 |
/// \return <tt>(*this)</tt> |
485 | 486 |
CapacityScaling& stSupply(const Node& s, const Node& t, Value k) { |
486 | 487 |
for (int i = 0; i != _node_num; ++i) { |
487 | 488 |
_supply[i] = 0; |
488 | 489 |
} |
489 | 490 |
_supply[_node_id[s]] = k; |
490 | 491 |
_supply[_node_id[t]] = -k; |
491 | 492 |
return *this; |
492 | 493 |
} |
493 | 494 |
|
494 | 495 |
/// @} |
495 | 496 |
|
496 | 497 |
/// \name Execution control |
497 | 498 |
/// The algorithm can be executed using \ref run(). |
498 | 499 |
|
499 | 500 |
/// @{ |
500 | 501 |
|
501 | 502 |
/// \brief Run the algorithm. |
502 | 503 |
/// |
503 | 504 |
/// This function runs the algorithm. |
504 | 505 |
/// The paramters can be specified using functions \ref lowerMap(), |
505 | 506 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
506 | 507 |
/// For example, |
507 | 508 |
/// \code |
508 | 509 |
/// CapacityScaling<ListDigraph> cs(graph); |
509 | 510 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
510 | 511 |
/// .supplyMap(sup).run(); |
511 | 512 |
/// \endcode |
512 | 513 |
/// |
513 | 514 |
/// This function can be called more than once. All the parameters |
514 | 515 |
/// that have been given are kept for the next call, unless |
515 | 516 |
/// \ref reset() is called, thus only the modified parameters |
516 | 517 |
/// have to be set again. See \ref reset() for examples. |
517 | 518 |
/// However, the underlying digraph must not be modified after this |
518 | 519 |
/// class have been constructed, since it copies and extends the graph. |
519 | 520 |
/// |
520 | 521 |
/// \param factor The capacity scaling factor. It must be larger than |
521 | 522 |
/// one to use scaling. If it is less or equal to one, then scaling |
522 | 523 |
/// will be disabled. |
523 | 524 |
/// |
524 | 525 |
/// \return \c INFEASIBLE if no feasible flow exists, |
525 | 526 |
/// \n \c OPTIMAL if the problem has optimal solution |
526 | 527 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
527 | 528 |
/// optimal flow and node potentials (primal and dual solutions), |
528 | 529 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
529 | 530 |
/// and infinite upper bound. It means that the objective function |
530 | 531 |
/// is unbounded on that arc, however, note that it could actually be |
531 | 532 |
/// bounded over the feasible flows, but this algroithm cannot handle |
532 | 533 |
/// these cases. |
533 | 534 |
/// |
534 | 535 |
/// \see ProblemType |
535 | 536 |
ProblemType run(int factor = 4) { |
536 | 537 |
_factor = factor; |
537 | 538 |
ProblemType pt = init(); |
538 | 539 |
if (pt != OPTIMAL) return pt; |
539 | 540 |
return start(); |
540 | 541 |
} |
541 | 542 |
|
542 | 543 |
/// \brief Reset all the parameters that have been given before. |
543 | 544 |
/// |
544 | 545 |
/// This function resets all the paramaters that have been given |
545 | 546 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
546 | 547 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
547 | 548 |
/// |
548 | 549 |
/// It is useful for multiple run() calls. If this function is not |
549 | 550 |
/// used, all the parameters given before are kept for the next |
550 | 551 |
/// \ref run() call. |
551 | 552 |
/// However, the underlying digraph must not be modified after this |
552 | 553 |
/// class have been constructed, since it copies and extends the graph. |
553 | 554 |
/// |
554 | 555 |
/// For example, |
555 | 556 |
/// \code |
556 | 557 |
/// CapacityScaling<ListDigraph> cs(graph); |
557 | 558 |
/// |
558 | 559 |
/// // First run |
559 | 560 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
560 | 561 |
/// .supplyMap(sup).run(); |
561 | 562 |
/// |
562 | 563 |
/// // Run again with modified cost map (reset() is not called, |
563 | 564 |
/// // so only the cost map have to be set again) |
564 | 565 |
/// cost[e] += 100; |
565 | 566 |
/// cs.costMap(cost).run(); |
566 | 567 |
/// |
567 | 568 |
/// // Run again from scratch using reset() |
568 | 569 |
/// // (the lower bounds will be set to zero on all arcs) |
569 | 570 |
/// cs.reset(); |
570 | 571 |
/// cs.upperMap(capacity).costMap(cost) |
571 | 572 |
/// .supplyMap(sup).run(); |
572 | 573 |
/// \endcode |
573 | 574 |
/// |
574 | 575 |
/// \return <tt>(*this)</tt> |
575 | 576 |
CapacityScaling& reset() { |
576 | 577 |
for (int i = 0; i != _node_num; ++i) { |
577 | 578 |
_supply[i] = 0; |
578 | 579 |
} |
579 | 580 |
for (int j = 0; j != _res_arc_num; ++j) { |
580 | 581 |
_lower[j] = 0; |
581 | 582 |
_upper[j] = INF; |
582 | 583 |
_cost[j] = _forward[j] ? 1 : -1; |
583 | 584 |
} |
584 | 585 |
_have_lower = false; |
585 | 586 |
return *this; |
586 | 587 |
} |
587 | 588 |
|
588 | 589 |
/// @} |
589 | 590 |
|
590 | 591 |
/// \name Query Functions |
591 | 592 |
/// The results of the algorithm can be obtained using these |
592 | 593 |
/// functions.\n |
593 | 594 |
/// The \ref run() function must be called before using them. |
594 | 595 |
|
595 | 596 |
/// @{ |
596 | 597 |
|
597 | 598 |
/// \brief Return the total cost of the found flow. |
598 | 599 |
/// |
599 | 600 |
/// This function returns the total cost of the found flow. |
600 | 601 |
/// Its complexity is O(e). |
601 | 602 |
/// |
602 | 603 |
/// \note The return type of the function can be specified as a |
603 | 604 |
/// template parameter. For example, |
604 | 605 |
/// \code |
605 | 606 |
/// cs.totalCost<double>(); |
606 | 607 |
/// \endcode |
607 | 608 |
/// It is useful if the total cost cannot be stored in the \c Cost |
608 | 609 |
/// type of the algorithm, which is the default return type of the |
609 | 610 |
/// function. |
610 | 611 |
/// |
611 | 612 |
/// \pre \ref run() must be called before using this function. |
612 | 613 |
template <typename Number> |
613 | 614 |
Number totalCost() const { |
614 | 615 |
Number c = 0; |
615 | 616 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
616 | 617 |
int i = _arc_idb[a]; |
617 | 618 |
c += static_cast<Number>(_res_cap[i]) * |
618 | 619 |
(-static_cast<Number>(_cost[i])); |
619 | 620 |
} |
620 | 621 |
return c; |
621 | 622 |
} |
622 | 623 |
|
623 | 624 |
#ifndef DOXYGEN |
624 | 625 |
Cost totalCost() const { |
625 | 626 |
return totalCost<Cost>(); |
626 | 627 |
} |
627 | 628 |
#endif |
628 | 629 |
|
629 | 630 |
/// \brief Return the flow on the given arc. |
630 | 631 |
/// |
631 | 632 |
/// This function returns the flow on the given arc. |
632 | 633 |
/// |
633 | 634 |
/// \pre \ref run() must be called before using this function. |
634 | 635 |
Value flow(const Arc& a) const { |
635 | 636 |
return _res_cap[_arc_idb[a]]; |
636 | 637 |
} |
637 | 638 |
|
638 | 639 |
/// \brief Return the flow map (the primal solution). |
639 | 640 |
/// |
640 | 641 |
/// This function copies the flow value on each arc into the given |
641 | 642 |
/// map. The \c Value type of the algorithm must be convertible to |
642 | 643 |
/// the \c Value type of the map. |
643 | 644 |
/// |
644 | 645 |
/// \pre \ref run() must be called before using this function. |
645 | 646 |
template <typename FlowMap> |
646 | 647 |
void flowMap(FlowMap &map) const { |
647 | 648 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
648 | 649 |
map.set(a, _res_cap[_arc_idb[a]]); |
649 | 650 |
} |
650 | 651 |
} |
651 | 652 |
|
652 | 653 |
/// \brief Return the potential (dual value) of the given node. |
653 | 654 |
/// |
654 | 655 |
/// This function returns the potential (dual value) of the |
655 | 656 |
/// given node. |
656 | 657 |
/// |
657 | 658 |
/// \pre \ref run() must be called before using this function. |
658 | 659 |
Cost potential(const Node& n) const { |
659 | 660 |
return _pi[_node_id[n]]; |
660 | 661 |
} |
661 | 662 |
|
662 | 663 |
/// \brief Return the potential map (the dual solution). |
663 | 664 |
/// |
664 | 665 |
/// This function copies the potential (dual value) of each node |
665 | 666 |
/// into the given map. |
666 | 667 |
/// The \c Cost type of the algorithm must be convertible to the |
667 | 668 |
/// \c Value type of the map. |
668 | 669 |
/// |
669 | 670 |
/// \pre \ref run() must be called before using this function. |
670 | 671 |
template <typename PotentialMap> |
671 | 672 |
void potentialMap(PotentialMap &map) const { |
672 | 673 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
673 | 674 |
map.set(n, _pi[_node_id[n]]); |
674 | 675 |
} |
675 | 676 |
} |
676 | 677 |
|
677 | 678 |
/// @} |
678 | 679 |
|
679 | 680 |
private: |
680 | 681 |
|
681 | 682 |
// Initialize the algorithm |
682 | 683 |
ProblemType init() { |
683 | 684 |
if (_node_num == 0) return INFEASIBLE; |
684 | 685 |
|
685 | 686 |
// Check the sum of supply values |
686 | 687 |
_sum_supply = 0; |
687 | 688 |
for (int i = 0; i != _root; ++i) { |
688 | 689 |
_sum_supply += _supply[i]; |
689 | 690 |
} |
690 | 691 |
if (_sum_supply > 0) return INFEASIBLE; |
691 | 692 |
|
692 | 693 |
// Initialize vectors |
693 | 694 |
for (int i = 0; i != _root; ++i) { |
694 | 695 |
_pi[i] = 0; |
695 | 696 |
_excess[i] = _supply[i]; |
696 | 697 |
} |
697 | 698 |
|
698 | 699 |
// Remove non-zero lower bounds |
699 | 700 |
const Value MAX = std::numeric_limits<Value>::max(); |
700 | 701 |
int last_out; |
701 | 702 |
if (_have_lower) { |
702 | 703 |
for (int i = 0; i != _root; ++i) { |
703 | 704 |
last_out = _first_out[i+1]; |
704 | 705 |
for (int j = _first_out[i]; j != last_out; ++j) { |
705 | 706 |
if (_forward[j]) { |
706 | 707 |
Value c = _lower[j]; |
707 | 708 |
if (c >= 0) { |
708 | 709 |
_res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF; |
709 | 710 |
} else { |
710 | 711 |
_res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF; |
711 | 712 |
} |
712 | 713 |
_excess[i] -= c; |
713 | 714 |
_excess[_target[j]] += c; |
714 | 715 |
} else { |
715 | 716 |
_res_cap[j] = 0; |
716 | 717 |
} |
717 | 718 |
} |
718 | 719 |
} |
719 | 720 |
} else { |
720 | 721 |
for (int j = 0; j != _res_arc_num; ++j) { |
721 | 722 |
_res_cap[j] = _forward[j] ? _upper[j] : 0; |
722 | 723 |
} |
723 | 724 |
} |
724 | 725 |
|
725 | 726 |
// Handle negative costs |
726 | 727 |
for (int i = 0; i != _root; ++i) { |
727 | 728 |
last_out = _first_out[i+1] - 1; |
728 | 729 |
for (int j = _first_out[i]; j != last_out; ++j) { |
729 | 730 |
Value rc = _res_cap[j]; |
730 | 731 |
if (_cost[j] < 0 && rc > 0) { |
731 | 732 |
if (rc >= MAX) return UNBOUNDED; |
732 | 733 |
_excess[i] -= rc; |
733 | 734 |
_excess[_target[j]] += rc; |
734 | 735 |
_res_cap[j] = 0; |
735 | 736 |
_res_cap[_reverse[j]] += rc; |
736 | 737 |
} |
737 | 738 |
} |
738 | 739 |
} |
739 | 740 |
|
740 | 741 |
// Handle GEQ supply type |
741 | 742 |
if (_sum_supply < 0) { |
742 | 743 |
_pi[_root] = 0; |
743 | 744 |
_excess[_root] = -_sum_supply; |
744 | 745 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
745 | 746 |
int ra = _reverse[a]; |
746 | 747 |
_res_cap[a] = -_sum_supply + 1; |
747 | 748 |
_res_cap[ra] = 0; |
748 | 749 |
_cost[a] = 0; |
749 | 750 |
_cost[ra] = 0; |
750 | 751 |
} |
751 | 752 |
} else { |
752 | 753 |
_pi[_root] = 0; |
753 | 754 |
_excess[_root] = 0; |
754 | 755 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
755 | 756 |
int ra = _reverse[a]; |
756 | 757 |
_res_cap[a] = 1; |
757 | 758 |
_res_cap[ra] = 0; |
758 | 759 |
_cost[a] = 0; |
759 | 760 |
_cost[ra] = 0; |
760 | 761 |
} |
761 | 762 |
} |
762 | 763 |
|
763 | 764 |
// Initialize delta value |
764 | 765 |
if (_factor > 1) { |
765 | 766 |
// With scaling |
766 | 767 |
Value max_sup = 0, max_dem = 0; |
767 | 768 |
for (int i = 0; i != _node_num; ++i) { |
768 | 769 |
Value ex = _excess[i]; |
769 | 770 |
if ( ex > max_sup) max_sup = ex; |
770 | 771 |
if (-ex > max_dem) max_dem = -ex; |
771 | 772 |
} |
772 | 773 |
Value max_cap = 0; |
773 | 774 |
for (int j = 0; j != _res_arc_num; ++j) { |
774 | 775 |
if (_res_cap[j] > max_cap) max_cap = _res_cap[j]; |
775 | 776 |
} |
776 | 777 |
max_sup = std::min(std::min(max_sup, max_dem), max_cap); |
777 | 778 |
for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ; |
778 | 779 |
} else { |
779 | 780 |
// Without scaling |
780 | 781 |
_delta = 1; |
781 | 782 |
} |
782 | 783 |
|
783 | 784 |
return OPTIMAL; |
784 | 785 |
} |
785 | 786 |
|
786 | 787 |
ProblemType start() { |
787 | 788 |
// Execute the algorithm |
788 | 789 |
ProblemType pt; |
789 | 790 |
if (_delta > 1) |
790 | 791 |
pt = startWithScaling(); |
791 | 792 |
else |
792 | 793 |
pt = startWithoutScaling(); |
793 | 794 |
|
794 | 795 |
// Handle non-zero lower bounds |
795 | 796 |
if (_have_lower) { |
796 | 797 |
int limit = _first_out[_root]; |
797 | 798 |
for (int j = 0; j != limit; ++j) { |
798 | 799 |
if (!_forward[j]) _res_cap[j] += _lower[j]; |
799 | 800 |
} |
800 | 801 |
} |
801 | 802 |
|
802 | 803 |
// Shift potentials if necessary |
803 | 804 |
Cost pr = _pi[_root]; |
804 | 805 |
if (_sum_supply < 0 || pr > 0) { |
805 | 806 |
for (int i = 0; i != _node_num; ++i) { |
806 | 807 |
_pi[i] -= pr; |
807 | 808 |
} |
808 | 809 |
} |
809 | 810 |
|
810 | 811 |
return pt; |
811 | 812 |
} |
812 | 813 |
|
813 | 814 |
// Execute the capacity scaling algorithm |
814 | 815 |
ProblemType startWithScaling() { |
815 | 816 |
// Perform capacity scaling phases |
816 | 817 |
int s, t; |
817 | 818 |
ResidualDijkstra _dijkstra(*this); |
818 | 819 |
while (true) { |
819 | 820 |
// Saturate all arcs not satisfying the optimality condition |
820 | 821 |
int last_out; |
821 | 822 |
for (int u = 0; u != _node_num; ++u) { |
822 | 823 |
last_out = _sum_supply < 0 ? |
823 | 824 |
_first_out[u+1] : _first_out[u+1] - 1; |
824 | 825 |
for (int a = _first_out[u]; a != last_out; ++a) { |
825 | 826 |
int v = _target[a]; |
826 | 827 |
Cost c = _cost[a] + _pi[u] - _pi[v]; |
827 | 828 |
Value rc = _res_cap[a]; |
828 | 829 |
if (c < 0 && rc >= _delta) { |
829 | 830 |
_excess[u] -= rc; |
830 | 831 |
_excess[v] += rc; |
831 | 832 |
_res_cap[a] = 0; |
832 | 833 |
_res_cap[_reverse[a]] += rc; |
833 | 834 |
} |
834 | 835 |
} |
835 | 836 |
} |
836 | 837 |
|
837 | 838 |
// Find excess nodes and deficit nodes |
838 | 839 |
_excess_nodes.clear(); |
839 | 840 |
_deficit_nodes.clear(); |
840 | 841 |
for (int u = 0; u != _node_num; ++u) { |
841 | 842 |
Value ex = _excess[u]; |
842 | 843 |
if (ex >= _delta) _excess_nodes.push_back(u); |
843 | 844 |
if (ex <= -_delta) _deficit_nodes.push_back(u); |
844 | 845 |
} |
845 | 846 |
int next_node = 0, next_def_node = 0; |
846 | 847 |
|
847 | 848 |
// Find augmenting shortest paths |
848 | 849 |
while (next_node < int(_excess_nodes.size())) { |
849 | 850 |
// Check deficit nodes |
850 | 851 |
if (_delta > 1) { |
851 | 852 |
bool delta_deficit = false; |
852 | 853 |
for ( ; next_def_node < int(_deficit_nodes.size()); |
853 | 854 |
++next_def_node ) { |
854 | 855 |
if (_excess[_deficit_nodes[next_def_node]] <= -_delta) { |
855 | 856 |
delta_deficit = true; |
856 | 857 |
break; |
857 | 858 |
} |
858 | 859 |
} |
859 | 860 |
if (!delta_deficit) break; |
860 | 861 |
} |
861 | 862 |
|
862 | 863 |
// Run Dijkstra in the residual network |
863 | 864 |
s = _excess_nodes[next_node]; |
864 | 865 |
if ((t = _dijkstra.run(s, _delta)) == -1) { |
865 | 866 |
if (_delta > 1) { |
866 | 867 |
++next_node; |
867 | 868 |
continue; |
868 | 869 |
} |
869 | 870 |
return INFEASIBLE; |
870 | 871 |
} |
871 | 872 |
|
872 | 873 |
// Augment along a shortest path from s to t |
873 | 874 |
Value d = std::min(_excess[s], -_excess[t]); |
874 | 875 |
int u = t; |
875 | 876 |
int a; |
876 | 877 |
if (d > _delta) { |
877 | 878 |
while ((a = _pred[u]) != -1) { |
878 | 879 |
if (_res_cap[a] < d) d = _res_cap[a]; |
879 | 880 |
u = _source[a]; |
880 | 881 |
} |
881 | 882 |
} |
882 | 883 |
u = t; |
883 | 884 |
while ((a = _pred[u]) != -1) { |
884 | 885 |
_res_cap[a] -= d; |
885 | 886 |
_res_cap[_reverse[a]] += d; |
886 | 887 |
u = _source[a]; |
887 | 888 |
} |
888 | 889 |
_excess[s] -= d; |
889 | 890 |
_excess[t] += d; |
890 | 891 |
|
891 | 892 |
if (_excess[s] < _delta) ++next_node; |
892 | 893 |
} |
893 | 894 |
|
894 | 895 |
if (_delta == 1) break; |
895 | 896 |
_delta = _delta <= _factor ? 1 : _delta / _factor; |
896 | 897 |
} |
897 | 898 |
|
898 | 899 |
return OPTIMAL; |
899 | 900 |
} |
900 | 901 |
|
901 | 902 |
// Execute the successive shortest path algorithm |
902 | 903 |
ProblemType startWithoutScaling() { |
903 | 904 |
// Find excess nodes |
904 | 905 |
_excess_nodes.clear(); |
905 | 906 |
for (int i = 0; i != _node_num; ++i) { |
906 | 907 |
if (_excess[i] > 0) _excess_nodes.push_back(i); |
907 | 908 |
} |
908 | 909 |
if (_excess_nodes.size() == 0) return OPTIMAL; |
909 | 910 |
int next_node = 0; |
910 | 911 |
|
911 | 912 |
// Find shortest paths |
912 | 913 |
int s, t; |
913 | 914 |
ResidualDijkstra _dijkstra(*this); |
914 | 915 |
while ( _excess[_excess_nodes[next_node]] > 0 || |
915 | 916 |
++next_node < int(_excess_nodes.size()) ) |
916 | 917 |
{ |
917 | 918 |
// Run Dijkstra in the residual network |
918 | 919 |
s = _excess_nodes[next_node]; |
919 | 920 |
if ((t = _dijkstra.run(s)) == -1) return INFEASIBLE; |
920 | 921 |
|
921 | 922 |
// Augment along a shortest path from s to t |
922 | 923 |
Value d = std::min(_excess[s], -_excess[t]); |
923 | 924 |
int u = t; |
924 | 925 |
int a; |
925 | 926 |
if (d > 1) { |
926 | 927 |
while ((a = _pred[u]) != -1) { |
927 | 928 |
if (_res_cap[a] < d) d = _res_cap[a]; |
928 | 929 |
u = _source[a]; |
929 | 930 |
} |
930 | 931 |
} |
931 | 932 |
u = t; |
932 | 933 |
while ((a = _pred[u]) != -1) { |
933 | 934 |
_res_cap[a] -= d; |
934 | 935 |
_res_cap[_reverse[a]] += d; |
935 | 936 |
u = _source[a]; |
936 | 937 |
} |
937 | 938 |
_excess[s] -= d; |
938 | 939 |
_excess[t] += d; |
939 | 940 |
} |
940 | 941 |
|
941 | 942 |
return OPTIMAL; |
942 | 943 |
} |
943 | 944 |
|
944 | 945 |
}; //class CapacityScaling |
945 | 946 |
|
946 | 947 |
///@} |
947 | 948 |
|
948 | 949 |
} //namespace lemon |
949 | 950 |
|
950 | 951 |
#endif //LEMON_CAPACITY_SCALING_H |
1 | 1 |
/* -*- C++ -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2008 |
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_COST_SCALING_H |
20 | 20 |
#define LEMON_COST_SCALING_H |
21 | 21 |
|
22 | 22 |
/// \ingroup min_cost_flow_algs |
23 | 23 |
/// \file |
24 | 24 |
/// \brief Cost scaling algorithm for finding a minimum cost flow. |
25 | 25 |
|
26 | 26 |
#include <vector> |
27 | 27 |
#include <deque> |
28 | 28 |
#include <limits> |
29 | 29 |
|
30 | 30 |
#include <lemon/core.h> |
31 | 31 |
#include <lemon/maps.h> |
32 | 32 |
#include <lemon/math.h> |
33 | 33 |
#include <lemon/static_graph.h> |
34 | 34 |
#include <lemon/circulation.h> |
35 | 35 |
#include <lemon/bellman_ford.h> |
36 | 36 |
|
37 | 37 |
namespace lemon { |
38 | 38 |
|
39 | 39 |
/// \brief Default traits class of CostScaling algorithm. |
40 | 40 |
/// |
41 | 41 |
/// Default traits class of CostScaling algorithm. |
42 | 42 |
/// \tparam GR Digraph type. |
43 | 43 |
/// \tparam V The number type used for flow amounts, capacity bounds |
44 | 44 |
/// and supply values. By default it is \c int. |
45 | 45 |
/// \tparam C The number type used for costs and potentials. |
46 | 46 |
/// By default it is the same as \c V. |
47 | 47 |
#ifdef DOXYGEN |
48 | 48 |
template <typename GR, typename V = int, typename C = V> |
49 | 49 |
#else |
50 | 50 |
template < typename GR, typename V = int, typename C = V, |
51 | 51 |
bool integer = std::numeric_limits<C>::is_integer > |
52 | 52 |
#endif |
53 | 53 |
struct CostScalingDefaultTraits |
54 | 54 |
{ |
55 | 55 |
/// The type of the digraph |
56 | 56 |
typedef GR Digraph; |
57 | 57 |
/// The type of the flow amounts, capacity bounds and supply values |
58 | 58 |
typedef V Value; |
59 | 59 |
/// The type of the arc costs |
60 | 60 |
typedef C Cost; |
61 | 61 |
|
62 | 62 |
/// \brief The large cost type used for internal computations |
63 | 63 |
/// |
64 | 64 |
/// The large cost type used for internal computations. |
65 | 65 |
/// It is \c long \c long if the \c Cost type is integer, |
66 | 66 |
/// otherwise it is \c double. |
67 | 67 |
/// \c Cost must be convertible to \c LargeCost. |
68 | 68 |
typedef double LargeCost; |
69 | 69 |
}; |
70 | 70 |
|
71 | 71 |
// Default traits class for integer cost types |
72 | 72 |
template <typename GR, typename V, typename C> |
73 | 73 |
struct CostScalingDefaultTraits<GR, V, C, true> |
74 | 74 |
{ |
75 | 75 |
typedef GR Digraph; |
76 | 76 |
typedef V Value; |
77 | 77 |
typedef C Cost; |
78 | 78 |
#ifdef LEMON_HAVE_LONG_LONG |
79 | 79 |
typedef long long LargeCost; |
80 | 80 |
#else |
81 | 81 |
typedef long LargeCost; |
82 | 82 |
#endif |
83 | 83 |
}; |
84 | 84 |
|
85 | 85 |
|
86 | 86 |
/// \addtogroup min_cost_flow_algs |
87 | 87 |
/// @{ |
88 | 88 |
|
89 | 89 |
/// \brief Implementation of the Cost Scaling algorithm for |
90 | 90 |
/// finding a \ref min_cost_flow "minimum cost flow". |
91 | 91 |
/// |
92 | 92 |
/// \ref CostScaling implements a cost scaling algorithm that performs |
93 |
/// push/augment and relabel operations for finding a minimum cost |
|
94 |
/// flow. It is an efficient primal-dual solution method, which |
|
93 |
/// push/augment and relabel operations for finding a \ref min_cost_flow |
|
94 |
/// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation, |
|
95 |
/// \ref goldberg97efficient, \ref bunnagel98efficient. |
|
96 |
/// It is a highly efficient primal-dual solution method, which |
|
95 | 97 |
/// can be viewed as the generalization of the \ref Preflow |
96 | 98 |
/// "preflow push-relabel" algorithm for the maximum flow problem. |
97 | 99 |
/// |
98 | 100 |
/// Most of the parameters of the problem (except for the digraph) |
99 | 101 |
/// can be given using separate functions, and the algorithm can be |
100 | 102 |
/// executed using the \ref run() function. If some parameters are not |
101 | 103 |
/// specified, then default values will be used. |
102 | 104 |
/// |
103 | 105 |
/// \tparam GR The digraph type the algorithm runs on. |
104 | 106 |
/// \tparam V The number type used for flow amounts, capacity bounds |
105 | 107 |
/// and supply values in the algorithm. By default it is \c int. |
106 | 108 |
/// \tparam C The number type used for costs and potentials in the |
107 | 109 |
/// algorithm. By default it is the same as \c V. |
108 | 110 |
/// |
109 | 111 |
/// \warning Both number types must be signed and all input data must |
110 | 112 |
/// be integer. |
111 | 113 |
/// \warning This algorithm does not support negative costs for such |
112 | 114 |
/// arcs that have infinite upper bound. |
113 | 115 |
/// |
114 | 116 |
/// \note %CostScaling provides three different internal methods, |
115 | 117 |
/// from which the most efficient one is used by default. |
116 | 118 |
/// For more information, see \ref Method. |
117 | 119 |
#ifdef DOXYGEN |
118 | 120 |
template <typename GR, typename V, typename C, typename TR> |
119 | 121 |
#else |
120 | 122 |
template < typename GR, typename V = int, typename C = V, |
121 | 123 |
typename TR = CostScalingDefaultTraits<GR, V, C> > |
122 | 124 |
#endif |
123 | 125 |
class CostScaling |
124 | 126 |
{ |
125 | 127 |
public: |
126 | 128 |
|
127 | 129 |
/// The type of the digraph |
128 | 130 |
typedef typename TR::Digraph Digraph; |
129 | 131 |
/// The type of the flow amounts, capacity bounds and supply values |
130 | 132 |
typedef typename TR::Value Value; |
131 | 133 |
/// The type of the arc costs |
132 | 134 |
typedef typename TR::Cost Cost; |
133 | 135 |
|
134 | 136 |
/// \brief The large cost type |
135 | 137 |
/// |
136 | 138 |
/// The large cost type used for internal computations. |
137 | 139 |
/// Using the \ref CostScalingDefaultTraits "default traits class", |
138 | 140 |
/// it is \c long \c long if the \c Cost type is integer, |
139 | 141 |
/// otherwise it is \c double. |
140 | 142 |
typedef typename TR::LargeCost LargeCost; |
141 | 143 |
|
142 | 144 |
/// The \ref CostScalingDefaultTraits "traits class" of the algorithm |
143 | 145 |
typedef TR Traits; |
144 | 146 |
|
145 | 147 |
public: |
146 | 148 |
|
147 | 149 |
/// \brief Problem type constants for the \c run() function. |
148 | 150 |
/// |
149 | 151 |
/// Enum type containing the problem type constants that can be |
150 | 152 |
/// returned by the \ref run() function of the algorithm. |
151 | 153 |
enum ProblemType { |
152 | 154 |
/// The problem has no feasible solution (flow). |
153 | 155 |
INFEASIBLE, |
154 | 156 |
/// The problem has optimal solution (i.e. it is feasible and |
155 | 157 |
/// bounded), and the algorithm has found optimal flow and node |
156 | 158 |
/// potentials (primal and dual solutions). |
157 | 159 |
OPTIMAL, |
158 | 160 |
/// The digraph contains an arc of negative cost and infinite |
159 | 161 |
/// upper bound. It means that the objective function is unbounded |
160 | 162 |
/// on that arc, however, note that it could actually be bounded |
161 | 163 |
/// over the feasible flows, but this algroithm cannot handle |
162 | 164 |
/// these cases. |
163 | 165 |
UNBOUNDED |
164 | 166 |
}; |
165 | 167 |
|
166 | 168 |
/// \brief Constants for selecting the internal method. |
167 | 169 |
/// |
168 | 170 |
/// Enum type containing constants for selecting the internal method |
169 | 171 |
/// for the \ref run() function. |
170 | 172 |
/// |
171 | 173 |
/// \ref CostScaling provides three internal methods that differ mainly |
172 | 174 |
/// in their base operations, which are used in conjunction with the |
173 | 175 |
/// relabel operation. |
174 | 176 |
/// By default, the so called \ref PARTIAL_AUGMENT |
175 | 177 |
/// "Partial Augment-Relabel" method is used, which proved to be |
176 | 178 |
/// the most efficient and the most robust on various test inputs. |
177 | 179 |
/// However, the other methods can be selected using the \ref run() |
178 | 180 |
/// function with the proper parameter. |
179 | 181 |
enum Method { |
180 | 182 |
/// Local push operations are used, i.e. flow is moved only on one |
181 | 183 |
/// admissible arc at once. |
182 | 184 |
PUSH, |
183 | 185 |
/// Augment operations are used, i.e. flow is moved on admissible |
184 | 186 |
/// paths from a node with excess to a node with deficit. |
185 | 187 |
AUGMENT, |
186 | 188 |
/// Partial augment operations are used, i.e. flow is moved on |
187 | 189 |
/// admissible paths started from a node with excess, but the |
188 | 190 |
/// lengths of these paths are limited. This method can be viewed |
189 | 191 |
/// as a combined version of the previous two operations. |
190 | 192 |
PARTIAL_AUGMENT |
191 | 193 |
}; |
192 | 194 |
|
193 | 195 |
private: |
194 | 196 |
|
195 | 197 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
196 | 198 |
|
197 | 199 |
typedef std::vector<int> IntVector; |
198 | 200 |
typedef std::vector<char> BoolVector; |
199 | 201 |
typedef std::vector<Value> ValueVector; |
200 | 202 |
typedef std::vector<Cost> CostVector; |
201 | 203 |
typedef std::vector<LargeCost> LargeCostVector; |
202 | 204 |
|
203 | 205 |
private: |
204 | 206 |
|
205 | 207 |
template <typename KT, typename VT> |
206 | 208 |
class VectorMap { |
207 | 209 |
public: |
208 | 210 |
typedef KT Key; |
209 | 211 |
typedef VT Value; |
210 | 212 |
|
211 | 213 |
VectorMap(std::vector<Value>& v) : _v(v) {} |
212 | 214 |
|
213 | 215 |
const Value& operator[](const Key& key) const { |
214 | 216 |
return _v[StaticDigraph::id(key)]; |
215 | 217 |
} |
216 | 218 |
|
217 | 219 |
Value& operator[](const Key& key) { |
218 | 220 |
return _v[StaticDigraph::id(key)]; |
219 | 221 |
} |
220 | 222 |
|
221 | 223 |
void set(const Key& key, const Value& val) { |
222 | 224 |
_v[StaticDigraph::id(key)] = val; |
223 | 225 |
} |
224 | 226 |
|
225 | 227 |
private: |
226 | 228 |
std::vector<Value>& _v; |
227 | 229 |
}; |
228 | 230 |
|
229 | 231 |
typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap; |
230 | 232 |
typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap; |
231 | 233 |
|
232 | 234 |
private: |
233 | 235 |
|
234 | 236 |
// Data related to the underlying digraph |
235 | 237 |
const GR &_graph; |
236 | 238 |
int _node_num; |
237 | 239 |
int _arc_num; |
238 | 240 |
int _res_node_num; |
239 | 241 |
int _res_arc_num; |
240 | 242 |
int _root; |
241 | 243 |
|
242 | 244 |
// Parameters of the problem |
243 | 245 |
bool _have_lower; |
244 | 246 |
Value _sum_supply; |
245 | 247 |
|
246 | 248 |
// Data structures for storing the digraph |
247 | 249 |
IntNodeMap _node_id; |
248 | 250 |
IntArcMap _arc_idf; |
249 | 251 |
IntArcMap _arc_idb; |
250 | 252 |
IntVector _first_out; |
251 | 253 |
BoolVector _forward; |
252 | 254 |
IntVector _source; |
253 | 255 |
IntVector _target; |
254 | 256 |
IntVector _reverse; |
255 | 257 |
|
256 | 258 |
// Node and arc data |
257 | 259 |
ValueVector _lower; |
258 | 260 |
ValueVector _upper; |
259 | 261 |
CostVector _scost; |
260 | 262 |
ValueVector _supply; |
261 | 263 |
|
262 | 264 |
ValueVector _res_cap; |
263 | 265 |
LargeCostVector _cost; |
264 | 266 |
LargeCostVector _pi; |
265 | 267 |
ValueVector _excess; |
266 | 268 |
IntVector _next_out; |
267 | 269 |
std::deque<int> _active_nodes; |
268 | 270 |
|
269 | 271 |
// Data for scaling |
270 | 272 |
LargeCost _epsilon; |
271 | 273 |
int _alpha; |
272 | 274 |
|
273 | 275 |
// Data for a StaticDigraph structure |
274 | 276 |
typedef std::pair<int, int> IntPair; |
275 | 277 |
StaticDigraph _sgr; |
276 | 278 |
std::vector<IntPair> _arc_vec; |
277 | 279 |
std::vector<LargeCost> _cost_vec; |
278 | 280 |
LargeCostArcMap _cost_map; |
279 | 281 |
LargeCostNodeMap _pi_map; |
280 | 282 |
|
281 | 283 |
public: |
282 | 284 |
|
283 | 285 |
/// \brief Constant for infinite upper bounds (capacities). |
284 | 286 |
/// |
285 | 287 |
/// Constant for infinite upper bounds (capacities). |
286 | 288 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
287 | 289 |
/// \c std::numeric_limits<Value>::max() otherwise. |
288 | 290 |
const Value INF; |
289 | 291 |
|
290 | 292 |
public: |
291 | 293 |
|
292 | 294 |
/// \name Named Template Parameters |
293 | 295 |
/// @{ |
294 | 296 |
|
295 | 297 |
template <typename T> |
296 | 298 |
struct SetLargeCostTraits : public Traits { |
297 | 299 |
typedef T LargeCost; |
298 | 300 |
}; |
299 | 301 |
|
300 | 302 |
/// \brief \ref named-templ-param "Named parameter" for setting |
301 | 303 |
/// \c LargeCost type. |
302 | 304 |
/// |
303 | 305 |
/// \ref named-templ-param "Named parameter" for setting \c LargeCost |
304 | 306 |
/// type, which is used for internal computations in the algorithm. |
305 | 307 |
/// \c Cost must be convertible to \c LargeCost. |
306 | 308 |
template <typename T> |
307 | 309 |
struct SetLargeCost |
308 | 310 |
: public CostScaling<GR, V, C, SetLargeCostTraits<T> > { |
309 | 311 |
typedef CostScaling<GR, V, C, SetLargeCostTraits<T> > Create; |
310 | 312 |
}; |
311 | 313 |
|
312 | 314 |
/// @} |
313 | 315 |
|
314 | 316 |
public: |
315 | 317 |
|
316 | 318 |
/// \brief Constructor. |
317 | 319 |
/// |
318 | 320 |
/// The constructor of the class. |
319 | 321 |
/// |
320 | 322 |
/// \param graph The digraph the algorithm runs on. |
321 | 323 |
CostScaling(const GR& graph) : |
322 | 324 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
323 | 325 |
_cost_map(_cost_vec), _pi_map(_pi), |
324 | 326 |
INF(std::numeric_limits<Value>::has_infinity ? |
325 | 327 |
std::numeric_limits<Value>::infinity() : |
326 | 328 |
std::numeric_limits<Value>::max()) |
327 | 329 |
{ |
328 | 330 |
// Check the number types |
329 | 331 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
330 | 332 |
"The flow type of CostScaling must be signed"); |
331 | 333 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
332 | 334 |
"The cost type of CostScaling must be signed"); |
333 | 335 |
|
334 | 336 |
// Resize vectors |
335 | 337 |
_node_num = countNodes(_graph); |
336 | 338 |
_arc_num = countArcs(_graph); |
337 | 339 |
_res_node_num = _node_num + 1; |
338 | 340 |
_res_arc_num = 2 * (_arc_num + _node_num); |
339 | 341 |
_root = _node_num; |
340 | 342 |
|
341 | 343 |
_first_out.resize(_res_node_num + 1); |
342 | 344 |
_forward.resize(_res_arc_num); |
343 | 345 |
_source.resize(_res_arc_num); |
344 | 346 |
_target.resize(_res_arc_num); |
345 | 347 |
_reverse.resize(_res_arc_num); |
346 | 348 |
|
347 | 349 |
_lower.resize(_res_arc_num); |
348 | 350 |
_upper.resize(_res_arc_num); |
349 | 351 |
_scost.resize(_res_arc_num); |
350 | 352 |
_supply.resize(_res_node_num); |
351 | 353 |
|
352 | 354 |
_res_cap.resize(_res_arc_num); |
353 | 355 |
_cost.resize(_res_arc_num); |
354 | 356 |
_pi.resize(_res_node_num); |
355 | 357 |
_excess.resize(_res_node_num); |
356 | 358 |
_next_out.resize(_res_node_num); |
357 | 359 |
|
358 | 360 |
_arc_vec.reserve(_res_arc_num); |
359 | 361 |
_cost_vec.reserve(_res_arc_num); |
360 | 362 |
|
361 | 363 |
// Copy the graph |
362 | 364 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num; |
363 | 365 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
364 | 366 |
_node_id[n] = i; |
365 | 367 |
} |
366 | 368 |
i = 0; |
367 | 369 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
368 | 370 |
_first_out[i] = j; |
369 | 371 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
370 | 372 |
_arc_idf[a] = j; |
371 | 373 |
_forward[j] = true; |
372 | 374 |
_source[j] = i; |
373 | 375 |
_target[j] = _node_id[_graph.runningNode(a)]; |
374 | 376 |
} |
375 | 377 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
376 | 378 |
_arc_idb[a] = j; |
377 | 379 |
_forward[j] = false; |
378 | 380 |
_source[j] = i; |
379 | 381 |
_target[j] = _node_id[_graph.runningNode(a)]; |
380 | 382 |
} |
381 | 383 |
_forward[j] = false; |
382 | 384 |
_source[j] = i; |
383 | 385 |
_target[j] = _root; |
384 | 386 |
_reverse[j] = k; |
385 | 387 |
_forward[k] = true; |
386 | 388 |
_source[k] = _root; |
387 | 389 |
_target[k] = i; |
388 | 390 |
_reverse[k] = j; |
389 | 391 |
++j; ++k; |
390 | 392 |
} |
391 | 393 |
_first_out[i] = j; |
392 | 394 |
_first_out[_res_node_num] = k; |
393 | 395 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
394 | 396 |
int fi = _arc_idf[a]; |
395 | 397 |
int bi = _arc_idb[a]; |
396 | 398 |
_reverse[fi] = bi; |
397 | 399 |
_reverse[bi] = fi; |
398 | 400 |
} |
399 | 401 |
|
400 | 402 |
// Reset parameters |
401 | 403 |
reset(); |
402 | 404 |
} |
403 | 405 |
|
404 | 406 |
/// \name Parameters |
405 | 407 |
/// The parameters of the algorithm can be specified using these |
406 | 408 |
/// functions. |
407 | 409 |
|
408 | 410 |
/// @{ |
409 | 411 |
|
410 | 412 |
/// \brief Set the lower bounds on the arcs. |
411 | 413 |
/// |
412 | 414 |
/// This function sets the lower bounds on the arcs. |
413 | 415 |
/// If it is not used before calling \ref run(), the lower bounds |
414 | 416 |
/// will be set to zero on all arcs. |
415 | 417 |
/// |
416 | 418 |
/// \param map An arc map storing the lower bounds. |
417 | 419 |
/// Its \c Value type must be convertible to the \c Value type |
418 | 420 |
/// of the algorithm. |
419 | 421 |
/// |
420 | 422 |
/// \return <tt>(*this)</tt> |
421 | 423 |
template <typename LowerMap> |
422 | 424 |
CostScaling& lowerMap(const LowerMap& map) { |
423 | 425 |
_have_lower = true; |
424 | 426 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
425 | 427 |
_lower[_arc_idf[a]] = map[a]; |
426 | 428 |
_lower[_arc_idb[a]] = map[a]; |
427 | 429 |
} |
428 | 430 |
return *this; |
429 | 431 |
} |
430 | 432 |
|
431 | 433 |
/// \brief Set the upper bounds (capacities) on the arcs. |
432 | 434 |
/// |
433 | 435 |
/// This function sets the upper bounds (capacities) on the arcs. |
434 | 436 |
/// If it is not used before calling \ref run(), the upper bounds |
435 | 437 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
436 | 438 |
/// unbounded from above). |
437 | 439 |
/// |
438 | 440 |
/// \param map An arc map storing the upper bounds. |
439 | 441 |
/// Its \c Value type must be convertible to the \c Value type |
440 | 442 |
/// of the algorithm. |
441 | 443 |
/// |
442 | 444 |
/// \return <tt>(*this)</tt> |
443 | 445 |
template<typename UpperMap> |
444 | 446 |
CostScaling& upperMap(const UpperMap& map) { |
445 | 447 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
446 | 448 |
_upper[_arc_idf[a]] = map[a]; |
447 | 449 |
} |
448 | 450 |
return *this; |
449 | 451 |
} |
450 | 452 |
|
451 | 453 |
/// \brief Set the costs of the arcs. |
452 | 454 |
/// |
453 | 455 |
/// This function sets the costs of the arcs. |
454 | 456 |
/// If it is not used before calling \ref run(), the costs |
455 | 457 |
/// will be set to \c 1 on all arcs. |
456 | 458 |
/// |
457 | 459 |
/// \param map An arc map storing the costs. |
458 | 460 |
/// Its \c Value type must be convertible to the \c Cost type |
459 | 461 |
/// of the algorithm. |
460 | 462 |
/// |
461 | 463 |
/// \return <tt>(*this)</tt> |
462 | 464 |
template<typename CostMap> |
463 | 465 |
CostScaling& costMap(const CostMap& map) { |
464 | 466 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
465 | 467 |
_scost[_arc_idf[a]] = map[a]; |
466 | 468 |
_scost[_arc_idb[a]] = -map[a]; |
467 | 469 |
} |
468 | 470 |
return *this; |
469 | 471 |
} |
470 | 472 |
|
471 | 473 |
/// \brief Set the supply values of the nodes. |
472 | 474 |
/// |
473 | 475 |
/// This function sets the supply values of the nodes. |
474 | 476 |
/// If neither this function nor \ref stSupply() is used before |
475 | 477 |
/// calling \ref run(), the supply of each node will be set to zero. |
476 | 478 |
/// |
477 | 479 |
/// \param map A node map storing the supply values. |
478 | 480 |
/// Its \c Value type must be convertible to the \c Value type |
479 | 481 |
/// of the algorithm. |
480 | 482 |
/// |
481 | 483 |
/// \return <tt>(*this)</tt> |
482 | 484 |
template<typename SupplyMap> |
483 | 485 |
CostScaling& supplyMap(const SupplyMap& map) { |
484 | 486 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
485 | 487 |
_supply[_node_id[n]] = map[n]; |
486 | 488 |
} |
487 | 489 |
return *this; |
488 | 490 |
} |
489 | 491 |
|
490 | 492 |
/// \brief Set single source and target nodes and a supply value. |
491 | 493 |
/// |
492 | 494 |
/// This function sets a single source node and a single target node |
493 | 495 |
/// and the required flow value. |
494 | 496 |
/// If neither this function nor \ref supplyMap() is used before |
495 | 497 |
/// calling \ref run(), the supply of each node will be set to zero. |
496 | 498 |
/// |
497 | 499 |
/// Using this function has the same effect as using \ref supplyMap() |
498 | 500 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
499 | 501 |
/// assigned to \c t and all other nodes have zero supply value. |
500 | 502 |
/// |
501 | 503 |
/// \param s The source node. |
502 | 504 |
/// \param t The target node. |
503 | 505 |
/// \param k The required amount of flow from node \c s to node \c t |
504 | 506 |
/// (i.e. the supply of \c s and the demand of \c t). |
505 | 507 |
/// |
506 | 508 |
/// \return <tt>(*this)</tt> |
507 | 509 |
CostScaling& stSupply(const Node& s, const Node& t, Value k) { |
508 | 510 |
for (int i = 0; i != _res_node_num; ++i) { |
509 | 511 |
_supply[i] = 0; |
510 | 512 |
} |
511 | 513 |
_supply[_node_id[s]] = k; |
512 | 514 |
_supply[_node_id[t]] = -k; |
513 | 515 |
return *this; |
514 | 516 |
} |
515 | 517 |
|
516 | 518 |
/// @} |
517 | 519 |
|
518 | 520 |
/// \name Execution control |
519 | 521 |
/// The algorithm can be executed using \ref run(). |
520 | 522 |
|
521 | 523 |
/// @{ |
522 | 524 |
|
523 | 525 |
/// \brief Run the algorithm. |
524 | 526 |
/// |
525 | 527 |
/// This function runs the algorithm. |
526 | 528 |
/// The paramters can be specified using functions \ref lowerMap(), |
527 | 529 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
528 | 530 |
/// For example, |
529 | 531 |
/// \code |
530 | 532 |
/// CostScaling<ListDigraph> cs(graph); |
531 | 533 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
532 | 534 |
/// .supplyMap(sup).run(); |
533 | 535 |
/// \endcode |
534 | 536 |
/// |
535 | 537 |
/// This function can be called more than once. All the parameters |
536 | 538 |
/// that have been given are kept for the next call, unless |
537 | 539 |
/// \ref reset() is called, thus only the modified parameters |
538 | 540 |
/// have to be set again. See \ref reset() for examples. |
539 | 541 |
/// However, the underlying digraph must not be modified after this |
540 | 542 |
/// class have been constructed, since it copies and extends the graph. |
541 | 543 |
/// |
542 | 544 |
/// \param method The internal method that will be used in the |
543 | 545 |
/// algorithm. For more information, see \ref Method. |
544 | 546 |
/// \param factor The cost scaling factor. It must be larger than one. |
545 | 547 |
/// |
546 | 548 |
/// \return \c INFEASIBLE if no feasible flow exists, |
547 | 549 |
/// \n \c OPTIMAL if the problem has optimal solution |
548 | 550 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
549 | 551 |
/// optimal flow and node potentials (primal and dual solutions), |
550 | 552 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
551 | 553 |
/// and infinite upper bound. It means that the objective function |
552 | 554 |
/// is unbounded on that arc, however, note that it could actually be |
553 | 555 |
/// bounded over the feasible flows, but this algroithm cannot handle |
554 | 556 |
/// these cases. |
555 | 557 |
/// |
556 | 558 |
/// \see ProblemType, Method |
557 | 559 |
ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) { |
558 | 560 |
_alpha = factor; |
559 | 561 |
ProblemType pt = init(); |
560 | 562 |
if (pt != OPTIMAL) return pt; |
561 | 563 |
start(method); |
562 | 564 |
return OPTIMAL; |
563 | 565 |
} |
564 | 566 |
|
565 | 567 |
/// \brief Reset all the parameters that have been given before. |
566 | 568 |
/// |
567 | 569 |
/// This function resets all the paramaters that have been given |
568 | 570 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
569 | 571 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
570 | 572 |
/// |
571 | 573 |
/// It is useful for multiple run() calls. If this function is not |
572 | 574 |
/// used, all the parameters given before are kept for the next |
573 | 575 |
/// \ref run() call. |
574 | 576 |
/// However, the underlying digraph must not be modified after this |
575 | 577 |
/// class have been constructed, since it copies and extends the graph. |
576 | 578 |
/// |
577 | 579 |
/// For example, |
578 | 580 |
/// \code |
579 | 581 |
/// CostScaling<ListDigraph> cs(graph); |
580 | 582 |
/// |
581 | 583 |
/// // First run |
582 | 584 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
583 | 585 |
/// .supplyMap(sup).run(); |
584 | 586 |
/// |
585 | 587 |
/// // Run again with modified cost map (reset() is not called, |
586 | 588 |
/// // so only the cost map have to be set again) |
587 | 589 |
/// cost[e] += 100; |
588 | 590 |
/// cs.costMap(cost).run(); |
589 | 591 |
/// |
590 | 592 |
/// // Run again from scratch using reset() |
591 | 593 |
/// // (the lower bounds will be set to zero on all arcs) |
592 | 594 |
/// cs.reset(); |
593 | 595 |
/// cs.upperMap(capacity).costMap(cost) |
594 | 596 |
/// .supplyMap(sup).run(); |
595 | 597 |
/// \endcode |
596 | 598 |
/// |
597 | 599 |
/// \return <tt>(*this)</tt> |
598 | 600 |
CostScaling& reset() { |
599 | 601 |
for (int i = 0; i != _res_node_num; ++i) { |
600 | 602 |
_supply[i] = 0; |
601 | 603 |
} |
602 | 604 |
int limit = _first_out[_root]; |
603 | 605 |
for (int j = 0; j != limit; ++j) { |
604 | 606 |
_lower[j] = 0; |
605 | 607 |
_upper[j] = INF; |
606 | 608 |
_scost[j] = _forward[j] ? 1 : -1; |
607 | 609 |
} |
608 | 610 |
for (int j = limit; j != _res_arc_num; ++j) { |
609 | 611 |
_lower[j] = 0; |
610 | 612 |
_upper[j] = INF; |
611 | 613 |
_scost[j] = 0; |
612 | 614 |
_scost[_reverse[j]] = 0; |
613 | 615 |
} |
614 | 616 |
_have_lower = false; |
615 | 617 |
return *this; |
616 | 618 |
} |
617 | 619 |
|
618 | 620 |
/// @} |
619 | 621 |
|
620 | 622 |
/// \name Query Functions |
621 | 623 |
/// The results of the algorithm can be obtained using these |
622 | 624 |
/// functions.\n |
623 | 625 |
/// The \ref run() function must be called before using them. |
624 | 626 |
|
625 | 627 |
/// @{ |
626 | 628 |
|
627 | 629 |
/// \brief Return the total cost of the found flow. |
628 | 630 |
/// |
629 | 631 |
/// This function returns the total cost of the found flow. |
630 | 632 |
/// Its complexity is O(e). |
631 | 633 |
/// |
632 | 634 |
/// \note The return type of the function can be specified as a |
633 | 635 |
/// template parameter. For example, |
634 | 636 |
/// \code |
635 | 637 |
/// cs.totalCost<double>(); |
636 | 638 |
/// \endcode |
637 | 639 |
/// It is useful if the total cost cannot be stored in the \c Cost |
638 | 640 |
/// type of the algorithm, which is the default return type of the |
639 | 641 |
/// function. |
640 | 642 |
/// |
641 | 643 |
/// \pre \ref run() must be called before using this function. |
642 | 644 |
template <typename Number> |
643 | 645 |
Number totalCost() const { |
644 | 646 |
Number c = 0; |
645 | 647 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
646 | 648 |
int i = _arc_idb[a]; |
647 | 649 |
c += static_cast<Number>(_res_cap[i]) * |
648 | 650 |
(-static_cast<Number>(_scost[i])); |
649 | 651 |
} |
650 | 652 |
return c; |
651 | 653 |
} |
652 | 654 |
|
653 | 655 |
#ifndef DOXYGEN |
654 | 656 |
Cost totalCost() const { |
655 | 657 |
return totalCost<Cost>(); |
656 | 658 |
} |
657 | 659 |
#endif |
658 | 660 |
|
659 | 661 |
/// \brief Return the flow on the given arc. |
660 | 662 |
/// |
661 | 663 |
/// This function returns the flow on the given arc. |
662 | 664 |
/// |
663 | 665 |
/// \pre \ref run() must be called before using this function. |
664 | 666 |
Value flow(const Arc& a) const { |
665 | 667 |
return _res_cap[_arc_idb[a]]; |
666 | 668 |
} |
667 | 669 |
|
668 | 670 |
/// \brief Return the flow map (the primal solution). |
669 | 671 |
/// |
670 | 672 |
/// This function copies the flow value on each arc into the given |
671 | 673 |
/// map. The \c Value type of the algorithm must be convertible to |
672 | 674 |
/// the \c Value type of the map. |
673 | 675 |
/// |
674 | 676 |
/// \pre \ref run() must be called before using this function. |
675 | 677 |
template <typename FlowMap> |
676 | 678 |
void flowMap(FlowMap &map) const { |
677 | 679 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
678 | 680 |
map.set(a, _res_cap[_arc_idb[a]]); |
679 | 681 |
} |
680 | 682 |
} |
681 | 683 |
|
682 | 684 |
/// \brief Return the potential (dual value) of the given node. |
683 | 685 |
/// |
684 | 686 |
/// This function returns the potential (dual value) of the |
685 | 687 |
/// given node. |
686 | 688 |
/// |
687 | 689 |
/// \pre \ref run() must be called before using this function. |
688 | 690 |
Cost potential(const Node& n) const { |
689 | 691 |
return static_cast<Cost>(_pi[_node_id[n]]); |
690 | 692 |
} |
691 | 693 |
|
692 | 694 |
/// \brief Return the potential map (the dual solution). |
693 | 695 |
/// |
694 | 696 |
/// This function copies the potential (dual value) of each node |
695 | 697 |
/// into the given map. |
696 | 698 |
/// The \c Cost type of the algorithm must be convertible to the |
697 | 699 |
/// \c Value type of the map. |
698 | 700 |
/// |
699 | 701 |
/// \pre \ref run() must be called before using this function. |
700 | 702 |
template <typename PotentialMap> |
701 | 703 |
void potentialMap(PotentialMap &map) const { |
702 | 704 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
703 | 705 |
map.set(n, static_cast<Cost>(_pi[_node_id[n]])); |
704 | 706 |
} |
705 | 707 |
} |
706 | 708 |
|
707 | 709 |
/// @} |
708 | 710 |
|
709 | 711 |
private: |
710 | 712 |
|
711 | 713 |
// Initialize the algorithm |
712 | 714 |
ProblemType init() { |
713 | 715 |
if (_res_node_num == 0) return INFEASIBLE; |
714 | 716 |
|
715 | 717 |
// Check the sum of supply values |
716 | 718 |
_sum_supply = 0; |
717 | 719 |
for (int i = 0; i != _root; ++i) { |
718 | 720 |
_sum_supply += _supply[i]; |
719 | 721 |
} |
720 | 722 |
if (_sum_supply > 0) return INFEASIBLE; |
721 | 723 |
|
722 | 724 |
|
723 | 725 |
// Initialize vectors |
724 | 726 |
for (int i = 0; i != _res_node_num; ++i) { |
725 | 727 |
_pi[i] = 0; |
726 | 728 |
_excess[i] = _supply[i]; |
727 | 729 |
} |
728 | 730 |
|
729 | 731 |
// Remove infinite upper bounds and check negative arcs |
730 | 732 |
const Value MAX = std::numeric_limits<Value>::max(); |
731 | 733 |
int last_out; |
732 | 734 |
if (_have_lower) { |
733 | 735 |
for (int i = 0; i != _root; ++i) { |
734 | 736 |
last_out = _first_out[i+1]; |
735 | 737 |
for (int j = _first_out[i]; j != last_out; ++j) { |
736 | 738 |
if (_forward[j]) { |
737 | 739 |
Value c = _scost[j] < 0 ? _upper[j] : _lower[j]; |
738 | 740 |
if (c >= MAX) return UNBOUNDED; |
739 | 741 |
_excess[i] -= c; |
740 | 742 |
_excess[_target[j]] += c; |
741 | 743 |
} |
742 | 744 |
} |
743 | 745 |
} |
744 | 746 |
} else { |
745 | 747 |
for (int i = 0; i != _root; ++i) { |
746 | 748 |
last_out = _first_out[i+1]; |
747 | 749 |
for (int j = _first_out[i]; j != last_out; ++j) { |
748 | 750 |
if (_forward[j] && _scost[j] < 0) { |
749 | 751 |
Value c = _upper[j]; |
750 | 752 |
if (c >= MAX) return UNBOUNDED; |
751 | 753 |
_excess[i] -= c; |
752 | 754 |
_excess[_target[j]] += c; |
753 | 755 |
} |
754 | 756 |
} |
755 | 757 |
} |
756 | 758 |
} |
757 | 759 |
Value ex, max_cap = 0; |
758 | 760 |
for (int i = 0; i != _res_node_num; ++i) { |
759 | 761 |
ex = _excess[i]; |
760 | 762 |
_excess[i] = 0; |
761 | 763 |
if (ex < 0) max_cap -= ex; |
762 | 764 |
} |
763 | 765 |
for (int j = 0; j != _res_arc_num; ++j) { |
764 | 766 |
if (_upper[j] >= MAX) _upper[j] = max_cap; |
765 | 767 |
} |
766 | 768 |
|
767 | 769 |
// Initialize the large cost vector and the epsilon parameter |
768 | 770 |
_epsilon = 0; |
769 | 771 |
LargeCost lc; |
770 | 772 |
for (int i = 0; i != _root; ++i) { |
771 | 773 |
last_out = _first_out[i+1]; |
772 | 774 |
for (int j = _first_out[i]; j != last_out; ++j) { |
773 | 775 |
lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha; |
774 | 776 |
_cost[j] = lc; |
775 | 777 |
if (lc > _epsilon) _epsilon = lc; |
776 | 778 |
} |
777 | 779 |
} |
778 | 780 |
_epsilon /= _alpha; |
779 | 781 |
|
780 | 782 |
// Initialize maps for Circulation and remove non-zero lower bounds |
781 | 783 |
ConstMap<Arc, Value> low(0); |
782 | 784 |
typedef typename Digraph::template ArcMap<Value> ValueArcMap; |
783 | 785 |
typedef typename Digraph::template NodeMap<Value> ValueNodeMap; |
784 | 786 |
ValueArcMap cap(_graph), flow(_graph); |
785 | 787 |
ValueNodeMap sup(_graph); |
786 | 788 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
787 | 789 |
sup[n] = _supply[_node_id[n]]; |
788 | 790 |
} |
789 | 791 |
if (_have_lower) { |
790 | 792 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
791 | 793 |
int j = _arc_idf[a]; |
792 | 794 |
Value c = _lower[j]; |
793 | 795 |
cap[a] = _upper[j] - c; |
794 | 796 |
sup[_graph.source(a)] -= c; |
795 | 797 |
sup[_graph.target(a)] += c; |
796 | 798 |
} |
797 | 799 |
} else { |
798 | 800 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
799 | 801 |
cap[a] = _upper[_arc_idf[a]]; |
800 | 802 |
} |
801 | 803 |
} |
802 | 804 |
|
803 | 805 |
// Find a feasible flow using Circulation |
804 | 806 |
Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap> |
805 | 807 |
circ(_graph, low, cap, sup); |
806 | 808 |
if (!circ.flowMap(flow).run()) return INFEASIBLE; |
807 | 809 |
|
808 | 810 |
// Set residual capacities and handle GEQ supply type |
809 | 811 |
if (_sum_supply < 0) { |
810 | 812 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
811 | 813 |
Value fa = flow[a]; |
812 | 814 |
_res_cap[_arc_idf[a]] = cap[a] - fa; |
813 | 815 |
_res_cap[_arc_idb[a]] = fa; |
814 | 816 |
sup[_graph.source(a)] -= fa; |
815 | 817 |
sup[_graph.target(a)] += fa; |
816 | 818 |
} |
817 | 819 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
818 | 820 |
_excess[_node_id[n]] = sup[n]; |
819 | 821 |
} |
820 | 822 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
821 | 823 |
int u = _target[a]; |
822 | 824 |
int ra = _reverse[a]; |
823 | 825 |
_res_cap[a] = -_sum_supply + 1; |
824 | 826 |
_res_cap[ra] = -_excess[u]; |
825 | 827 |
_cost[a] = 0; |
826 | 828 |
_cost[ra] = 0; |
827 | 829 |
_excess[u] = 0; |
828 | 830 |
} |
829 | 831 |
} else { |
830 | 832 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
831 | 833 |
Value fa = flow[a]; |
832 | 834 |
_res_cap[_arc_idf[a]] = cap[a] - fa; |
833 | 835 |
_res_cap[_arc_idb[a]] = fa; |
834 | 836 |
} |
835 | 837 |
for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
836 | 838 |
int ra = _reverse[a]; |
837 | 839 |
_res_cap[a] = 1; |
838 | 840 |
_res_cap[ra] = 0; |
839 | 841 |
_cost[a] = 0; |
840 | 842 |
_cost[ra] = 0; |
841 | 843 |
} |
842 | 844 |
} |
843 | 845 |
|
844 | 846 |
return OPTIMAL; |
845 | 847 |
} |
846 | 848 |
|
847 | 849 |
// Execute the algorithm and transform the results |
848 | 850 |
void start(Method method) { |
849 | 851 |
// Maximum path length for partial augment |
850 | 852 |
const int MAX_PATH_LENGTH = 4; |
851 | 853 |
|
852 | 854 |
// Execute the algorithm |
853 | 855 |
switch (method) { |
854 | 856 |
case PUSH: |
855 | 857 |
startPush(); |
856 | 858 |
break; |
857 | 859 |
case AUGMENT: |
858 | 860 |
startAugment(); |
859 | 861 |
break; |
860 | 862 |
case PARTIAL_AUGMENT: |
861 | 863 |
startAugment(MAX_PATH_LENGTH); |
862 | 864 |
break; |
863 | 865 |
} |
864 | 866 |
|
865 | 867 |
// Compute node potentials for the original costs |
866 | 868 |
_arc_vec.clear(); |
867 | 869 |
_cost_vec.clear(); |
868 | 870 |
for (int j = 0; j != _res_arc_num; ++j) { |
869 | 871 |
if (_res_cap[j] > 0) { |
870 | 872 |
_arc_vec.push_back(IntPair(_source[j], _target[j])); |
871 | 873 |
_cost_vec.push_back(_scost[j]); |
872 | 874 |
} |
873 | 875 |
} |
874 | 876 |
_sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); |
875 | 877 |
|
876 | 878 |
typename BellmanFord<StaticDigraph, LargeCostArcMap> |
877 | 879 |
::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map); |
878 | 880 |
bf.distMap(_pi_map); |
879 | 881 |
bf.init(0); |
880 | 882 |
bf.start(); |
881 | 883 |
|
882 | 884 |
// Handle non-zero lower bounds |
883 | 885 |
if (_have_lower) { |
884 | 886 |
int limit = _first_out[_root]; |
885 | 887 |
for (int j = 0; j != limit; ++j) { |
886 | 888 |
if (!_forward[j]) _res_cap[j] += _lower[j]; |
887 | 889 |
} |
888 | 890 |
} |
889 | 891 |
} |
890 | 892 |
|
891 | 893 |
/// Execute the algorithm performing augment and relabel operations |
892 | 894 |
void startAugment(int max_length = std::numeric_limits<int>::max()) { |
893 | 895 |
// Paramters for heuristics |
894 | 896 |
const int BF_HEURISTIC_EPSILON_BOUND = 1000; |
895 | 897 |
const int BF_HEURISTIC_BOUND_FACTOR = 3; |
896 | 898 |
|
897 | 899 |
// Perform cost scaling phases |
898 | 900 |
IntVector pred_arc(_res_node_num); |
899 | 901 |
std::vector<int> path_nodes; |
900 | 902 |
for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
901 | 903 |
1 : _epsilon / _alpha ) |
902 | 904 |
{ |
903 | 905 |
// "Early Termination" heuristic: use Bellman-Ford algorithm |
904 | 906 |
// to check if the current flow is optimal |
905 | 907 |
if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { |
906 | 908 |
_arc_vec.clear(); |
907 | 909 |
_cost_vec.clear(); |
908 | 910 |
for (int j = 0; j != _res_arc_num; ++j) { |
909 | 911 |
if (_res_cap[j] > 0) { |
910 | 912 |
_arc_vec.push_back(IntPair(_source[j], _target[j])); |
911 | 913 |
_cost_vec.push_back(_cost[j] + 1); |
912 | 914 |
} |
913 | 915 |
} |
914 | 916 |
_sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); |
915 | 917 |
|
916 | 918 |
BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map); |
917 | 919 |
bf.init(0); |
918 | 920 |
bool done = false; |
919 | 921 |
int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num)); |
920 | 922 |
for (int i = 0; i < K && !done; ++i) |
921 | 923 |
done = bf.processNextWeakRound(); |
922 | 924 |
if (done) break; |
923 | 925 |
} |
924 | 926 |
|
925 | 927 |
// Saturate arcs not satisfying the optimality condition |
926 | 928 |
for (int a = 0; a != _res_arc_num; ++a) { |
927 | 929 |
if (_res_cap[a] > 0 && |
928 | 930 |
_cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { |
929 | 931 |
Value delta = _res_cap[a]; |
930 | 932 |
_excess[_source[a]] -= delta; |
931 | 933 |
_excess[_target[a]] += delta; |
932 | 934 |
_res_cap[a] = 0; |
933 | 935 |
_res_cap[_reverse[a]] += delta; |
934 | 936 |
} |
935 | 937 |
} |
936 | 938 |
|
937 | 939 |
// Find active nodes (i.e. nodes with positive excess) |
938 | 940 |
for (int u = 0; u != _res_node_num; ++u) { |
939 | 941 |
if (_excess[u] > 0) _active_nodes.push_back(u); |
940 | 942 |
} |
941 | 943 |
|
942 | 944 |
// Initialize the next arcs |
943 | 945 |
for (int u = 0; u != _res_node_num; ++u) { |
944 | 946 |
_next_out[u] = _first_out[u]; |
945 | 947 |
} |
946 | 948 |
|
947 | 949 |
// Perform partial augment and relabel operations |
948 | 950 |
while (true) { |
949 | 951 |
// Select an active node (FIFO selection) |
950 | 952 |
while (_active_nodes.size() > 0 && |
951 | 953 |
_excess[_active_nodes.front()] <= 0) { |
952 | 954 |
_active_nodes.pop_front(); |
953 | 955 |
} |
954 | 956 |
if (_active_nodes.size() == 0) break; |
955 | 957 |
int start = _active_nodes.front(); |
956 | 958 |
path_nodes.clear(); |
957 | 959 |
path_nodes.push_back(start); |
958 | 960 |
|
959 | 961 |
// Find an augmenting path from the start node |
960 | 962 |
int tip = start; |
961 | 963 |
while (_excess[tip] >= 0 && |
962 | 964 |
int(path_nodes.size()) <= max_length) { |
963 | 965 |
int u; |
964 | 966 |
LargeCost min_red_cost, rc; |
965 | 967 |
int last_out = _sum_supply < 0 ? |
966 | 968 |
_first_out[tip+1] : _first_out[tip+1] - 1; |
967 | 969 |
for (int a = _next_out[tip]; a != last_out; ++a) { |
968 | 970 |
if (_res_cap[a] > 0 && |
969 | 971 |
_cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { |
970 | 972 |
u = _target[a]; |
971 | 973 |
pred_arc[u] = a; |
972 | 974 |
_next_out[tip] = a; |
973 | 975 |
tip = u; |
974 | 976 |
path_nodes.push_back(tip); |
975 | 977 |
goto next_step; |
976 | 978 |
} |
977 | 979 |
} |
978 | 980 |
|
979 | 981 |
// Relabel tip node |
980 | 982 |
min_red_cost = std::numeric_limits<LargeCost>::max() / 2; |
981 | 983 |
for (int a = _first_out[tip]; a != last_out; ++a) { |
982 | 984 |
rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]]; |
983 | 985 |
if (_res_cap[a] > 0 && rc < min_red_cost) { |
984 | 986 |
min_red_cost = rc; |
985 | 987 |
} |
986 | 988 |
} |
987 | 989 |
_pi[tip] -= min_red_cost + _epsilon; |
988 | 990 |
|
989 | 991 |
// Reset the next arc of tip |
990 | 992 |
_next_out[tip] = _first_out[tip]; |
991 | 993 |
|
992 | 994 |
// Step back |
993 | 995 |
if (tip != start) { |
994 | 996 |
path_nodes.pop_back(); |
995 | 997 |
tip = path_nodes.back(); |
996 | 998 |
} |
997 | 999 |
|
998 | 1000 |
next_step: ; |
999 | 1001 |
} |
1000 | 1002 |
|
1001 | 1003 |
// Augment along the found path (as much flow as possible) |
1002 | 1004 |
Value delta; |
1003 | 1005 |
int u, v = path_nodes.front(), pa; |
1004 | 1006 |
for (int i = 1; i < int(path_nodes.size()); ++i) { |
1005 | 1007 |
u = v; |
1006 | 1008 |
v = path_nodes[i]; |
1007 | 1009 |
pa = pred_arc[v]; |
1008 | 1010 |
delta = std::min(_res_cap[pa], _excess[u]); |
1009 | 1011 |
_res_cap[pa] -= delta; |
1010 | 1012 |
_res_cap[_reverse[pa]] += delta; |
1011 | 1013 |
_excess[u] -= delta; |
1012 | 1014 |
_excess[v] += delta; |
1013 | 1015 |
if (_excess[v] > 0 && _excess[v] <= delta) |
1014 | 1016 |
_active_nodes.push_back(v); |
1015 | 1017 |
} |
1016 | 1018 |
} |
1017 | 1019 |
} |
1018 | 1020 |
} |
1019 | 1021 |
|
1020 | 1022 |
/// Execute the algorithm performing push and relabel operations |
1021 | 1023 |
void startPush() { |
1022 | 1024 |
// Paramters for heuristics |
1023 | 1025 |
const int BF_HEURISTIC_EPSILON_BOUND = 1000; |
1024 | 1026 |
const int BF_HEURISTIC_BOUND_FACTOR = 3; |
1025 | 1027 |
|
1026 | 1028 |
// Perform cost scaling phases |
1027 | 1029 |
BoolVector hyper(_res_node_num, false); |
1028 | 1030 |
for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? |
1029 | 1031 |
1 : _epsilon / _alpha ) |
1030 | 1032 |
{ |
1031 | 1033 |
// "Early Termination" heuristic: use Bellman-Ford algorithm |
1032 | 1034 |
// to check if the current flow is optimal |
1033 | 1035 |
if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { |
1034 | 1036 |
_arc_vec.clear(); |
1035 | 1037 |
_cost_vec.clear(); |
1036 | 1038 |
for (int j = 0; j != _res_arc_num; ++j) { |
1037 | 1039 |
if (_res_cap[j] > 0) { |
1038 | 1040 |
_arc_vec.push_back(IntPair(_source[j], _target[j])); |
1039 | 1041 |
_cost_vec.push_back(_cost[j] + 1); |
1040 | 1042 |
} |
1041 | 1043 |
} |
1042 | 1044 |
_sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); |
1043 | 1045 |
|
1044 | 1046 |
BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map); |
1045 | 1047 |
bf.init(0); |
1046 | 1048 |
bool done = false; |
1047 | 1049 |
int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num)); |
1048 | 1050 |
for (int i = 0; i < K && !done; ++i) |
1049 | 1051 |
done = bf.processNextWeakRound(); |
1050 | 1052 |
if (done) break; |
1051 | 1053 |
} |
1052 | 1054 |
|
1053 | 1055 |
// Saturate arcs not satisfying the optimality condition |
1054 | 1056 |
for (int a = 0; a != _res_arc_num; ++a) { |
1055 | 1057 |
if (_res_cap[a] > 0 && |
1056 | 1058 |
_cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { |
1057 | 1059 |
Value delta = _res_cap[a]; |
1058 | 1060 |
_excess[_source[a]] -= delta; |
1059 | 1061 |
_excess[_target[a]] += delta; |
1060 | 1062 |
_res_cap[a] = 0; |
1061 | 1063 |
_res_cap[_reverse[a]] += delta; |
1062 | 1064 |
} |
1063 | 1065 |
} |
1064 | 1066 |
|
1065 | 1067 |
// Find active nodes (i.e. nodes with positive excess) |
1066 | 1068 |
for (int u = 0; u != _res_node_num; ++u) { |
1067 | 1069 |
if (_excess[u] > 0) _active_nodes.push_back(u); |
1068 | 1070 |
} |
1069 | 1071 |
|
1070 | 1072 |
// Initialize the next arcs |
1071 | 1073 |
for (int u = 0; u != _res_node_num; ++u) { |
1072 | 1074 |
_next_out[u] = _first_out[u]; |
1073 | 1075 |
} |
1074 | 1076 |
|
1075 | 1077 |
// Perform push and relabel operations |
1076 | 1078 |
while (_active_nodes.size() > 0) { |
1077 | 1079 |
LargeCost min_red_cost, rc; |
1078 | 1080 |
Value delta; |
1079 | 1081 |
int n, t, a, last_out = _res_arc_num; |
1080 | 1082 |
|
1081 | 1083 |
// Select an active node (FIFO selection) |
1082 | 1084 |
next_node: |
1083 | 1085 |
n = _active_nodes.front(); |
1084 | 1086 |
last_out = _sum_supply < 0 ? |
1085 | 1087 |
_first_out[n+1] : _first_out[n+1] - 1; |
1086 | 1088 |
|
1087 | 1089 |
// Perform push operations if there are admissible arcs |
1088 | 1090 |
if (_excess[n] > 0) { |
1089 | 1091 |
for (a = _next_out[n]; a != last_out; ++a) { |
1090 | 1092 |
if (_res_cap[a] > 0 && |
1091 | 1093 |
_cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { |
1092 | 1094 |
delta = std::min(_res_cap[a], _excess[n]); |
1093 | 1095 |
t = _target[a]; |
1094 | 1096 |
|
1095 | 1097 |
// Push-look-ahead heuristic |
1096 | 1098 |
Value ahead = -_excess[t]; |
1097 | 1099 |
int last_out_t = _sum_supply < 0 ? |
1098 | 1100 |
_first_out[t+1] : _first_out[t+1] - 1; |
1099 | 1101 |
for (int ta = _next_out[t]; ta != last_out_t; ++ta) { |
1100 | 1102 |
if (_res_cap[ta] > 0 && |
1101 | 1103 |
_cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0) |
1102 | 1104 |
ahead += _res_cap[ta]; |
1103 | 1105 |
if (ahead >= delta) break; |
1104 | 1106 |
} |
1105 | 1107 |
if (ahead < 0) ahead = 0; |
1106 | 1108 |
|
1107 | 1109 |
// Push flow along the arc |
1108 | 1110 |
if (ahead < delta) { |
1109 | 1111 |
_res_cap[a] -= ahead; |
1110 | 1112 |
_res_cap[_reverse[a]] += ahead; |
1111 | 1113 |
_excess[n] -= ahead; |
1112 | 1114 |
_excess[t] += ahead; |
1113 | 1115 |
_active_nodes.push_front(t); |
1114 | 1116 |
hyper[t] = true; |
1115 | 1117 |
_next_out[n] = a; |
1116 | 1118 |
goto next_node; |
1117 | 1119 |
} else { |
1118 | 1120 |
_res_cap[a] -= delta; |
1119 | 1121 |
_res_cap[_reverse[a]] += delta; |
1120 | 1122 |
_excess[n] -= delta; |
1121 | 1123 |
_excess[t] += delta; |
1122 | 1124 |
if (_excess[t] > 0 && _excess[t] <= delta) |
1123 | 1125 |
_active_nodes.push_back(t); |
1124 | 1126 |
} |
1125 | 1127 |
|
1126 | 1128 |
if (_excess[n] == 0) { |
1127 | 1129 |
_next_out[n] = a; |
1128 | 1130 |
goto remove_nodes; |
1129 | 1131 |
} |
1130 | 1132 |
} |
1131 | 1133 |
} |
1132 | 1134 |
_next_out[n] = a; |
1133 | 1135 |
} |
1134 | 1136 |
|
1135 | 1137 |
// Relabel the node if it is still active (or hyper) |
1136 | 1138 |
if (_excess[n] > 0 || hyper[n]) { |
1137 | 1139 |
min_red_cost = std::numeric_limits<LargeCost>::max() / 2; |
1138 | 1140 |
for (int a = _first_out[n]; a != last_out; ++a) { |
1139 | 1141 |
rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]]; |
1140 | 1142 |
if (_res_cap[a] > 0 && rc < min_red_cost) { |
1141 | 1143 |
min_red_cost = rc; |
1142 | 1144 |
} |
1143 | 1145 |
} |
1144 | 1146 |
_pi[n] -= min_red_cost + _epsilon; |
1145 | 1147 |
hyper[n] = false; |
1146 | 1148 |
|
1147 | 1149 |
// Reset the next arc |
1148 | 1150 |
_next_out[n] = _first_out[n]; |
1149 | 1151 |
} |
1150 | 1152 |
|
1151 | 1153 |
// Remove nodes that are not active nor hyper |
1152 | 1154 |
remove_nodes: |
1153 | 1155 |
while ( _active_nodes.size() > 0 && |
1154 | 1156 |
_excess[_active_nodes.front()] <= 0 && |
1155 | 1157 |
!hyper[_active_nodes.front()] ) { |
1156 | 1158 |
_active_nodes.pop_front(); |
1157 | 1159 |
} |
1158 | 1160 |
} |
1159 | 1161 |
} |
1160 | 1162 |
} |
1161 | 1163 |
|
1162 | 1164 |
}; //class CostScaling |
1163 | 1165 |
|
1164 | 1166 |
///@} |
1165 | 1167 |
|
1166 | 1168 |
} //namespace lemon |
1167 | 1169 |
|
1168 | 1170 |
#endif //LEMON_COST_SCALING_H |
0 comments (0 inline)